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