JSirene.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <limits>
#include <numeric>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "TRandom3.h"

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/MultiHead.hh"
#include "km3net-dataformat/offline/Evt.hh" 
#include "km3net-dataformat/definitions/applications.hh"
#include "km3net-dataformat/definitions/trkmembers.hh"

#include "JLang/JSharedPointer.hh"
#include "JLang/Jpp.hh"
#include "JLang/JPredicate.hh"
#include "JSystem/JDateAndTime.hh"

#include "JPhysics/JCDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JGeane.hh"
#include "JPhysics/JGeanz.hh"
#include "JPhysics/JRadiation.hh"
#include "JPhysics/JRadiationFunction.hh"
#include "JPhysics/JRadiationSource.hh"
#include "JPhysics/JACoeffSource.hh"
#include "JPhysics/JPhysicsSupportkit.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"

#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JPDB.hh"

#include "JSirene/JSirene.hh"
#include "JSirene/pythia.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
int    numberOfBins   = 200;               //!< number of bins for average CDF integral of optical module
double safetyFactor   = 1.7;               //!< safety factor  for average CDF integral of optical module


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

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

      this->type = type;

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

      new (&integral) integral_type(function, numberOfBins, safetyFactor);

      STATUS("OK" << endl);
    }

    JPDFType_t    type;
    function_type function;
    integral_type integral;
  };

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


  /**
   * Get random direction.
   *
   * \param  t2                  square of theta RMS [rad^2]
   * \return                     direction
   */
  inline JVersor3Z getRandomDirection(const double t2)
  {
    const double tv  = sqrt(gRandom->Exp(1.0) * t2);
    const double phi = gRandom->Uniform(-PI, +PI);
            
    return JVersor3Z(tv*cos(phi), tv*sin(phi));
  }
}


/**
 * \file
 *
 * Main program to simulate detector response to muons and showers.
 *
 * \image html sirene.png "Picture by Claudine Colnard"
 *
 * Note that CDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD;\n
 * 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.
 *
 * More accuracy can be achieved by setting compile option RADITION but it will run slower.
 *
 * 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>)\n
 * which may take a considerable amount of time to completion.\n
 * 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.\n
 * The tabulated PDFs can be plotted using the <tt>JPlotPDF</tt> or <tt>JPlotPDG</tt> applications.\n 
 * 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, JMetaTypes_t, JRootTypes_t>::typelist> outputFile;
  JLimit_t&      numberOfEvents = inputFile.getLimit();
  string         detectorFile;
  JParameters    parameters;
  bool           writeEMShowers;
  bool           keep;
  size_t         numberOfHits;
  double         factor;
  UInt_t         seed;

  try { 

    JProperties properties;

    properties.insert(gmake_property(parameters.Ecut_GeV));
    properties.insert(gmake_property(parameters.Emin_GeV));
    properties.insert(gmake_property(parameters.Dmin_m));
    properties.insert(gmake_property(parameters.Emax_GeV));
    properties.insert(gmake_property(parameters.Dmax_m));
    properties.insert(gmake_property(parameters.Tmax_ns));
    properties.insert(gmake_property(parameters.Nmax_NPE));
    properties.insert(gmake_property(parameters.Nmax_PMT));

    properties.insert(gmake_property(numberOfBins));
    properties.insert(gmake_property(safetyFactor));

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

  // remove empty modules

  for (JDetector::iterator module = detector.begin(); module != detector.end(); ) {
    if (!module->empty())
      ++module;
    else
      module = detector.erase(module);
  }

  vector< JSharedPointer<JRadiationInterface> > radiation;
  vector< JSharedPointer<JACoeffInterface>    > ionization;
  
  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.5, 40, 0.01, 0.1, 0.1);
    const JRadiation sodium   (11.0, 23.0, 40, 0.01, 0.1, 0.1);
#ifdef RADIATION
    const JRadiation calcium  (20.0, 40.0, 40, 0.01, 0.1, 0.1);
    const JRadiation magnesium(12.0, 24.3, 40, 0.01, 0.1, 0.1);
    const JRadiation potassium(19.0, 39.0, 40, 0.01, 0.1, 0.1);
    const JRadiation sulphur  (16.0, 32.0, 40, 0.01, 0.1, 0.1);
#endif

    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));
    JSharedPointer<JRadiation> Sodium   (new JRadiationFunction(sodium,   300, 0.2, 1.0e11));
#ifdef RADIATION
    JSharedPointer<JRadiation> Calcium  (new JRadiationFunction(calcium,  300, 0.2, 1.0e11));
    JSharedPointer<JRadiation> Magnesium(new JRadiationFunction(magnesium,300, 0.2, 1.0e11));
    JSharedPointer<JRadiation> Potassium(new JRadiationFunction(potassium,300, 0.2, 1.0e11));
    JSharedPointer<JRadiation> Sulphur  (new JRadiationFunction(sulphur,  300, 0.2, 1.0e11));
#endif

    radiation.push_back(new JRadiationSource(11, Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  EErad_t));
    radiation.push_back(new JRadiationSource(12, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad_t));
    radiation.push_back(new JRadiationSource(13, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  EErad_t));
    radiation.push_back(new JRadiationSource(14, Sodium,   DENSITY_SEA_WATER * JSeaWater::Na(), EErad_t));
#ifdef RADIATION
    radiation.push_back(new JRadiationSource(15, Calcium,  DENSITY_SEA_WATER * JSeaWater::Ca(), EErad_t));
    radiation.push_back(new JRadiationSource(16, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), EErad_t));
    radiation.push_back(new JRadiationSource(17, Potassium,DENSITY_SEA_WATER * JSeaWater::K(),  EErad_t));
    radiation.push_back(new JRadiationSource(18, Sulphur,  DENSITY_SEA_WATER * JSeaWater::S(),  EErad_t));
#endif

    radiation.push_back(new JRadiationSource(21, Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  Brems_t));
    radiation.push_back(new JRadiationSource(22, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems_t));
    radiation.push_back(new JRadiationSource(23, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  Brems_t)); 
    radiation.push_back(new JRadiationSource(24, Sodium,   DENSITY_SEA_WATER * JSeaWater::Na(), Brems_t)); 
#ifdef RADIATION
    radiation.push_back(new JRadiationSource(25, Calcium,  DENSITY_SEA_WATER * JSeaWater::Ca(), Brems_t));
    radiation.push_back(new JRadiationSource(26, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), Brems_t));
    radiation.push_back(new JRadiationSource(27, Potassium,DENSITY_SEA_WATER * JSeaWater::K(),  Brems_t));
    radiation.push_back(new JRadiationSource(28, Sulphur,  DENSITY_SEA_WATER * JSeaWater::S(),  Brems_t));
#endif

    radiation.push_back(new JRadiationSource(31, Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  GNrad_t)); 
    radiation.push_back(new JRadiationSource(32, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad_t)); 
    radiation.push_back(new JRadiationSource(33, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  GNrad_t)); 
    radiation.push_back(new JRadiationSource(34, Sodium,   DENSITY_SEA_WATER * JSeaWater::Na(), GNrad_t)); 
#ifdef RADIATION
    radiation.push_back(new JRadiationSource(35, Calcium,  DENSITY_SEA_WATER * JSeaWater::Ca(), GNrad_t));
    radiation.push_back(new JRadiationSource(36, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), GNrad_t));
    radiation.push_back(new JRadiationSource(37, Potassium,DENSITY_SEA_WATER * JSeaWater::K(),  GNrad_t));
    radiation.push_back(new JRadiationSource(38, Sulphur,  DENSITY_SEA_WATER * JSeaWater::S(),  GNrad_t));
#endif

    radiation.push_back(new JDISSource(100, DENSITY_SEA_WATER));

    ionization.push_back(new JACoeffSource(Oxygen,   DENSITY_SEA_WATER * JSeaWater::O()));
    ionization.push_back(new JACoeffSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl()));
    ionization.push_back(new JACoeffSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H()));
    ionization.push_back(new JACoeffSource(Sodium,   DENSITY_SEA_WATER * JSeaWater::Na()));
#ifdef RADIATION
    ionization.push_back(new JACoeffSource(Calcium,  DENSITY_SEA_WATER * JSeaWater::Ca()));
    ionization.push_back(new JACoeffSource(Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg()));
    ionization.push_back(new JACoeffSource(Potassium,DENSITY_SEA_WATER * JSeaWater::K()));
    ionization.push_back(new JACoeffSource(Sulphur,  DENSITY_SEA_WATER * JSeaWater::S()));
#endif
    
    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 = APPLICATION_JSIRENE;
    buffer.simul.rbegin()->version = getGITVersion();
    buffer.simul.rbegin()->date    = getDate();
    buffer.simul.rbegin()->time    = getTime();

    buffer.detector.push_back(JAANET::detector());

    buffer.detector.rbegin()->program  = APPLICATION_JSIRENE;
    buffer.detector.rbegin()->filename = detectorFile;

    buffer.push(&JHead::simul);
    buffer.push(&JHead::detector);
    
    if (!keep) {

      if        (buffer.is_valid(&JHead::fixedcan)) {
        
        buffer.fixedcan.xcenter += center.x;
        buffer.fixedcan.ycenter += center.y;
        buffer.fixedcan.zmin    += center.z;
        buffer.fixedcan.zmax    += center.z;

        buffer.push(&JHead::fixedcan);
        
      } else if (buffer.is_valid(&JHead::can)) {

        buffer.can.zmin         += center.z;
        buffer.can.zmax         += center.z;

        buffer.push(&JHead::can);
      }

      buffer.coord_origin = coord_origin(0.0, 0.0, 0.0);
      
      buffer.push(&JHead::coord_origin);
    }
  
    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("Light generation volume: " << cylinder << endl);

  TH1D       job("job", NULL, 400,  0.5, 400.5);
  TProfile   cpu("cpu", NULL,  16,  0.0,   8.0);
  TProfile2D rms("rms", NULL,  16,  0.0,   8.0, 201, -0.5, 200.5);
  TProfile2D rad("rad", NULL,  16,  0.0,   8.0, 201, -0.5, 200.5);


  outputFile.open();

  if (!outputFile.is_open()) {
    FATAL("Error opening file " << outputFile << endl);
  }

  outputFile.put(meta);
  outputFile.put(header);
  outputFile.put(*gRandom);  

  const double         epsilon = 1.0e-6;                 // precision angle [rad]
  const JRange<double> pi(epsilon, PI - epsilon);        // constrain angle

  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 (!track->is_finalstate()) {
        continue;                                      // only final state particles produce light
      }

      if (is_muon(*track)) {

        // -----------------------------------------------
        //                      muon
        // -----------------------------------------------

        job.Fill(2.0);

        typedef JDetectorSubset<JModuleComparator3Z_t>   JDetectorSubset_t;

        const JCylinder3D::intersection_type intersection = cylinder.getIntersection(getAxis(*track));

        double Zmin = intersection.first;
        double Zmax = intersection.second;

        if (Zmax - Zmin <= parameters.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() <= parameters.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() >= parameters.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 [m^-1]
        
          for (int i = 0; i != N; ++i) {
            ls += li[i] = radiation[i]->getInverseInteractionLength(vertex.getE());
          }

          double As = 0.0;                                              // ionization energy loss

          for (size_t i = 0; i != ionization.size(); ++i) {
            As += ionization[i]->getA(vertex.getE());
          }

          double step  = gRandom->Exp(1.0) / ls;                        // distance to next radiation process
          double range = vertex.getRange(As);                           // range of muon

          if (vertex.getE() < parameters.Emax_GeV) {                    // limited step size
            if (parameters.Dmax_m < range) {
              range = parameters.Dmax_m;
            }
          }

          double ts = getThetaMCS(vertex.getE(), min(step,range));      // multiple Coulomb scattering angle [rad]
          double T2 = ts*ts;                                            // 

          rms.Fill(log10(vertex.getE()), (Double_t) 0, ts*ts);
          
          vertex.getDirection() += getRandomDirection(T2/3.0);          // multiple Coulomb planar scattering

          vertex.step(As, min(step,range));                             // ionization energy loss

          double Es = 0.0;                                              // shower energy [GeV]

          if (step < range) {

            if (vertex.getE() >= parameters.Emin_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]
                  ts = radiation[i]->getThetaRMS(vertex.getE(), Es);    // scattering angle [rad]

                  T2 += ts*ts;

                  rms.Fill(log10(vertex.getE()), (Double_t) radiation[i]->getID(), ts*ts);
                  rad.Fill(log10(vertex.getE()), (Double_t) radiation[i]->getID(), Es);
                  
                  break;
                }
              }
            }
          }

          vertex.applyEloss(getRandomDirection(T2), Es);

          muon.push_back(vertex);
        }

        // add muon end point

        if (vertex.getZ() < Zmax && 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 Z  = module->getZ() - module->getX() / getTanThetaC();

          if (Z >= muon.begin()->getZ() && Z <= muon.rbegin()->getZ()) {
          
            const JVector2D pos = muon.getPosition(Z);
            const double    R   = hypot(module->getX() - pos.getX(), 
                                        module->getY() - pos.getY());

            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) {
                    
                    vector<double> npe;

                    for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {

                      const double R     = hypot(pmt->getX() - pos.getX(), 
                                                 pmt->getY() - pos.getY());
                      const double theta = pi.constrain(pmt->getTheta());
                      const double phi   = pi.constrain(fabs(pmt->getPhi()));
                    
                      npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
                    }

                    const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);

                    for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {

                      const double R     = hypot(pmt->getX() - pos.getX(), 
                                                 pmt->getY() - pos.getY());
                      const double Z     = pmt->getZ() - z0;
                      const double theta = pi.constrain(pmt->getTheta());
                      const double phi   = pi.constrain(fabs(pmt->getPhi()));
                    
                      size_t n0 = min(ns[distance(module->begin(),pmt)], parameters.Nmax_PMT);
                    
                      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 (std::accumulate(npe.begin(), npe.end(), 0.0) > 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 Es = vertex->getEs();
          
          if (Es >= parameters.Ecut_GeV) {

            const double z0 = vertex->getZ();
            const double t0 = vertex->getT();
            const double DZ = geanz.getMaximum(Es);

            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 R  = hypot(module->getX() - vertex->getX(),
                                      module->getY() - vertex->getY());
              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) {
                      
                      vector<double> npe;

                      for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {

                        const double R     = hypot(pmt->getX() - vertex->getX(), 
                                                   pmt->getY() - vertex->getY());
                        const double Z     = pmt->getZ() - z0 - DZ;
                        const double D     = sqrt(R*R + Z*Z);
                        const double cd    = Z / D; 
                        const double theta = pi.constrain(pmt->getTheta());
                        const double phi   = pi.constrain(fabs(pmt->getPhi()));
                        
                        npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor);
                      }

                      const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
                        
                      for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {

                        const double R     = hypot(pmt->getX() - vertex->getX(), 
                                                   pmt->getY() - vertex->getY());
                        const double theta = pi.constrain(pmt->getTheta());
                        const double phi   = pi.constrain(fabs(pmt->getPhi()));
                        
                        size_t n0 = min(ns[distance(module->begin(),pmt)], parameters.Nmax_PMT);

                        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 (std::accumulate(npe.begin(), npe.end(), 0.0) > 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 (track->len > 0.0) {

        // -----------------------------------------------
        //      decayed particles treated as mip (tau includes mip+deltaray) 
        // -----------------------------------------------

        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) {
                
            double W = 1.0;                                   // mip

            if (is_deltarays(CDF[i].type)) {

              if (is_tau(*track)) 
                W = getDeltaRaysFromTau(E);                   // delta-rays
              else
                continue;
            }    

            if (R < CDF[i].integral.getXmax()) {
              
              try {
                
                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);

                  vector<double> npe;
                  
                  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {

                    const double R     = pmt->getX();
                    const double theta = pi.constrain(pmt->getTheta());
                    const double phi   = pi.constrain(fabs(pmt->getPhi()));
                    
                    npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
                  }

                  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);

                  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 = pi.constrain(pmt->getTheta());
                    const double phi   = pi.constrain(fabs(pmt->getPhi()));
                    
                    size_t n0 = min(ns[distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
                    
                    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 (std::accumulate(npe.begin(), npe.end(), 0.0) > 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 >= parameters.Ecut_GeV && cylinder.getDistance(getPosition(*track)) < parameters.Dmin_m) {
        
            const double z0 = 0.0;
            const double t0 = track->t;
            const double DZ = geanz.getMaximum(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() - 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) * E * module->size() * factor;
                    const int    N   = gRandom->Poisson(NPE); 

                    if (N != 0) {

                      buffer = *module;
            
                      buffer.transform(transformation);
          
                      vector<double> npe;
                        
                      for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
                    
                        const double R     = pmt->getX();
                        const double Z     = pmt->getZ() - z0 - DZ;
                        const double D     = sqrt(R*R + Z*Z);
                        const double cd    = Z / D; 
                        const double theta = pi.constrain(pmt->getTheta());
                        const double phi   = pi.constrain(fabs(pmt->getPhi()));
                          
                        npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * E * factor);
                      }
                          
                      const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);

                      for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
                    
                        const double theta = pi.constrain(pmt->getTheta());
                        const double phi   = pi.constrain(fabs(pmt->getPhi()));
                          
                        size_t n0 = min(ns[distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
                 
                        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 (std::accumulate(npe.begin(), npe.end(), 0.0) > 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(parameters.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(job);
  outputFile.put(cpu);
  outputFile.put(rms);
  outputFile.put(rad);
  outputFile.put(*gRandom);

  JMultipleFileScanner<JTYPELIST<JMetaTypes_t, MultiHead>::typelist> io(inputFile);

  io >> outputFile;

  outputFile.close();
}