JSirene.cc

Go to the documentation of this file.

#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 "evt/Head.hh"
#include "evt/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/pythia.hh"
//#include "JSirene/JPythia.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"


int debug;                              //!< debug level

namespace {

  using namespace JPP;

  typedef JHermiteSplineFunction1D_t                                             JFunction1D_t;
  typedef JMAPLIST<JPolint1FunctionalMap,
                   JPolint1FunctionalGridMap,
                   JPolint1FunctionalGridMap>::maplist                           J3DMap_t;
  typedef JMAPLIST<JPolint1FunctionalMap,
                   JPolint1FunctionalMap,
                   JPolint1FunctionalGridMap,
                   JPolint1FunctionalGridMap>::maplist                           J4DMap_t;

  typedef JCDFTable<JFunction1D_t, J3DMap_t>                                     JCDF4D_t;  // muon
  typedef JCDFTable<JFunction1D_t, J4DMap_t>                                     JCDF5D_t;  // shower

  typedef JCDFTable1D<JFunction1D_t::argument_type, JFunction1D_t::result_type>  JCDF1D_t;  // muon
  typedef JCDFTable2D<JFunction1D_t::argument_type, JFunction1D_t::result_type>  JCDF2D_t;  // shower

  
  const int    NUMBER_OF_BINS = 200;    //!< number of bins for CDF integral
  const double SAFETY_FACTOR  = 1.7;    //!< safety factor  for CDF integral

  
  /**
   * Auxiliary class for CDF tables.
   */
  template<class function_type,         // CDF function
           class integral_type>         // CDF integral
  struct JABC {
    /**
     * Constructor.
     *
     * \param  file_descriptor     file name descriptor
     * \param  type                PDF type
     */
    JABC(const std::string& file_descriptor, const JPDFType_t type) 
    {
      this->type = type;

      const std::string file_name = getFilename(file_descriptor, this->type);

      STATUS("loading input from file " << file_name << "... " << flush);
      
      try {
        function.load(file_name.c_str());
      }
      catch(const JException& error) {
        FATAL(error.what() << endl);
      }

      integral = integral_type(function, NUMBER_OF_BINS, SAFETY_FACTOR);

      STATUS("OK" << endl);
    }

    JPDFType_t    type;
    function_type function;
    integral_type integral;
  };

  
  typedef JABC<JCDF4D_t, JCDF1D_t>  JCDF_t;    //!< muon   light
  typedef JABC<JCDF5D_t, JCDF2D_t>  JCDG_t;    //!< shower light
}


/**
 * \file
 *
 * Main program to simulate detector response to muons and showers.
 *
 * 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 <tt>JMakePDF.sh</tt>:
 * <pre>
 * JMakePDF.sh -P
 * </pre>
 * (option <tt>-h</tt> will print all available options).
 * Note that the script will launch a large number of processes (<tt>JMakePDF</tt> and <tt>JMakePDG</tt>)
 * 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 <tt>-M</tt> to merge the PDF files, i.e:
 * <pre>
 * JMakePDF.sh -M
 * </pre>
 *
 * CDF tables are obtained by running the same script with option <tt>-C</tt>, i.e:
 * <pre>
 * JMakePDF.sh -C
 * </pre>
 *
 * The various PDFs can be drawn using the <tt>JDrawPDF</tt> or <tt>JDrawPDG</tt> applications.
 * The tabulated PDFs can be plotted using the <tt>JPlotPDF</tt> or <tt>JPlotPDG</tt> applications. 
 * The tabulated CDFs can be plotted using the <tt>JPlotCDF</tt> or <tt>JPlotCDG</tt> applications. 
 * \author mdejong
 */
int main(int argc, char **argv)
{
  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;
  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['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);
    
    buffer.simul.program = "JSirene";
    buffer.simul.version = getGITVersion();
    buffer.simul.date    = getDate();
    buffer.simul.time    = getTime();
    
    center = get<Vec>(buffer);

    buffer.coord_origin  = coord_origin(0.0, 0.0, 0.0);
    buffer.can.zmin     += center.z;
    buffer.can.zmax     += center.z;
      
    const JCircle2D circle(detector.begin(), detector.end());

    center += Vec(circle.getX(), circle.getY(), 0.0);

    push(buffer, &JHead::simul);
    push(buffer, &JHead::coord_origin);
    push(buffer, &JHead::can);
    copy(buffer, header);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  NOTICE("Offset applied to true tracks is: " << center << endl);
  
  JCylinder3D cylinder(detector.begin(), detector.end());

  cylinder.addMargin(maximal_distance);

  if (cylinder.getZmin() < Zbed) {
    cylinder.setZmin(Zbed);
  }


  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();                              // input

    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();
}