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

#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/JDate.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.
 *
 * 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(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);

    //More accuracy can be achieved uncommenting the commented elements and its initializations, remember to do the same in JSeawater.hh (The code will run slower).

    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);
    /*  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);
    */

    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));
    /*  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));
    */   

    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));
    /*  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)); 
    */

    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)); 
    /*  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));
    */

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

      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("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;                 //
  const JRange<double> python(epsilon, PI - epsilon);    // constrain PMT angles

  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) {
		    
		    double ns = 0.0;
		  
		    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 = python.constrain(pmt->getTheta());
		      const double phi   = python.constrain(fabs(pmt->getPhi()));
		    
		      const double npe   = CDF[i].function.getNPE(R, theta, phi) * factor * W;
		    
		      ns += npe;
		    
		      int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), parameters.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 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) {
		      
		      double ns = 0.0;

		      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 = python.constrain(pmt->getTheta());
			const double phi   = python.constrain(fabs(pmt->getPhi()));
			
			const double npe   = CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor;
		      
			ns += npe;
			
			int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), parameters.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 (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);

		  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 = python.constrain(pmt->getTheta());
		    const double phi   = python.constrain(fabs(pmt->getPhi()));
		    
		    const double npe   = CDF[i].function.getNPE(R, theta, phi) * factor * W;
		    
		    ns += npe;
		    
		    int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), parameters.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 >= 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);
	  
		      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 - DZ;
			const double D     = sqrt(R*R + Z*Z);
			const double cd    = Z / D; 
			const double theta = python.constrain(pmt->getTheta());
			const double phi   = python.constrain(fabs(pmt->getPhi()));
			  
			const double npe   = CDG[i].function.getNPE(D, cd, theta, phi) * E * factor;

			ns += npe;
			  
			int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), parameters.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(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();
}