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