JSirene.cc File Reference

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <limits>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TProfile.h"
#include "TRandom3.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "JLang/JSharedPointer.hh"
#include "JLang/Jpp.hh"
#include "JLang/JPredicate.hh"
#include "JSystem/JDate.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JCDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JRadiationSource.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JAAnet/JHit_t.hh"
#include "JAAnet/JTrk_t.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JPDB.hh"
#include "JSirene/km3-opa_efrac.hh"
#include "JSirene/JSeaWater.hh"
#include "JSirene/JCDFTable1D.hh"
#include "JSirene/JCDFTable2D.hh"
#include "JSirene/JSireneToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSubset.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Variables
int	debug
	debug level More...

Detailed Description

Main program to simulate detector response to muons and showers.

Picture by Claudine Colnard

Note that CDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD;
The file names are obtained by replacing JPHYSICS::WILD_CARD; with

• JPHYSICS::DIRECT_LIGHT_FROM_MUON;
• JPHYSICS::SCATTERED_LIGHT_FROM_MUON;
• JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS;
• JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS;
• JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER; and
• JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER, respectively.

The CDF tables can be produced with the script JMakePDF.sh:

JMakePDF.sh -P

(option -h will print all available options). Note that the script will launch a large number of processes (JMakePDF and JMakePDG)
which may take a considerable amount of time to completion.
On a standard desktop, all jobs should be finished within 1/2 a day or so.

The same script should then be run with option -M to merge the PDF files, i.e:

JMakePDF.sh -M

CDF tables are obtained by running the same script with option -C, i.e:

JMakePDF.sh -C

The various PDFs can be drawn using the JDrawPDF or JDrawPDG applications.
The tabulated PDFs can be plotted using the JPlotPDF or JPlotPDG applications.
The tabulated CDFs can be plotted using the JPlotCDF or JPlotCDG applications.

Author

mdejong

Definition in file JSirene.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 171 of file JSirene.cc.

{
  using namespace std;
  using namespace JPP;
 
  string         fileDescriptor;
  JMultipleFileScanner<JAAnetTypes_t>                                           inputFile;
  JFileRecorder       <JTYPELIST<JAAnetTypes_t, JMeta, JRootTypes_t>::typelist> outputFile;
  JLimit_t&      numberOfEvents = inputFile.getLimit();
  string         detectorFile;
  double         Ecut_GeV  =  0.1;                            // minimal energy for generation of light from shower [GeV]
  double         Emin_GeV  =  1.0;                            // minimal energy of muon for shower generation [GeV]
  double         Dmin_m    =  0.1;                            // minimal distance for positioning [m]
  double         Tmax_ns   =  0.1;                            // maximal time between hits on same PMT to be merged
  int            Nmax_npe  =  numeric_limits<int>::max();     // maximal number of photo-electrons
  bool           writeEMShowers;
  bool           keep;
  size_t         numberOfHits;
  double         factor;
  UInt_t         seed;
 
  try { 
 
    JProperties properties;
 
    properties.insert(gmake_property(Ecut_GeV));
    properties.insert(gmake_property(Emin_GeV));
    properties.insert(gmake_property(Dmin_m));
    properties.insert(gmake_property(Tmax_ns));
    properties.insert(gmake_property(Nmax_npe));
 
    JParser<> zap("Main program to simulate detector response to muons and showers.");
 
    zap['@'] = make_field(properties)        = JPARSER::initialised();
    zap['F'] = make_field(fileDescriptor,    "file name descriptor for CDF tables");
    zap['f'] = make_field(inputFile)         = JPARSER::initialised();
    zap['o'] = make_field(outputFile)        = "sirene.root";
    zap['n'] = make_field(numberOfEvents)    = JLimit::max();
    zap['a'] = make_field(detectorFile)      = "";
    zap['s'] = make_field(writeEMShowers,    "store generated EM showers in event");
    zap['k'] = make_field(keep,              "keep position of tracks");
    zap['N'] = make_field(numberOfHits,      "minimum number of hits to output event") = 1;
    zap['U'] = make_field(factor,            "scaling factor applied to light yields") = 1.0;
    zap['S'] = make_field(seed)              = 0;
    zap['d'] = make_field(debug)             = 1;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  gRandom->SetSeed(seed);
 
 
  const JMeta meta(argc, argv);
 
  const double Zbed = 0.0;                             // level of seabed  [m]
 
  vector<JCDF_t> CDF;
  vector<JCDG_t> CDG;
 
  CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_MUON));
  CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_MUON));
  CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_DELTARAYS));
  CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_DELTARAYS));
  
  CDG.push_back(JCDG_t(fileDescriptor, DIRECT_LIGHT_FROM_EMSHOWER));
  CDG.push_back(JCDG_t(fileDescriptor, SCATTERED_LIGHT_FROM_EMSHOWER));
 
  
  double maximal_road_width = 0.0;                     // road width [m]
  double maximal_distance   = 0.0;                     // road width [m]
 
  for (size_t i = 0; i != CDF.size(); ++i) {
    
    DEBUG("Range CDF["<< CDF[i].type << "] " << CDF[i].function.intensity.getXmax() << " m" << endl);
 
    maximal_road_width = max(maximal_road_width, CDF[i].function.intensity.getXmax());
  }
 
  for (size_t i = 0; i != CDG.size(); ++i) {
 
    DEBUG("Range CDG["<< CDG[i].type << "] " << CDG[i].function.intensity.getXmax() << " m" << endl);
 
    if (!is_scattered(CDF[i].type)) {
      maximal_road_width = max(maximal_road_width, CDG[i].function.intensity.getXmax());
    }
    
    maximal_distance   = max(maximal_distance,   CDG[i].function.intensity.getXmax());
  }
  
  NOTICE("Maximal road width [m] " << maximal_road_width << endl);
  NOTICE("Maximal distance   [m] " << maximal_distance   << endl);
 
 
  if (detectorFile == "" || inputFile.empty()) {
    STATUS("Nothing to be done." << endl);
    return 0;
  }
 
  JDetector detector;
 
  try {
 
    STATUS("Load detector... " << flush);
 
    load(detectorFile, detector);
 
    STATUS("OK" << endl);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
 
  vector< JSharedPointer<JRadiationInterface> > radiation;
 
  if (true) {
 
    STATUS("Setting up radiation tables... " << flush);
 
    const JRadiation hydrogen( 1.0,  1.0, 40, 0.01, 0.1, 0.1);
    const JRadiation oxygen  ( 8.0, 16.0, 40, 0.01, 0.1, 0.1);
    const JRadiation chlorine(17.0, 35.0, 40, 0.01, 0.1, 0.1);
 
    JSharedPointer<JRadiation> Hydrogen(new JRadiationFunction(hydrogen, 300, 0.2, 1.0e11));
    JSharedPointer<JRadiation> Oxygen  (new JRadiationFunction(oxygen,   300, 0.2, 1.0e11));
    JSharedPointer<JRadiation> Chlorine(new JRadiationFunction(chlorine, 300, 0.2, 1.0e11));
    
    JRadiationSource::source_type EErad(&JRadiation::TotalCrossSectionEErad, &JRadiation::EfromEErad);
    JRadiationSource::source_type Brems(&JRadiation::TotalCrossSectionBrems, &JRadiation::EfromBrems);
    JRadiationSource::source_type GNrad(&JRadiation::TotalCrossSectionGNrad, &JRadiation::EfromGNrad);
 
    radiation.push_back(new JRadiationSource(Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  EErad));
    radiation.push_back(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad));
    radiation.push_back(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  EErad));
 
    radiation.push_back(new JRadiationSource(Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  Brems));
    radiation.push_back(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems));
    radiation.push_back(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  Brems));
 
    radiation.push_back(new JRadiationSource(Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  GNrad));
    radiation.push_back(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad));
    radiation.push_back(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  GNrad));
 
    STATUS("OK" << endl);
  }
 
 
  Vec  center(0,0,0);
  Head header;
 
  try {
 
    header = inputFile.getHeader();
 
    JHead buffer(header);
 
    center = get<Vec>(buffer);
 
    buffer.simul.push_back(JAANET::simul());
 
    buffer.simul.rbegin()->program = "JSirene";
    buffer.simul.rbegin()->version = getGITVersion();
    buffer.simul.rbegin()->date    = getDate();
    buffer.simul.rbegin()->time    = getTime();
 
    buffer.push(&JHead::simul);
 
    if (!keep) {
 
      buffer.coord_origin  = coord_origin(0.0, 0.0, 0.0);
      buffer.can.zmin     += center.z;
      buffer.can.zmax     += center.z;
 
      buffer.push(&JHead::coord_origin);
      buffer.push(&JHead::can);
    }
  
    const JCircle2D circle(detector.begin(), detector.end());
 
    center += Vec(circle.getX(), circle.getY(), 0.0);
 
    copy(buffer, header);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  if (!keep)
    NOTICE("Offset applied to true tracks is: " << center << endl);
  else
    NOTICE("No offset applied to true tracks." << endl);
 
  JCylinder3D cylinder(detector.begin(), detector.end());
 
  cylinder.addMargin(maximal_distance);
 
  if (cylinder.getZmin() < Zbed) {
    cylinder.setZmin(Zbed);
  }
 
  NOTICE("Ligth generation volume: " << cylinder << endl);
 
  TProfile cpu("cpu", NULL,  14,  1.0,   8.0);
  TH1D     job("job", NULL, 400,  0.5, 400.5);
 
 
  outputFile.open();
 
  if (!outputFile.is_open()) {
    FATAL("Error opening file " << outputFile << endl);
  }
 
  outputFile.put(meta);
  outputFile.put(header);
  outputFile.put(*gRandom);  
 
  
  JTimer timer;
 
  for (JMultipleFileScanner<Evt>& in = inputFile; in.hasNext(); ) {
 
    STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
 
    job.Fill(1.0);
 
    Evt* evt = in.next(); 
 
    if (!keep) {
      for (vector<Trk>::iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
        track->pos += center;
      }
    }
    
    Evt event(*evt);                                   // output
 
    event.mc_hits.clear();
 
    JHits_t mc_hits;                                   // temporary buffer
 
    timer.reset();
    timer.start();
    
    for (vector<Trk>::const_iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
 
      if (is_muon(*track)) {
 
        // -----------------------------------------------
        //                      muon
        // -----------------------------------------------
 
        job.Fill(2.0);
 
        typedef JDetectorSubset<JModuleComparator3Z_t>   JDetectorSubset_t;
 
        const pair<double, double> intersection = cylinder.getIntersection(getAxis(*track));
 
        double Zmin = intersection.first;
        double Zmax = intersection.second;
 
        if (Zmax - Zmin <= Dmin_m) {
          continue;
        }
        
        JVertex vertex(0.0, track->t, track->E);                        // start of muon
      
        if (track->pos.z < Zbed) {                                      // propagate muon through rock
            
          if (track->dir.z > 0.0)
            vertex.step(gRock, (Zbed - track->pos.z) / track->dir.z); 
          else
            continue;
        }
 
        if (vertex.getZ() < Zmin) {                                     // propagate muon through water
          vertex.step(gWater, Zmin - vertex.getZ());
        }
 
        if (vertex.getRange() <= Dmin_m) {
          continue;
        }
 
        job.Fill(3.0);
 
        const JDetectorSubset_t subdetector(detector, getAxis(*track), maximal_road_width);
 
        if (subdetector.empty()) {
          continue;
        }
 
        job.Fill(4.0);
        
        JTrack muon(vertex);                                            // propagate muon trough detector
 
        while (vertex.getE() >= Emin_GeV && vertex.getZ() < Zmax) {
        
          const int N = radiation.size();
        
          double li[N];                                                 // inverse interaction lengths
          double ls = 1.0e-5;                                           // minimal total inverse interaction length (100 km)^-1
        
          for (int i = 0; i != N; ++i) {
            ls += li[i] = radiation[i]->getInverseInteractionLength(vertex.getE());
          }
 
          const double ds = min(gRandom->Exp(1.0) / ls, vertex.getRange());
 
          vertex.step(ds);
 
          if (vertex.getE() >= Emin_GeV) {
 
            double Es = 0.0;                                            // shower energy [GeV]
            double y  = gRandom->Uniform(ls);
            
            for (int i = 0; i != N; ++i) {
              
              y -= li[i];
              
              if (y < 0.0) {
                Es = radiation[i]->getEnergyOfShower(vertex.getE());    // shower energy [GeV]
                break;
              }
            }
            
            vertex.applyEloss(Es);
 
            if (Es >= Ecut_GeV) {
              muon.push_back(vertex);
            }
          }
        }
 
        // add muon end point
        
        if (vertex.getE() < Emin_GeV && vertex.getRange() > 0.0) {
 
          vertex.step(vertex.getRange());
 
          muon.push_back(vertex);
        }
 
        // add information to output muon
        
        vector<Trk>::iterator trk = find_if(event.mc_trks.begin(),
                                            event.mc_trks.end(),
                                            make_predicate(&Trk::id, track->id));
 
        if (trk != event.mc_trks.end()) {
          trk->len = (muon.rbegin()->getZ() < Zmax ? +1 : -1) * (muon.rbegin()->getZ() - muon.begin()->getZ());
          trk->setusr(mc_usr_keys::energy_lost_in_can, muon.begin()->getE() - muon.rbegin()->getE());
        }
        
        for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
 
          const double z0 = muon.begin()->getZ();
          const double t0 = muon.begin()->getT();
          const double R  = module->getX();
          const double Z  = module->getZ() - R / getTanThetaC();
 
          if (Z >= muon.begin()->getZ() && Z <= muon.rbegin()->getZ()) {
          
            for (size_t i = 0; i != CDF.size(); ++i) {
                      
              if (R < CDF[i].integral.getXmax()) {
              
                try {
 
                  double W = 1.0;                                       // mip
                  
                  if (is_deltarays(CDF[i].type)) {
                    W = getDeltaRaysFromMuon(muon.getE(Z));             // delta-rays
                  }
 
                  const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
                  const int    N   = gRandom->Poisson(NPE); 
 
                  if (N != 0) {
                    
                    double ns = 0.0;
                  
                    for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
 
                      const double R     = pmt->getX();
                      const double Z     = pmt->getZ() - z0;
                      const double theta = pmt->getTheta();
                      const double phi   = fabs(pmt->getPhi());
                    
                      const double npe   = CDF[i].function.getNPE(R, theta, phi) * factor * W;
                    
                      ns += npe;
                    
                      int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), Nmax_npe);
                    
                      job.Fill((double) (100 + CDF[i].type), (double) n0);
                    
                      while (n0 != 0) {
                      
                        const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
                        const int    n1 = getNumberOfPhotoElectrons(n0);
                     
                        mc_hits.push_back(JHit_t(mc_hits.size() + 1,
                                                 pmt->getID(),
                                                 getHitType(CDF[i].type),
                                                 track->id,
                                                 t0  +  (R * getTanThetaC() + Z) / C  +  t1,
                                                 n1));
 
                        n0 -= n1;
                      }
                    }
 
                    if (ns > NPE) {
                      job.Fill((double) (300 + CDF[i].type));
                    }
                  }
                }
                catch(const exception& error) {
                  job.Fill((double) (200 + CDF[i].type));
                }
              }
            }
          }
        }
 
        for (JTrack::const_iterator vertex = muon.begin(); vertex != muon.end(); ++vertex) {
 
          const double z0 = vertex->getZ();
          const double t0 = vertex->getT();
          const double Es = vertex->getEs();
          
          int origin = track->id;
 
          if (writeEMShowers) {
            origin = event.mc_trks.size() + 1;
          }
 
          int number_of_hits = 0;
            
          JDetectorSubset_t::range_type range = subdetector.getRange(z0 - maximal_distance,
                                                                     z0 + maximal_distance);
 
          for (JDetector::const_iterator module = range.begin(); module != range.end(); ++module) {
 
            const double dz = geanz.getMaximum(Es);
            const double R  = module->getX();
            const double Z  = module->getZ() - z0 - dz;
            const double D  = sqrt(R*R + Z*Z);
            const double cd = Z / D; 
            
            for (size_t i = 0; i != CDG.size(); ++i) {
              
              if (D < CDG[i].integral.getXmax()) {
                
                try {
 
                  const double NPE = CDG[i].integral.getNPE(D, cd) * Es * module->size() * factor;
                  const int    N   = gRandom->Poisson(NPE); 
 
                  if (N != 0) {
                      
                    double ns = 0.0;
 
                    for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
 
                      const double dz    = geanz.getMaximum(Es);
                      const double R     = pmt->getX();
                      const double Z     = pmt->getZ() - z0 - dz;
                      const double theta = pmt->getTheta();
                      const double phi   = fabs(pmt->getPhi());
                      const double D     = sqrt(R*R + Z*Z);
                      const double cd    = Z / D; 
                        
                      const double npe   = CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor;
                      
                      ns += npe;
                        
                      int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), Nmax_npe);
 
                      job.Fill((double) (100 + CDG[i].type), (double) n0);
                      
                      while (n0 != 0) {
 
                        const double dz = geanz.getLength(Es, gRandom->Rndm());
                        const double Z  = pmt->getZ() - z0 - dz;
                        const double D  = sqrt(R*R + Z*Z);
                        const double cd = Z / D; 
                          
                        const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
                        const int    n1 = getNumberOfPhotoElectrons(n0);
 
                        mc_hits.push_back(JHit_t(mc_hits.size() + 1,
                                                 pmt->getID(),
                                                 getHitType(CDG[i].type),
                                                 origin,
                                                 t0  +  (dz +  D * getIndexOfRefraction()) / C  +  t1,
                                                 n1));
 
                        n0 -= n1;
                        
                        number_of_hits += n1;
                      }
                    }
                    
                    if (ns > NPE) {
                      job.Fill((double) (300 + CDG[i].type));
                    }
                  }
                }
                catch(const exception& error) {
                  job.Fill((double) (200 + CDG[i].type));
                }
              }
            }
          }
          
          if (writeEMShowers && number_of_hits != 0) {
 
            event.mc_trks.push_back(JTrk_t(origin,
                                           TRACK_TYPE_PHOTON,
                                           track->id,
                                           track->pos + track->dir * vertex->getZ(),
                                           track->dir,
                                           vertex->getT(),
                                           Es));
          }
        }
 
      } else if (is_tau(*track)) {
 
        // -----------------------------------------------
        //                       tau
        // -----------------------------------------------
 
        job.Fill(6.0);
 
        const double z0 = 0.0;
        const double z1 = z0 + track->len;
        const double t0 = track->t;
        const double E  = track->E;
 
        const JTransformation3D transformation = getTransformation(*track);
        
        JModule buffer;
        
        for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
 
          const JPosition3D pos = transformation.transform(module->getPosition());
 
          const double R  = pos.getX();
          const double Z  = pos.getZ() - R / getTanThetaC();
 
          if (Z < z0  ||
              Z > z1  ||
              R > maximal_road_width) {
            continue;
          }
          
          for (size_t i = 0; i != CDF.size(); ++i) {
                
            if (R < CDF[i].integral.getXmax()) {
              
              try {
 
                double W = 1.0;                                       // mip
                    
                if (is_deltarays(CDF[i].type)) {
                  W = getDeltaRaysFromTau(E);                         // delta-rays
                }
                
                const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
                const int    N   = gRandom->Poisson(NPE); 
                
                if (N != 0) {
                    
                  buffer = *module;
            
                  buffer.transform(transformation);
 
                  double ns = 0.0;
                  
                  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
 
                    const double R     = pmt->getX();
                    const double Z     = pmt->getZ() - z0;
                    const double theta = pmt->getTheta();
                    const double phi   = fabs(pmt->getPhi());
                    
                    const double npe   = CDF[i].function.getNPE(R, theta, phi) * factor * W;
                    
                    ns += npe;
                    
                    int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), Nmax_npe);
                    
                    job.Fill((double) (120 + CDF[i].type), (double) n0);
                    
                    while (n0 != 0) {
                      
                      const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
                      const int    n1 = getNumberOfPhotoElectrons(n0);
                      
                      mc_hits.push_back(JHit_t(mc_hits.size() + 1,
                                               pmt->getID(),
                                               getHitType(CDF[i].type),
                                               track->id,
                                               t0  +  (R * getTanThetaC() + Z) / C  +  t1,
                                               n1));
                      
                      n0 -= n1;
                    }
                  }
 
                  if (ns > NPE) {
                    job.Fill((double) (320 + CDF[i].type));
                  }
                }
              }
              catch(const exception& error) {
                job.Fill((double) (220 + CDF[i].type));
              }
 
            }
          }
        }
        
        if (!buffer.empty()) {
          job.Fill(7.0);
        }
          
      } else if (!is_neutrino(*track)) {
 
        if (JPDB::getInstance().hasPDG(track->type)) {
    
          // -----------------------------------------------
          //               electron or hadron
          // -----------------------------------------------
          
          job.Fill(8.0);
        
          double E = track->E;
 
          try {
            E = getKineticEnergy(E, JPDB::getInstance().getPDG(track->type).mass);
          }
          catch(const exception& error) {
            ERROR(error.what() << endl);
          }
 
          E = pythia(track->type, E);
              
          if (E < Ecut_GeV || cylinder.getDistance(getPosition(*track)) > Dmin_m) {
            continue;
          }
        
          const double z0 = 0.0;
          const double t0 = track->t;
 
          const JTransformation3D transformation = getTransformation(*track);
        
          JModule buffer;
        
          for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
 
            const JPosition3D pos = transformation.transform(module->getPosition());
 
            const double dz = geanz.getMaximum(E);
            const double R  = pos.getX();
            const double Z  = pos.getZ() - z0 - dz;
            const double D  = sqrt(R*R + Z*Z);
            const double cd = Z / D; 
            
            if (D > maximal_distance) {
              continue;
            }
            
            for (size_t i = 0; i != CDG.size(); ++i) {
                    
              if (D < CDG[i].integral.getXmax()) {
                
                try {
                        
                  const double NPE = CDG[i].integral.getNPE(D, cd) * E * module->size() * factor;
                  const int    N   = gRandom->Poisson(NPE); 
 
                  if (N != 0) {
 
                    buffer = *module;
            
                    buffer.transform(transformation);
          
                    double ns = 0.0;
                        
                    for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
                    
                      const double dz    = geanz.getMaximum(E);
                      const double R     = pmt->getX();
                      const double Z     = pmt->getZ() - z0 - dz;
                      const double theta = pmt->getTheta();
                      const double phi   = fabs(pmt->getPhi());
                      const double D     = sqrt(R*R + Z*Z);
                      const double cd    = Z / D; 
                          
                      const double npe   = CDG[i].function.getNPE(D, cd, theta, phi) * E * factor;
 
                      ns += npe;
                          
                      int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), Nmax_npe);
                 
                      job.Fill((double) (140 + CDG[i].type), (double) n0);
                          
                      while (n0 != 0) {
                      
                        const double dz = geanz.getLength(E, gRandom->Rndm());
                        const double Z  = pmt->getZ() - z0 - dz;
                        const double D  = sqrt(R*R + Z*Z);
                        const double cd = Z / D;
                      
                        const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
                        const int    n1 = getNumberOfPhotoElectrons(n0);
 
                        mc_hits.push_back(JHit_t(mc_hits.size() + 1,
                                                 pmt->getID(),
                                                 getHitType(CDG[i].type, true),
                                                 track->id,
                                                 t0  +  (dz +  D * getIndexOfRefraction()) / C  +  t1,
                                                 n1));
 
                        n0 -= n1;
                      }
                    }
                  
                    if (ns > NPE) {
                      job.Fill((double) (340 + CDG[i].type));
                    }
                  }
                }
                catch(const exception& error) {
                  job.Fill((double) (240 + CDG[i].type));
                }
              }
            }
          }
 
          if (!buffer.empty()) {
            job.Fill(9.0);
          }
 
        } else {
 
          job.Fill(21.0);
        }
      }
    }
 
    if (!mc_hits.empty()) {
 
      mc_hits.merge(Tmax_ns);
 
      event.mc_hits.resize(mc_hits.size());
 
      copy(mc_hits.begin(), mc_hits.end(), event.mc_hits.begin());
    }
 
    timer.stop();
 
    if (has_neutrino(event)) {
      cpu.Fill(log10(get_neutrino(event).E), (double) timer.usec_ucpu * 1.0e-3);
    }
 
    if (event.mc_hits.size() >= numberOfHits) {
 
      outputFile.put(event);
 
      job.Fill(10.0);
    }
  }
  STATUS(endl); 
 
  outputFile.put(cpu);
  outputFile.put(job);
  outputFile.put(*gRandom);
  outputFile.close();
}

Variable Documentation

debug

int debug

debug level

Definition at line 59 of file JSirene.cc.