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