Jpp  master_rocky-40-g5f0272dcd
the software that should make you happy
JSirene.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <limits>
6 #include <numeric>
7 
8 #include "TROOT.h"
9 #include "TFile.h"
10 #include "TH1D.h"
11 #include "TH2D.h"
12 #include "TProfile.h"
13 #include "TProfile2D.h"
14 #include "TRandom3.h"
15 
16 #include "km3net-dataformat/offline/Head.hh"
17 #include "km3net-dataformat/offline/MultiHead.hh"
18 #include "km3net-dataformat/offline/Evt.hh"
19 #include "km3net-dataformat/definitions/applications.hh"
20 #include "km3net-dataformat/definitions/trkmembers.hh"
21 
22 #include "JLang/JSharedPointer.hh"
23 #include "JLang/Jpp.hh"
24 #include "JLang/JPredicate.hh"
25 #include "JSystem/JDateAndTime.hh"
26 
27 #include "JPhysics/JCDFTable.hh"
28 #include "JPhysics/JPDFTypes.hh"
29 #include "JPhysics/JPDFToolkit.hh"
30 #include "JPhysics/JGeane.hh"
31 #include "JPhysics/JGeanz.hh"
32 #include "JPhysics/JRadiation.hh"
33 #include "JPhysics/JRadiationFunction.hh"
34 #include "JPhysics/JRadiationSource.hh"
35 #include "JPhysics/JACoeffSource.hh"
36 #include "JPhysics/JPhysicsToolkit.hh"
37 #include "JPhysics/JPhysicsSupportkit.hh"
38 #include "JTools/JFunction1D_t.hh"
39 #include "JTools/JFunctionalMap_t.hh"
40 
41 #include "JAAnet/JAAnetToolkit.hh"
42 #include "JAAnet/JPDB.hh"
43 
44 #include "JSirene/JSirene.hh"
45 #include "JSirene/pythia.hh"
46 #include "JSirene/JSeaWater.hh"
47 #include "JSirene/JCDFTable1D.hh"
48 #include "JSirene/JCDFTable2D.hh"
49 #include "JSirene/JSireneToolkit.hh"
50 
51 #include "JDetector/JDetector.hh"
52 #include "JDetector/JDetectorSubset.hh"
53 #include "JDetector/JDetectorToolkit.hh"
54 #include "JGeometry3D/JCylinder3D.hh"
55 
56 #include "JSupport/JMultipleFileScanner.hh"
57 #include "JSupport/JFileRecorder.hh"
58 #include "JSupport/JMonteCarloFileSupportkit.hh"
59 #include "JSupport/JSupport.hh"
60 #include "JSupport/JMeta.hh"
61 
62 #include "Jeep/JProperties.hh"
63 #include "Jeep/JPrint.hh"
64 #include "Jeep/JTimer.hh"
65 #include "Jeep/JParser.hh"
66 #include "Jeep/JMessage.hh"
67 
68 
69 int debug; //!< debug level
70 int numberOfBins = 200; //!< number of bins for average CDF integral of optical module
71 double safetyFactor = 1.7; //!< safety factor for average CDF integral of optical module
72 
73 
74 namespace {
75 
76  using namespace JPP;
77 
78  typedef JHermiteSplineFunction1D_t JFunction1D_t;
79  typedef JMAPLIST<JPolint1FunctionalMap,
80  JPolint1FunctionalGridMap,
81  JPolint1FunctionalGridMap>::maplist J3DMap_t;
82  typedef JMAPLIST<JPolint1FunctionalMap,
83  JPolint1FunctionalMap,
84  JPolint1FunctionalGridMap,
85  JPolint1FunctionalGridMap>::maplist J4DMap_t;
86 
87  typedef JCDFTable<JFunction1D_t, J3DMap_t> JCDF4D_t; // muon
88  typedef JCDFTable<JFunction1D_t, J4DMap_t> JCDF5D_t; // shower
89 
90  typedef JCDFTable1D<JFunction1D_t::argument_type, JFunction1D_t::result_type> JCDF1D_t; // muon
91  typedef JCDFTable2D<JFunction1D_t::argument_type, JFunction1D_t::result_type> JCDF2D_t; // shower
92 
93 
94 
95  /**
96  * Auxiliary class for CDF tables.
97  */
98  template<class function_type, // CDF function
99  class integral_type> // CDF integral
100  struct JCDFHelper {
101  /**
102  * Constructor.
103  *
104  * \param file_descriptor file name descriptor
105  * \param type PDF type
106  */
107  JCDFHelper(const std::string& file_descriptor, const JPDFType_t type)
108  {
109  using namespace std;
110 
111  this->type = type;
112 
113  const string file_name = getFilename(file_descriptor, this->type);
114 
115  STATUS("loading input from file " << file_name << "... " << flush);
116 
117  try {
118  function.load(file_name.c_str());
119  }
120  catch(const JException& error) {
121  FATAL(error.what() << endl);
122  }
123 
124  new (&integral) integral_type(function, numberOfBins, safetyFactor);
125 
126  STATUS("OK" << endl);
127  }
128 
129  JPDFType_t type;
130  function_type function;
131  integral_type integral;
132  };
133 
134 
135  typedef JCDFHelper<JCDF4D_t, JCDF1D_t> JCDF_t; //!< muon light
136  typedef JCDFHelper<JCDF5D_t, JCDF2D_t> JCDG_t; //!< shower light
137 
138 
139  /**
140  * Get random direction.
141  *
142  * \param t2 square of theta RMS [rad^2]
143  * \return direction
144  */
145  inline JVersor3Z getRandomDirection(const double t2)
146  {
147  const double tv = sqrt(gRandom->Exp(1.0) * t2);
148  const double phi = gRandom->Uniform(-PI, +PI);
149 
150  return JVersor3Z(tv*cos(phi), tv*sin(phi));
151  }
152 
153 
154  /**
155  * Get Poisson number.
156  *
157  * \param x expectation value
158  * \return number
159  */
160  inline size_t getPoisson(const double x)
161  {
162  using namespace std;
163 
164  if (x < numeric_limits<int32_t>::max())
165  return (size_t) gRandom->Poisson(x);
166  else
167  return (size_t) gRandom->Gaus(x, sqrt(x));
168  }
169 }
170 
171 
172 /**
173  * \file
174  *
175  * Main program to simulate detector response to muons and showers.
176  *
177  * \image html sirene.png "Picture by Claudine Colnard"
178  *
179  * Note that CDF file descriptor should contain the wild card character JPHYSICS::WILDCARD;\n
180  * The file names are obtained by replacing JPHYSICS::WILDCARD; with
181  * - JPHYSICS::DIRECT_LIGHT_FROM_MUON;
182  * - JPHYSICS::SCATTERED_LIGHT_FROM_MUON;
183  * - JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS;
184  * - JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS;
185  * - JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER; and
186  * - JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER,
187  * respectively.
188  *
189  * More accuracy can be achieved by setting compile option RADITION but it will run slower.
190  *
191  * The CDF tables can be produced with the script <tt>JMakePDF.sh</tt>:
192  * <pre>
193  * JMakePDF.sh -P
194  * </pre>
195  * (option <tt>-h</tt> will print all available options).
196  * Note that the script will launch a large number of processes (<tt>JMakePDF</tt> and <tt>JMakePDG</tt>)\n
197  * which may take a considerable amount of time to completion.\n
198  * On a standard desktop, all jobs should be finished within 1/2 a day or so.
199  *
200  * The same script should then be run with option <tt>-M</tt> to merge the PDF files, i.e:
201  * <pre>
202  * JMakePDF.sh -M
203  * </pre>
204  *
205  * CDF tables are obtained by running the same script with option <tt>-C</tt>, i.e:
206  * <pre>
207  * JMakePDF.sh -C
208  * </pre>
209  *
210  * The various PDFs can be drawn using the <tt>JDrawPDF</tt> or <tt>JDrawPDG</tt> applications.\n
211  * The tabulated PDFs can be plotted using the <tt>JPlotPDF</tt> or <tt>JPlotPDG</tt> applications.\n
212  * The tabulated CDFs can be plotted using the <tt>JPlotCDF</tt> or <tt>JPlotCDG</tt> applications.
213  * \author mdejong
214  */
215 int main(int argc, char **argv)
216 {
217  using namespace std;
218  using namespace JPP;
219 
220  string fileDescriptor;
221  JMultipleFileScanner<JAAnetTypes_t> inputFile;
222  JFileRecorder <JTYPELIST<JAAnetTypes_t, JMetaTypes_t, JRootTypes_t>::typelist> outputFile;
223  JLimit_t& numberOfEvents = inputFile.getLimit();
224  string detectorFile;
225  JParameters parameters;
226  bool writeEMShowers;
227  size_t numberOfHits;
228  double factor;
229  UInt_t seed;
230 
231  try {
232 
233  JProperties properties;
234 
235  properties.insert(gmake_property(parameters.Ecut_GeV));
236  properties.insert(gmake_property(parameters.Emin_GeV));
237  properties.insert(gmake_property(parameters.Dmin_m));
238  properties.insert(gmake_property(parameters.Emax_GeV));
239  properties.insert(gmake_property(parameters.Dmax_m));
240  properties.insert(gmake_property(parameters.Tmax_ns));
241  properties.insert(gmake_property(parameters.Nmax_NPE));
242  properties.insert(gmake_property(parameters.Nmax_PMT));
243 
244  properties.insert(gmake_property(numberOfBins));
245  properties.insert(gmake_property(safetyFactor));
246 
247  JParser<> zap("Main program to simulate detector response to muons and showers.");
248 
249  zap['@'] = make_field(properties) = JPARSER::initialised();
250  zap['F'] = make_field(fileDescriptor, "file name descriptor for CDF tables");
251  zap['f'] = make_field(inputFile) = JPARSER::initialised();
252  zap['o'] = make_field(outputFile) = "sirene.root";
253  zap['n'] = make_field(numberOfEvents) = JLimit::max();
254  zap['a'] = make_field(detectorFile) = "";
255  zap['s'] = make_field(writeEMShowers, "store generated EM showers in event");
256  zap['N'] = make_field(numberOfHits, "minimum number of hits to output event") = 1;
257  zap['U'] = make_field(factor, "scaling factor applied to light yields") = 1.0;
258  zap['S'] = make_field(seed) = 0;
259  zap['d'] = make_field(debug) = 1;
260 
261  zap(argc, argv);
262  }
263  catch(const exception &error) {
264  FATAL(error.what() << endl);
265  }
266 
267 
268  gRandom->SetSeed(seed);
269 
270 
271  const JMeta meta(argc, argv);
272 
273  const double Zbed = 0.0; // level of seabed [m]
274 
275  vector<JCDF_t> CDF;
276  vector<JCDG_t> CDG;
277 
278  CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_MUON));
279  CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_MUON));
280  CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_DELTARAYS));
281  CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_DELTARAYS));
282 
283  CDG.push_back(JCDG_t(fileDescriptor, DIRECT_LIGHT_FROM_EMSHOWER));
284  CDG.push_back(JCDG_t(fileDescriptor, SCATTERED_LIGHT_FROM_EMSHOWER));
285 
286 
287  double maximal_road_width = 0.0; // road width [m]
288  double maximal_distance = 0.0; // road width [m]
289 
290  for (size_t i = 0; i != CDF.size(); ++i) {
291 
292  DEBUG("Range CDF["<< CDF[i].type << "] " << CDF[i].function.intensity.getXmax() << " m" << endl);
293 
294  maximal_road_width = max(maximal_road_width, CDF[i].function.intensity.getXmax());
295  }
296 
297  for (size_t i = 0; i != CDG.size(); ++i) {
298 
299  DEBUG("Range CDG["<< CDG[i].type << "] " << CDG[i].function.intensity.getXmax() << " m" << endl);
300 
301  if (!is_scattered(CDF[i].type)) {
302  maximal_road_width = max(maximal_road_width, CDG[i].function.intensity.getXmax());
303  }
304 
305  maximal_distance = max(maximal_distance, CDG[i].function.intensity.getXmax());
306  }
307 
308  NOTICE("Maximal road width [m] " << maximal_road_width << endl);
309  NOTICE("Maximal distance [m] " << maximal_distance << endl);
310 
311 
312  if (detectorFile == "" || inputFile.empty()) {
313  STATUS("Nothing to be done." << endl);
314  return 0;
315  }
316 
317  JDetector detector;
318 
319  try {
320 
321  STATUS("Load detector... " << flush);
322 
323  load(detectorFile, detector);
324 
325  STATUS("OK" << endl);
326  }
327  catch(const JException& error) {
328  FATAL(error);
329  }
330 
331  // remove empty modules
332 
333  for (JDetector::iterator module = detector.begin(); module != detector.end(); ) {
334  if (!module->empty())
335  ++module;
336  else
337  module = detector.erase(module);
338  }
339 
340  vector< JSharedPointer<JRadiationInterface> > radiation;
341  vector< JSharedPointer<JACoeffInterface> > ionization;
342 
343  if (true) {
344 
345  STATUS("Setting up radiation tables... " << flush);
346 
347  const JRadiation hydrogen ( 1.0, 1.0, 40, 0.01, 0.1, 0.1);
348  const JRadiation oxygen ( 8.0, 16.0, 40, 0.01, 0.1, 0.1);
349  const JRadiation chlorine (17.0, 35.5, 40, 0.01, 0.1, 0.1);
350  const JRadiation sodium (11.0, 23.0, 40, 0.01, 0.1, 0.1);
351 #ifdef RADIATION
352  const JRadiation calcium (20.0, 40.0, 40, 0.01, 0.1, 0.1);
353  const JRadiation magnesium(12.0, 24.3, 40, 0.01, 0.1, 0.1);
354  const JRadiation potassium(19.0, 39.0, 40, 0.01, 0.1, 0.1);
355  const JRadiation sulphur (16.0, 32.0, 40, 0.01, 0.1, 0.1);
356 #endif
357 
358  JSharedPointer<JRadiation> Hydrogen (new JRadiationFunction(hydrogen, 300, 0.2, 1.0e11));
359  JSharedPointer<JRadiation> Oxygen (new JRadiationFunction(oxygen, 300, 0.2, 1.0e11));
360  JSharedPointer<JRadiation> Chlorine (new JRadiationFunction(chlorine, 300, 0.2, 1.0e11));
361  JSharedPointer<JRadiation> Sodium (new JRadiationFunction(sodium, 300, 0.2, 1.0e11));
362 #ifdef RADIATION
363  JSharedPointer<JRadiation> Calcium (new JRadiationFunction(calcium, 300, 0.2, 1.0e11));
364  JSharedPointer<JRadiation> Magnesium(new JRadiationFunction(magnesium,300, 0.2, 1.0e11));
365  JSharedPointer<JRadiation> Potassium(new JRadiationFunction(potassium,300, 0.2, 1.0e11));
366  JSharedPointer<JRadiation> Sulphur (new JRadiationFunction(sulphur, 300, 0.2, 1.0e11));
367 #endif
368 
369  radiation.push_back(new JRadiationSource(11, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), EErad_t));
370  radiation.push_back(new JRadiationSource(12, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad_t));
371  radiation.push_back(new JRadiationSource(13, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), EErad_t));
372  radiation.push_back(new JRadiationSource(14, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), EErad_t));
373 #ifdef RADIATION
374  radiation.push_back(new JRadiationSource(15, Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), EErad_t));
375  radiation.push_back(new JRadiationSource(16, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), EErad_t));
376  radiation.push_back(new JRadiationSource(17, Potassium,DENSITY_SEA_WATER * JSeaWater::K(), EErad_t));
377  radiation.push_back(new JRadiationSource(18, Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), EErad_t));
378 #endif
379 
380  radiation.push_back(new JRadiationSource(21, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), Brems_t));
381  radiation.push_back(new JRadiationSource(22, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems_t));
382  radiation.push_back(new JRadiationSource(23, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), Brems_t));
383  radiation.push_back(new JRadiationSource(24, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), Brems_t));
384 #ifdef RADIATION
385  radiation.push_back(new JRadiationSource(25, Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), Brems_t));
386  radiation.push_back(new JRadiationSource(26, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), Brems_t));
387  radiation.push_back(new JRadiationSource(27, Potassium,DENSITY_SEA_WATER * JSeaWater::K(), Brems_t));
388  radiation.push_back(new JRadiationSource(28, Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), Brems_t));
389 #endif
390 
391  radiation.push_back(new JRadiationSource(31, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), GNrad_t));
392  radiation.push_back(new JRadiationSource(32, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad_t));
393  radiation.push_back(new JRadiationSource(33, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), GNrad_t));
394  radiation.push_back(new JRadiationSource(34, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), GNrad_t));
395 #ifdef RADIATION
396  radiation.push_back(new JRadiationSource(35, Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), GNrad_t));
397  radiation.push_back(new JRadiationSource(36, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), GNrad_t));
398  radiation.push_back(new JRadiationSource(37, Potassium,DENSITY_SEA_WATER * JSeaWater::K(), GNrad_t));
399  radiation.push_back(new JRadiationSource(38, Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), GNrad_t));
400 #endif
401 
402  radiation.push_back(new JDISSource(100, DENSITY_SEA_WATER));
403 
404  ionization.push_back(new JACoeffSource(Oxygen, DENSITY_SEA_WATER * JSeaWater::O()));
405  ionization.push_back(new JACoeffSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl()));
406  ionization.push_back(new JACoeffSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H()));
407  ionization.push_back(new JACoeffSource(Sodium, DENSITY_SEA_WATER * JSeaWater::Na()));
408 #ifdef RADIATION
409  ionization.push_back(new JACoeffSource(Calcium, DENSITY_SEA_WATER * JSeaWater::Ca()));
410  ionization.push_back(new JACoeffSource(Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg()));
411  ionization.push_back(new JACoeffSource(Potassium,DENSITY_SEA_WATER * JSeaWater::K()));
412  ionization.push_back(new JACoeffSource(Sulphur, DENSITY_SEA_WATER * JSeaWater::S()));
413 #endif
414 
415  STATUS("OK" << endl);
416  }
417 
418 
419  JCylinder3D cylinder(detector.begin(), detector.end());
420 
421  cylinder.addMargin(maximal_distance);
422 
423  if (cylinder.getZmin() < Zbed) {
424  cylinder.setZmin(Zbed);
425  }
426 
427  NOTICE("Light generation volume: " << cylinder << endl);
428 
429 
430  Vec offset(0,0,0);
431  Head header;
432 
433  try {
434 
435  header = inputFile.getHeader();
436 
437  JHead buffer(header);
438 
439  buffer.simul.push_back(JAANET::simul());
440 
441  buffer.simul.rbegin()->program = APPLICATION_JSIRENE;
442  buffer.simul.rbegin()->version = getGITVersion();
443  buffer.simul.rbegin()->date = getDate();
444  buffer.simul.rbegin()->time = getTime();
445 
446  buffer.push(&JHead::simul);
447 
448  buffer.detector.push_back(JAANET::detector());
449 
450  buffer.detector.rbegin()->program = APPLICATION_JSIRENE;
451  buffer.detector.rbegin()->filename = detectorFile;
452 
453  buffer.push(&JHead::detector);
454 
455  offset += Vec(cylinder.getX(), cylinder.getY(), 0.0);
456  offset -= getOrigin(buffer);
457 
458  if (buffer.is_valid(&JHead::fixedcan)) {
459 
460  buffer.fixedcan.xcenter += offset.x;
461  buffer.fixedcan.ycenter += offset.y;
462  buffer.fixedcan.zmin += offset.z;
463  buffer.fixedcan.zmax += offset.z;
464 
465  } else {
466 
467  buffer.fixedcan.xcenter = cylinder.getX();
468  buffer.fixedcan.ycenter = cylinder.getY();
469 
470  if (buffer.is_valid(&JHead::can)) {
471 
472  buffer.fixedcan.radius = buffer.can.r;
473  buffer.fixedcan.zmin = buffer.can.zmin + offset.z;
474  buffer.fixedcan.zmax = buffer.can.zmax + offset.z;
475  } else {
476 
477  buffer.fixedcan.radius = cylinder.getRadius();
478  buffer.fixedcan.zmin = cylinder.getZmin();
479  buffer.fixedcan.zmax = cylinder.getZmax();
480  }
481  }
482 
483  buffer.push(&JHead::fixedcan);
484 
485  if (buffer.is_valid(&JHead::coord_origin)) {
486 
487  buffer.coord_origin = coord_origin(0.0, 0.0, 0.0);
488 
489  buffer.push(&JHead::coord_origin);
490  }
491 
492  copy(buffer, header);
493  }
494  catch(const JException& error) {
495  FATAL(error);
496  }
497 
498  NOTICE("Offset applied to true tracks is: " << offset << endl);
499 
500  TH1D job("job", NULL, 400, 0.5, 400.5);
501  TProfile cpu("cpu", NULL, 16, 0.0, 8.0);
502  TProfile2D rms("rms", NULL, 16, 0.0, 8.0, 201, -0.5, 200.5);
503  TProfile2D rad("rad", NULL, 16, 0.0, 8.0, 201, -0.5, 200.5);
504 
505 
506  outputFile.open();
507 
508  if (!outputFile.is_open()) {
509  FATAL("Error opening file " << outputFile << endl);
510  }
511 
512  outputFile.put(meta);
513  outputFile.put(header);
514  outputFile.put(*gRandom);
515 
516  const double epsilon = 1.0e-6; // precision angle [rad]
517  const JRange<double> pi(epsilon, PI - epsilon); // constrain angle
518 
519  JTimer timer;
520 
521  for (JMultipleFileScanner<Evt>& in = inputFile; in.hasNext(); ) {
522 
523  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
524 
525  job.Fill(1.0);
526 
527  Evt* evt = in.next();
528 
529  for (vector<Trk>::iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
530  track->pos += offset;
531  }
532 
533  Evt event(*evt); // output
534 
535  event.mc_hits.clear();
536 
537  JHits_t mc_hits; // temporary buffer
538 
539  timer.reset();
540  timer.start();
541 
542  for (vector<Trk>::const_iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
543 
544  if (!track->is_finalstate()) {
545  continue; // only final state particles produce light
546  }
547 
548  if (is_muon(*track)) {
549 
550  // -----------------------------------------------
551  // muon
552  // -----------------------------------------------
553 
554  job.Fill(2.0);
555 
556  typedef JDetectorSubset<JModuleComparator3Z_t> JDetectorSubset_t;
557 
558  const JCylinder3D::intersection_type intersection = cylinder.getIntersection(getAxis(*track));
559 
560  double Zmin = intersection.first;
561  double Zmax = intersection.second;
562 
563  if (Zmax - Zmin <= parameters.Dmin_m) {
564  continue;
565  }
566 
567  JVertex vertex(0.0, track->t, track->E); // start of muon
568 
569  if (track->pos.z < Zbed) { // propagate muon through rock
570 
571  if (track->dir.z > 0.0)
572  vertex.step(gRock, (Zbed - track->pos.z) / track->dir.z);
573  else
574  continue;
575  }
576 
577  if (vertex.getZ() < Zmin) { // propagate muon through water
578  vertex.step(gWater, Zmin - vertex.getZ());
579  }
580 
581  if (vertex.getRange() <= parameters.Dmin_m) {
582  continue;
583  }
584 
585  job.Fill(3.0);
586 
587  const JDetectorSubset_t subdetector(detector, getAxis(*track), maximal_road_width);
588 
589  if (subdetector.empty()) {
590  continue;
591  }
592 
593  job.Fill(4.0);
594 
595  JTrack muon(vertex); // propagate muon trough detector
596 
597  while (vertex.getE() >= parameters.Emin_GeV && vertex.getZ() < Zmax) {
598 
599  const int N = radiation.size();
600 
601  double li[N]; // inverse interaction lengths
602  double ls = 1.0e-5; // minimal total inverse interaction length [m^-1]
603 
604  for (int i = 0; i != N; ++i) {
605  ls += li[i] = radiation[i]->getInverseInteractionLength(vertex.getE());
606  }
607 
608  double As = 0.0; // ionization energy loss
609 
610  for (size_t i = 0; i != ionization.size(); ++i) {
611  As += ionization[i]->getA(vertex.getE());
612  }
613 
614  double step = gRandom->Exp(1.0) / ls; // distance to next radiation process
615  double range = vertex.getRange(As); // range of muon
616 
617  if (vertex.getE() < parameters.Emax_GeV) { // limited step size
618  if (parameters.Dmax_m < range) {
619  range = parameters.Dmax_m;
620  }
621  }
622 
623  double ts = getThetaMCS(vertex.getE(), min(step,range)); // multiple Coulomb scattering angle [rad]
624  double T2 = ts*ts; //
625 
626  rms.Fill(log10(vertex.getE()), (Double_t) 0, ts*ts);
627 
628  vertex.getDirection() += getRandomDirection(T2/3.0); // multiple Coulomb planar scattering
629 
630  vertex.step(As, min(step,range)); // ionization energy loss
631 
632  double Es = 0.0; // shower energy [GeV]
633 
634  if (step < range) {
635 
636  if (vertex.getE() >= parameters.Emin_GeV) {
637 
638  double y = gRandom->Uniform(ls);
639 
640  for (int i = 0; i != N; ++i) {
641 
642  y -= li[i];
643 
644  if (y < 0.0) {
645 
646  Es = radiation[i]->getEnergyOfShower(vertex.getE()); // shower energy [GeV]
647  ts = radiation[i]->getThetaRMS(vertex.getE(), Es); // scattering angle [rad]
648 
649  T2 += ts*ts;
650 
651  rms.Fill(log10(vertex.getE()), (Double_t) radiation[i]->getID(), ts*ts);
652  rad.Fill(log10(vertex.getE()), (Double_t) radiation[i]->getID(), Es);
653 
654  break;
655  }
656  }
657  }
658  }
659 
660  vertex.applyEloss(getRandomDirection(T2), Es);
661 
662  muon.push_back(vertex);
663  }
664 
665  // add muon end point
666 
667  if (vertex.getZ() < Zmax && vertex.getRange() > 0.0) {
668 
669  vertex.step(vertex.getRange());
670 
671  muon.push_back(vertex);
672  }
673 
674  // add information to output muon
675 
676  vector<Trk>::iterator trk = find_if(event.mc_trks.begin(),
677  event.mc_trks.end(),
678  make_predicate(&Trk::id, track->id));
679 
680  if (trk != event.mc_trks.end()) {
681  trk->len = (muon.rbegin()->getZ() < Zmax ? +1 : -1) * (muon.rbegin()->getZ() - muon.begin()->getZ());
682  trk->setusr(mc_usr_keys::energy_lost_in_can, muon.begin()->getE() - muon.rbegin()->getE());
683  }
684 
685  for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
686 
687  const double z0 = muon.begin()->getZ();
688  const double t0 = muon.begin()->getT();
689  const double Z = module->getZ() - module->getX() / getTanThetaC();
690 
691  if (Z >= muon.begin()->getZ() && Z <= muon.rbegin()->getZ()) {
692 
693  const JVector2D pos = muon.getPosition(Z);
694  const double R = hypot(module->getX() - pos.getX(),
695  module->getY() - pos.getY());
696 
697  for (size_t i = 0; i != CDF.size(); ++i) {
698 
699  if (R < CDF[i].integral.getXmax()) {
700 
701  try {
702 
703  double W = 1.0; // mip
704 
705  if (is_deltarays(CDF[i].type)) {
706  W = getDeltaRaysFromMuon(muon.getE(Z)); // delta-rays
707  }
708 
709  const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
710  const size_t N = getPoisson(NPE);
711 
712  if (N != 0) {
713 
714  vector<double> npe;
715 
716  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
717 
718  const double R = hypot(pmt->getX() - pos.getX(),
719  pmt->getY() - pos.getY());
720  const double theta = pi.constrain(pmt->getTheta());
721  const double phi = pi.constrain(fabs(pmt->getPhi()));
722 
723  npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
724  }
725 
726  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
727 
728  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
729 
730  const double R = hypot(pmt->getX() - pos.getX(),
731  pmt->getY() - pos.getY());
732  const double Z = pmt->getZ() - z0;
733  const double theta = pi.constrain(pmt->getTheta());
734  const double phi = pi.constrain(fabs(pmt->getPhi()));
735 
736  size_t n0 = min(ns[distance(module->begin(),pmt)], parameters.Nmax_PMT);
737 
738  job.Fill((double) (100 + CDF[i].type), (double) n0);
739 
740  while (n0 != 0) {
741 
742  const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
743  const int n1 = getNumberOfPhotoElectrons(n0);
744 
745  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
746  pmt->getID(),
747  getHitType(CDF[i].type),
748  track->id,
749  t0 + (R * getTanThetaC() + Z) / C + t1,
750  n1));
751 
752  n0 -= n1;
753  }
754  }
755 
756  if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
757  job.Fill((double) (300 + CDF[i].type));
758  }
759  }
760  }
761  catch(const exception& error) {
762  job.Fill((double) (200 + CDF[i].type));
763  }
764  }
765  }
766  }
767  }
768 
769  for (JTrack::const_iterator vertex = muon.begin(); vertex != muon.end(); ++vertex) {
770 
771  const double Es = vertex->getEs();
772 
773  if (Es >= parameters.Ecut_GeV) {
774 
775  const double z0 = vertex->getZ();
776  const double t0 = vertex->getT();
777  const double DZ = geanz.getMaximum(Es);
778 
779  int origin = track->id;
780 
781  if (writeEMShowers) {
782  origin = event.mc_trks.size() + 1;
783  }
784 
785  int number_of_hits = 0;
786 
787  JDetectorSubset_t::range_type range = subdetector.getRange(z0 - maximal_distance,
788  z0 + maximal_distance);
789 
790  for (JDetector::const_iterator module = range.begin(); module != range.end(); ++module) {
791 
792  const double R = hypot(module->getX() - vertex->getX(),
793  module->getY() - vertex->getY());
794  const double Z = module->getZ() - z0 - DZ;
795  const double D = sqrt(R*R + Z*Z);
796  const double cd = Z / D;
797 
798  for (size_t i = 0; i != CDG.size(); ++i) {
799 
800  if (D < CDG[i].integral.getXmax()) {
801 
802  try {
803 
804  const double NPE = CDG[i].integral.getNPE(D, cd) * Es * module->size() * factor;
805  const size_t N = getPoisson(NPE);
806 
807  if (N != 0) {
808 
809  vector<double> npe;
810 
811  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
812 
813  const double R = hypot(pmt->getX() - vertex->getX(),
814  pmt->getY() - vertex->getY());
815  const double Z = pmt->getZ() - z0 - DZ;
816  const double D = sqrt(R*R + Z*Z);
817  const double cd = Z / D;
818  const double theta = pi.constrain(pmt->getTheta());
819  const double phi = pi.constrain(fabs(pmt->getPhi()));
820 
821  npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor);
822  }
823 
824  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
825 
826  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
827 
828  const double R = hypot(pmt->getX() - vertex->getX(),
829  pmt->getY() - vertex->getY());
830  const double theta = pi.constrain(pmt->getTheta());
831  const double phi = pi.constrain(fabs(pmt->getPhi()));
832 
833  size_t n0 = min(ns[distance(module->begin(),pmt)], parameters.Nmax_PMT);
834 
835  job.Fill((double) (100 + CDG[i].type), (double) n0);
836 
837  while (n0 != 0) {
838 
839  const double dz = geanz.getLength(Es, gRandom->Rndm());
840  const double Z = pmt->getZ() - z0 - dz;
841  const double D = sqrt(R*R + Z*Z);
842  const double cd = Z / D;
843 
844  const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
845  const int n1 = getNumberOfPhotoElectrons(n0);
846 
847  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
848  pmt->getID(),
849  getHitType(CDG[i].type),
850  origin,
851  t0 + (dz + D * getIndexOfRefraction()) / C + t1,
852  n1));
853 
854  n0 -= n1;
855 
856  number_of_hits += n1;
857  }
858  }
859 
860  if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
861  job.Fill((double) (300 + CDG[i].type));
862  }
863  }
864  }
865  catch(const exception& error) {
866  job.Fill((double) (200 + CDG[i].type));
867  }
868  }
869  }
870  }
871 
872  if (writeEMShowers && number_of_hits != 0) {
873 
874  event.mc_trks.push_back(JTrk_t(origin,
875  TRACK_TYPE_PHOTON,
876  track->id,
877  track->pos + track->dir * vertex->getZ(),
878  track->dir,
879  vertex->getT(),
880  Es));
881  }
882  }
883  }
884 
885  } else if (track->len > 0.0) {
886 
887  // -----------------------------------------------
888  // decayed particles treated as mip (tau includes mip+deltaray)
889  // -----------------------------------------------
890 
891  job.Fill(6.0);
892 
893  const double z0 = 0.0;
894  const double z1 = z0 + track->len;
895  const double t0 = track->t;
896  const double E = track->E;
897 
898  const JTransformation3D transformation = getTransformation(*track);
899 
900  JModule buffer;
901 
902  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
903 
904  const JPosition3D pos = transformation.transform(module->getPosition());
905 
906  const double R = pos.getX();
907  const double Z = pos.getZ() - R / getTanThetaC();
908 
909  if (Z < z0 ||
910  Z > z1 ||
911  R > maximal_road_width) {
912  continue;
913  }
914 
915  for (size_t i = 0; i != CDF.size(); ++i) {
916 
917  double W = 1.0; // mip
918 
919  if (is_deltarays(CDF[i].type)) {
920 
921  if (is_tau(*track))
922  W = getDeltaRaysFromTau(E); // delta-rays
923  else
924  continue;
925  }
926 
927  if (R < CDF[i].integral.getXmax()) {
928 
929  try {
930 
931  const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
932  const size_t N = getPoisson(NPE);
933 
934  if (N != 0) {
935 
936  buffer = *module;
937 
938  buffer.transform(transformation);
939 
940  vector<double> npe;
941 
942  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
943 
944  const double R = pmt->getX();
945  const double theta = pi.constrain(pmt->getTheta());
946  const double phi = pi.constrain(fabs(pmt->getPhi()));
947 
948  npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
949  }
950 
951  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
952 
953  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
954 
955  const double R = pmt->getX();
956  const double Z = pmt->getZ() - z0;
957  const double theta = pi.constrain(pmt->getTheta());
958  const double phi = pi.constrain(fabs(pmt->getPhi()));
959 
960  size_t n0 = min(ns[distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
961 
962  job.Fill((double) (120 + CDF[i].type), (double) n0);
963 
964  while (n0 != 0) {
965 
966  const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
967  const int n1 = getNumberOfPhotoElectrons(n0);
968 
969  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
970  pmt->getID(),
971  getHitType(CDF[i].type),
972  track->id,
973  t0 + (R * getTanThetaC() + Z) / C + t1,
974  n1));
975 
976  n0 -= n1;
977  }
978  }
979 
980  if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
981  job.Fill((double) (320 + CDF[i].type));
982  }
983  }
984  }
985  catch(const exception& error) {
986  job.Fill((double) (220 + CDF[i].type));
987  }
988  }
989  }
990  }
991 
992  if (!buffer.empty()) {
993  job.Fill(7.0);
994  }
995 
996  } else if (!is_neutrino(*track)) {
997 
998  if (JPDB::getInstance().hasPDG(track->type)) {
999 
1000  // -----------------------------------------------
1001  // electron or hadron
1002  // -----------------------------------------------
1003 
1004  job.Fill(8.0);
1005 
1006  double E = track->E;
1007 
1008  try {
1009  E = getKineticEnergy(E, JPDB::getInstance().getPDG(track->type).mass);
1010  }
1011  catch(const exception& error) {
1012  ERROR(error.what() << endl);
1013  }
1014 
1015  E = pythia(track->type, E);
1016 
1017  if (E >= parameters.Ecut_GeV && cylinder.getDistance(getPosition(*track)) < parameters.Dmin_m) {
1018 
1019  const double z0 = 0.0;
1020  const double t0 = track->t;
1021  const double DZ = geanz.getMaximum(E);
1022 
1023  const JTransformation3D transformation = getTransformation(*track);
1024 
1025  JModule buffer;
1026 
1027  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
1028 
1029  const JPosition3D pos = transformation.transform(module->getPosition());
1030 
1031  const double R = pos.getX();
1032  const double Z = pos.getZ() - z0 - DZ;
1033  const double D = sqrt(R*R + Z*Z);
1034  const double cd = Z / D;
1035 
1036  for (size_t i = 0; i != CDG.size(); ++i) {
1037 
1038  if (D < CDG[i].integral.getXmax()) {
1039 
1040  try {
1041 
1042  const double NPE = CDG[i].integral.getNPE(D, cd) * E * module->size() * factor;
1043  const size_t N = getPoisson(NPE);
1044 
1045  if (N != 0) {
1046 
1047  buffer = *module;
1048 
1049  buffer.transform(transformation);
1050 
1051  vector<double> npe;
1052 
1053  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1054 
1055  const double R = pmt->getX();
1056  const double Z = pmt->getZ() - z0 - DZ;
1057  const double D = sqrt(R*R + Z*Z);
1058  const double cd = Z / D;
1059  const double theta = pi.constrain(pmt->getTheta());
1060  const double phi = pi.constrain(fabs(pmt->getPhi()));
1061 
1062  npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * E * factor);
1063  }
1064 
1065  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
1066 
1067  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1068 
1069  const double theta = pi.constrain(pmt->getTheta());
1070  const double phi = pi.constrain(fabs(pmt->getPhi()));
1071 
1072  size_t n0 = min(ns[distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
1073 
1074  job.Fill((double) (140 + CDG[i].type), (double) n0);
1075 
1076  while (n0 != 0) {
1077 
1078  const double dz = geanz.getLength(E, gRandom->Rndm());
1079  const double Z = pmt->getZ() - z0 - dz;
1080  const double D = sqrt(R*R + Z*Z);
1081  const double cd = Z / D;
1082 
1083  const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
1084  const int n1 = getNumberOfPhotoElectrons(n0);
1085 
1086  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
1087  pmt->getID(),
1088  getHitType(CDG[i].type, true),
1089  track->id,
1090  t0 + (dz + D * getIndexOfRefraction()) / C + t1,
1091  n1));
1092 
1093  n0 -= n1;
1094  }
1095  }
1096 
1097  if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
1098  job.Fill((double) (340 + CDG[i].type));
1099  }
1100  }
1101  }
1102  catch(const exception& error) {
1103  job.Fill((double) (240 + CDG[i].type));
1104  }
1105  }
1106  }
1107  }
1108 
1109  if (!buffer.empty()) {
1110  job.Fill(9.0);
1111  }
1112 
1113  } else {
1114  job.Fill(21.0);
1115  }
1116  }
1117  }
1118  }
1119 
1120  if (!mc_hits.empty()) {
1121 
1122  mc_hits.merge(parameters.Tmax_ns);
1123 
1124  event.mc_hits.resize(mc_hits.size());
1125 
1126  copy(mc_hits.begin(), mc_hits.end(), event.mc_hits.begin());
1127  }
1128 
1129  timer.stop();
1130 
1131  if (has_neutrino(event)) {
1132  cpu.Fill(log10(get_neutrino(event).E), (double) timer.usec_ucpu * 1.0e-3);
1133  }
1134 
1135  if (event.mc_hits.size() >= numberOfHits) {
1136 
1137  outputFile.put(event);
1138 
1139  job.Fill(10.0);
1140  }
1141  }
1142  STATUS(endl);
1143 
1144  outputFile.put(job);
1145  outputFile.put(cpu);
1146  outputFile.put(rms);
1147  outputFile.put(rad);
1148  outputFile.put(*gRandom);
1149 
1150  JMultipleFileScanner<JTYPELIST<JMetaTypes_t, MultiHead>::typelist> io(inputFile);
1151 
1152  io >> outputFile;
1153 
1154  outputFile.close();
1155 }
JAAnetToolkit.hh
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JDateAndTime.hh
Date and time functions.
JDetectorSubset.hh
Auxiliary class to extract a subset of optical modules from a detector.
JDetector.hh
Data structure for detector geometry and calibration.
JFileRecorder.hh
Recording of objects on file according a format that follows from the file name extension.
JFunction1D_t.hh
JFunctionalMap_t.hh
Various implementations of functional maps.
JGeane.hh
Energy loss of muon.
JGeanz.hh
Longitudinal emission profile EM-shower.
JMessage.hh
General purpose messaging.
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
JMeta.hh
ROOT I/O of application specific meta data.
JMultipleFileScanner.hh
Scanning of objects from multiple files according a format that follows from the extension of each fi...
JPDFToolkit.hh
Auxiliary methods for PDF calculations.
JPDFTypes.hh
Numbering scheme for PDF types.
JParser.hh
Utility class to parse command line options.
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
JPhysicsToolkit.hh
Auxiliary methods for physics calculations.
JPrint.hh
I/O formatting auxiliaries.
JProperties.hh
Utility class to parse parameter values.
gmake_property
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Definition: JProperties.hh:1293
JRadiation.hh
Muon radiative cross sections.
JSireneToolkit.hh
Toolkit for JSirene.
numberOfBins
int numberOfBins
number of bins for average CDF integral of optical module
Definition: JSirene.cc:70
safetyFactor
double safetyFactor
safety factor for average CDF integral of optical module
Definition: JSirene.cc:71
main
int main(int argc, char **argv)
Definition: JSirene.cc:215
debug
int debug
debug level
Definition: JSirene.cc:69
JSupport.hh
ROOT TTree parameter settings of various packages.
Jpp.hh
Jpp environment information.
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition: PhysicsEvent.hh:434
APPLICATION_JSIRENE
static const char *const APPLICATION_JSIRENE
detector simulation
Definition: applications.hh:19
JAANET::JHead
Monte Carlo run header.
Definition: JHead.hh:1236
JAANET::JHead::simul
std::vector< JAANET::simul > simul
Definition: JHead.hh:1591
JAANET::JHead::coord_origin
JAANET::coord_origin coord_origin
Definition: JHead.hh:1601
JAANET::JHead::push
void push(T JHead::*pd)
Push given data member to Head.
Definition: JHead.hh:1374
JAANET::JHead::detector
std::vector< JAANET::detector > detector
Definition: JHead.hh:1587
JAANET::JHead::can
JAANET::can can
Definition: JHead.hh:1598
JAANET::JHead::fixedcan
JAANET::fixedcan fixedcan
Definition: JHead.hh:1599
JAANET::JHead::is_valid
bool is_valid(T JHead::*pd) const
Check validity of given data member in JHead.
Definition: JHead.hh:1319
JDETECTOR::JDetectorSubset_t
Detector subset without binary search functionality.
Definition: JDetectorSubset.hh:40
JDETECTOR::JDetectorSubset
Detector subset with binary search functionality.
Definition: JDetectorSubset.hh:121
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:96
JDETECTOR::JModule
Data structure for a composite optical module.
Definition: JModule.hh:75
JDETECTOR::JModule::transform
void transform(const JRotation3D &R, const JVector3D &pos)
Transformation of geometry (see method JGEOMETRY3D::JPosition3D::transform(const JRotation3D&,...
Definition: JModule.hh:345
JEEP::JProperties
Utility class to parse parameter values.
Definition: JProperties.hh:501
JEEP::JTimer
Auxiliary class for CPU timing and usage.
Definition: JTimer.hh:33
JEEP::JTimer::usec_ucpu
unsigned long long usec_ucpu
Definition: JTimer.hh:239
JEEP::JTimer::stop
void stop()
Stop timer.
Definition: JTimer.hh:127
JEEP::JTimer::reset
void reset()
Reset timer.
Definition: JTimer.hh:93
JEEP::JTimer::start
void start()
Start timer.
Definition: JTimer.hh:106
JGEOMETRY2D::JCircle2D::getRadius
double getRadius() const
Get radius.
Definition: JCircle2D.hh:144
JGEOMETRY2D::JVector2D
Data structure for vector in two dimensions.
Definition: JVector2D.hh:34
JGEOMETRY2D::JVector2D::getY
double getY() const
Get y position.
Definition: JVector2D.hh:74
JGEOMETRY2D::JVector2D::getX
double getX() const
Get x position.
Definition: JVector2D.hh:63
JGEOMETRY3D::JCylinder3D
Cylinder object.
Definition: JCylinder3D.hh:41
JGEOMETRY3D::JCylinder3D::getZmin
double getZmin() const
Get minimal z position.
Definition: JCylinder3D.hh:111
JGEOMETRY3D::JCylinder3D::getIntersection
intersection_type getIntersection(const JAxis3D &axis) const
Get intersection points of axis with cylinder.
Definition: JCylinder3D.hh:277
JGEOMETRY3D::JCylinder3D::setZmin
void setZmin(const double zmin)
Set minimal z position.
Definition: JCylinder3D.hh:133
JGEOMETRY3D::JCylinder3D::addMargin
void addMargin(const double D)
Add (safety) margin.
Definition: JCylinder3D.hh:172
JGEOMETRY3D::JCylinder3D::getDistance
double getDistance(const JVector3D &pos) const
Get distance between cylinder wall and given position.
Definition: JCylinder3D.hh:222
JGEOMETRY3D::JCylinder3D::getZmax
double getZmax() const
Get maximal z position.
Definition: JCylinder3D.hh:122
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:38
JGEOMETRY3D::JTime::getT
double getT() const
Get time.
Definition: JGeometry3D/JTime.hh:133
JGEOMETRY3D::JVector3D::getY
double getY() const
Get y position.
Definition: JVector3D.hh:104
JGEOMETRY3D::JVector3D::getZ
double getZ() const
Get z position.
Definition: JVector3D.hh:115
JGEOMETRY3D::JVector3D::getX
double getX() const
Get x position.
Definition: JVector3D.hh:94
JGEOMETRY3D::JVersor3Z
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:41
JGEOMETRY3D::JVersor3Z::getDirection
const JVersor3Z & getDirection() const
Get direction.
Definition: JVersor3Z.hh:81
JLANG::JException
General exception.
Definition: JException.hh:24
JLANG::JException::what
virtual const char * what() const override
Get error message.
Definition: JException.hh:64
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1698
JPHYSICS::JACoeffSource
Implementation for calculation of ionization constant.
Definition: JACoeffSource.hh:51
JPHYSICS::JCDFTable1D
Custom class for CDF table in 1 dimension.
Definition: JCDFTable1D.hh:41
JPHYSICS::JCDFTable2D
Custom class for CDF table in 2 dimensions.
Definition: JCDFTable2D.hh:43
JPHYSICS::JCDFTable
Multi-dimensional CDF table for arrival time of Cherenkov light.
Definition: JCDFTable.hh:57
JPHYSICS::JDISSource
Implementation for calculation of inverse interaction length and shower energy due to deep-inelastic ...
Definition: JRadiationSource.hh:163
JPHYSICS::JGeanz::getLength
double getLength(const double E, const double P, const double eps=1.0e-3) const
Get shower length for a given integrated probability.
Definition: JGeanz.hh:122
JPHYSICS::JGeanz::getMaximum
double getMaximum(const double E) const
Get depth of shower maximum.
Definition: JGeanz.hh:162
JPHYSICS::JRadiationFunction
Fast implementation of class JRadiation.
Definition: JRadiationFunction.hh:34
JPHYSICS::JRadiationSource
Implementation for calculation of inverse interaction length and shower energy.
Definition: JRadiationSource.hh:93
JPHYSICS::JRadiation
Auxiliary class for the calculation of the muon radiative cross sections.
Definition: JRadiation.hh:36
JSIRENE::JSeaWater::Cl
static double Cl()
Definition: JSeaWater.hh:29
JSIRENE::JSeaWater::Na
static double Na()
Definition: JSeaWater.hh:28
JSIRENE::JSeaWater::O
static double O()
Estimated mass fractions of chemical elements in sea water.
Definition: JSeaWater.hh:26
JSIRENE::JSeaWater::H
static double H()
Definition: JSeaWater.hh:27
JSUPPORT::JFileRecorder
Object writing to file.
Definition: JFileRecorder.hh:44
JSUPPORT::JMultipleFileScanner
General purpose class for object reading from a list of file names.
Definition: JMultipleFileScanner.hh:306
JSUPPORT::JMultipleFileScanner::hasNext
virtual bool hasNext() override
Check availability of next element.
Definition: JMultipleFileScanner.hh:409
JTOOLS::JRange::constrain
T constrain(argument_type x) const
Constrain value to range.
Definition: JRange.hh:350
K
static const uint32_t K[64]
Definition: crypt.cc:77
GAUSS_LEGENDRE::epsilon
const double epsilon
Definition: JQuadrature.cc:21
JAANET::getAxis
JAxis3D getAxis(const Trk &track)
Get axis.
Definition: JAAnetToolkit.hh:250
JAANET::get_neutrino
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Definition: JAAnetToolkit.hh:513
JAANET::getOrigin
Vec getOrigin(const JHead &header)
Get origin.
Definition: JHeadToolkit.hh:309
JAANET::getKineticEnergy
double getKineticEnergy(const Trk &trk)
Get track kinetic energy.
Definition: JAAnetToolkit.hh:962
JAANET::getTransformation
JTransformation3D getTransformation(const Trk &track)
Get transformation.
Definition: JAAnetToolkit.hh:274
JAANET::has_neutrino
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Definition: JAAnetToolkit.hh:491
JAANET::is_neutrino
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
Definition: JAAnetToolkit.hh:320
JAANET::getTime
double getTime(const Hit &hit)
Get true time of hit.
Definition: JAAnetToolkit.hh:96
JAANET::copy
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:162
JAANET::getPosition
JPosition3D getPosition(const Vec &pos)
Get position.
Definition: JAAnetToolkit.hh:178
JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition: JAAnetToolkit.hh:338
JAANET::getNPE
double getNPE(const Hit &hit)
Get true charge of hit.
Definition: JAAnetToolkit.hh:108
JAANET::TRACK_TYPE_PHOTON
@ TRACK_TYPE_PHOTON
Definition: JParticleTypes.hh:102
JAANET::is_tau
bool is_tau(const Trk &track)
Test whether given track is a (anti-)tau.
Definition: JAAnetToolkit.hh:346
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:505
JEEP::getFilename
std::string getFilename(const std::string &file_name)
Get file name part, i.e. part after last JEEP::PATHNAME_SEPARATOR if any.
Definition: JeepToolkit.hh:128
JEEP::getGITVersion
std::string getGITVersion(const std::string &tag)
Get GIT version for given GIT tag.
Definition: Jeep/JGITTags.hh:34
JLANG::make_predicate
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Definition: JPredicate.hh:128
JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20
JPHYSICS::getDeltaRaysFromMuon
double getDeltaRaysFromMuon(const double E, const JRange< double > T_GeV=JRange< double >(DELTARAY_TMIN, DELTARAY_TMAX))
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:260
JPHYSICS::is_deltarays
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:151
JPHYSICS::GNrad_t
static const JRadiationSource_t GNrad_t
Definition: JRadiation.hh:514
JPHYSICS::DENSITY_SEA_WATER
static const double DENSITY_SEA_WATER
Fixed environment values.
Definition: JPhysics/JConstants.hh:44
JPHYSICS::Brems_t
static const JRadiationSource_t Brems_t
Definition: JRadiation.hh:513
JPHYSICS::geanz
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
JPHYSICS::getThetaMCS
double getThetaMCS(const double E, const double x, const double X0, const double M, const double Q)
Get multiple Coulomb scattering angle.
Definition: JPhysicsSupportkit.hh:40
JPHYSICS::getIndexOfRefraction
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
Definition: JPhysics/JConstants.hh:134
JPHYSICS::gRock
static const JGeane_t gRock(2.67e-1 *0.9 *DENSITY_ROCK, 3.40e-4 *1.2 *DENSITY_ROCK)
Function object for energy loss of muon in rock.
JPHYSICS::is_scattered
bool is_scattered(const int pdf)
Test if given PDF type corresponds to scattered light.
Definition: JPDFTypes.hh:165
JPHYSICS::gWater
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
JPHYSICS::getTanThetaC
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Definition: JPhysics/JConstants.hh:156
JPHYSICS::EErad_t
static const JRadiationSource_t EErad_t
Definition: JRadiation.hh:512
JPHYSICS::JPDFType_t
JPDFType_t
PDF types.
Definition: JPDFTypes.hh:24
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER
@ SCATTERED_LIGHT_FROM_EMSHOWER
scattered light from EM shower
Definition: JPDFTypes.hh:38
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
@ SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition: JPDFTypes.hh:33
JPHYSICS::SCATTERED_LIGHT_FROM_MUON
@ SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition: JPDFTypes.hh:27
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER
@ DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition: JPDFTypes.hh:37
JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS
@ DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition: JPDFTypes.hh:32
JPHYSICS::DIRECT_LIGHT_FROM_MUON
@ DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition: JPDFTypes.hh:26
JPHYSICS::getDeltaRaysFromTau
double getDeltaRaysFromTau(const double E, const JRange< double > T_GeV=JRange< double >(DELTARAY_TMIN, DELTARAY_TMAX))
Equivalent EM-shower energy due to delta-rays per unit tau track length.
Definition: JPDFToolkit.hh:322
JPHYSICS::C
static const double C
Physics constants.
Definition: JPhysics/JConstants.hh:25
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
JSIRENE::getHitType
JHitType_t getHitType(const JPDFType_t pdf, const bool shower=false)
Get hit type corresponding to given PDF type.
Definition: JSireneToolkit.hh:151
JSIRENE::ns
std::vector< size_t > ns
Definition: JSireneToolkit.hh:138
JSIRENE::getNumberOfPhotoElectrons
const struct JSIRENE::@64 getNumberOfPhotoElectrons
Auxiliary data structure for determination of number of photo-electrons.
JSIRENE::pythia
static const JPythia pythia
Function object for relative light yield as a function of GEANT particle code.
Definition: JPythia.hh:96
JSYSTEM::getDate
const char * getDate()
Get current local date conform ISO-8601 standard.
Definition: JDateAndTime.hh:36
JTOOLS::accumulate
void accumulate(T &value, JBool< false > option)
Accumulation of value.
Definition: JToolsToolkit.hh:321
KM3NETDAQ::range_type
JRange< int > range_type
Definition: JDQSimulator.cc:160
create_update_issues.origin
string origin
Definition: create_update_issues.py:4
mc_usr_keys::energy_lost_in_can
const char *const energy_lost_in_can
Definition: io_ascii.hh:46
std
Definition: JSTDTypes.hh:14
pythia.hh
This file contains converted Fortran code from km3.
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:21
Evt::mc_hits
std::vector< Hit > mc_hits
MC: list of MC truth hits.
Definition: Evt.hh:48
Evt::mc_trks
std::vector< Trk > mc_trks
MC: list of MC truth tracks.
Definition: Evt.hh:49
Head
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:65
JAANET::JPDB::getInstance
static const JPDB & getInstance()
Get particle data book.
Definition: JPDB.hh:131
JAANET::can::zmin
double zmin
Bottom [m].
Definition: JHead.hh:598
JAANET::can::zmax
double zmax
Top [m].
Definition: JHead.hh:599
JAANET::can::r
double r
Radius [m].
Definition: JHead.hh:600
JAANET::coord_origin
Coordinate origin.
Definition: JHead.hh:732
JAANET::detector
Detector file.
Definition: JHead.hh:227
JAANET::fixedcan::zmax
double zmax
Top [m].
Definition: JHead.hh:639
JAANET::fixedcan::radius
double radius
Radius [m].
Definition: JHead.hh:640
JAANET::fixedcan::zmin
double zmin
Bottom [m].
Definition: JHead.hh:638
JAANET::fixedcan::ycenter
double ycenter
y-center [m]
Definition: JHead.hh:637
JAANET::fixedcan::xcenter
double xcenter
x-center [m]
Definition: JHead.hh:636
JAANET::simul
Generator for simulation.
Definition: JHead.hh:528
JCALIBRATE::JParameters
Auxiliary class for PMT parameters including threshold.
Definition: JParameters.hh:21
JGEOMETRY3D::JTransformation3D::transform
JPosition3D transform(const JPosition3D &pos) const
Transform position.
Definition: JTransformation3D.hh:49
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:68
JSIRENE::JHit_t
Auxiliary class to set-up Hit.
Definition: JSirene.hh:58
JSIRENE::JHits_t
Auxiliary data structure for list of hits with hit merging capability.
Definition: JSirene.hh:139
JSIRENE::JHits_t::merge
void merge(const double Tmax_ns)
Merge hits on same PMT that are within given time window.
Definition: JSirene.hh:148
JSIRENE::JPoint::getE
double getE() const
Get muon energy.
Definition: JSireneToolkit.hh:222
JSIRENE::JPoint::getEs
double getEs() const
Get shower energy.
Definition: JSireneToolkit.hh:233
JSIRENE::JTrack
Muon trajectory.
Definition: JSireneToolkit.hh:249
JSIRENE::JTrack::getPosition
JVector2D getPosition(const double z) const
Get muon position at given position along trajectory.
Definition: JSireneToolkit.hh:290
JSIRENE::JTrack::getE
double getE(const double z) const
Get muon energy at given position along trajectory.
Definition: JSireneToolkit.hh:266
JSIRENE::JTrk_t
Auxiliary class to set-up Trk.
Definition: JSirene.hh:190
JSIRENE::JVertex
Vertex of energy loss of muon.
Definition: JSireneToolkit.hh:337
JSIRENE::JVertex::getRange
double getRange() const
Get visible range of muon using default ionisation energy loss.
Definition: JSireneToolkit.hh:368
JSIRENE::JVertex::step
JVertex & step(const double ds)
Step using default ionisation energy loss.
Definition: JSireneToolkit.hh:426
JSIRENE::JVertex::applyEloss
void applyEloss(const JVersor3Z &Ts, const double Es)
Apply shower energy loss.
Definition: JSireneToolkit.hh:410
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:45
JSUPPORT::JLimit::max
static counter_type max()
Get maximum counter value.
Definition: JLimit.hh:128
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:72
JSYSTEM::ls
Auxiliary data structure to list files in directory.
Definition: JFilesystem.hh:20
JTOOLS::JHermiteSplineFunction1D_t
Type definition of a spline interpolation method based on a JCollection with double result type.
Definition: JFunction1D_t.hh:94
JTOOLS::JMAPLIST
Auxiliary class for recursive map list generation.
Definition: JMapList.hh:109
JTOOLS::JPolint1FunctionalGridMap
Type definition of a 1st degree polynomial interpolation based on a JGridMap implementation.
Definition: JFunctionalMap_t.hh:93
JTOOLS::JPolint1FunctionalMap
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
Definition: JFunctionalMap_t.hh:57
Trk::id
int id
track identifier
Definition: Trk.hh:16
Vec
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition: Vec.hh:13
Vec::z
double z
Definition: Vec.hh:14
Vec::x
double x
Definition: Vec.hh:14
Vec::y
double y
Definition: Vec.hh:14
Generated by doxygen 1.9.1