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