Jpp 21.0.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
Functions | Variables
JSirene.cc File Reference

Main program to simulate detector response to muons and showers. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <limits>
#include <numeric>
#include <memory>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "TRandom3.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/MultiHead.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/definitions/applications.hh"
#include "km3net-dataformat/definitions/trkmembers.hh"
#include "JLang/Jpp.hh"
#include "JLang/JPredicate.hh"
#include "JLang/JClock.hh"
#include "JSystem/JDateAndTime.hh"
#include "JPhysics/JCDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JGeane.hh"
#include "JPhysics/JGeanz.hh"
#include "JPhysics/JRadiation.hh"
#include "JPhysics/JRadiationFunction.hh"
#include "JPhysics/JRadiationSource.hh"
#include "JPhysics/JACoeffSource.hh"
#include "JPhysics/JPhysicsToolkit.hh"
#include "JPhysics/JPhysicsSupportkit.hh"
#include "JPhysics/JSeaWater.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JPDB.hh"
#include "JSirene/JSirene.hh"
#include "JSirene/pythia.hh"
#include "JSirene/JCDFTable1D.hh"
#include "JSirene/JCDFTable2D.hh"
#include "JSirene/JSireneToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSubset.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "JDynamics/JDynamics.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JROOT/JRandom.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Variables

int debug
 debug level
 
int numberOfBins = 200
 number of bins for average CDF integral of optical module
 
double safetyFactor = 1.7
 safety factor for average CDF integral of optical module
 

Detailed Description

Main program to simulate detector response to muons and showers.

Picture by Claudine Colnard

Note that CDF file descriptor should contain the wild card character JEEP::FILENAME_WILDCARD;
The file names are obtained by replacing JEEP::FILENAME_WILDCARD; with

The list of chemical elements used for energy losses of the muon constitutes a compromise between accuracy and speed.
It can be increased setting compile option RADIATION.

The CDF tables can be produced with the script JMakePDF.sh: 

JMakePDF.sh -P

(option -h will print all available options). Note that the script will launch a large number of processes (JMakePDF and JMakePDG)
which may take a considerable amount of time to completion.
On a standard desktop, all jobs should be finished within 1/2 a day or so.

The same script should then be run with option -M to merge the PDF files, i.e: 

JMakePDF.sh -M

CDF tables are obtained by running the same script with option -C, i.e: 

JMakePDF.sh -C

The various PDFs can be drawn using the JDrawPDF or JDrawPDG applications.
The tabulated PDFs can be plotted using the JPlotPDF or JPlotPDG applications.
The tabulated CDFs can be plotted using the JPlotCDF or JPlotCDG applications.

Author
mdejong

Definition in file JSirene.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 221 of file JSirene.cc.

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}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
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
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2140
gmake_property
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Definition JProperties.hh:1139
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
debug
int debug
debug level
Definition JSirene.cc:74
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::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
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::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::JPosition3D
Data structure for position in three dimensions.
Definition JPosition3D.hh:38
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
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
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1697
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::JRange
Range of values.
Definition JRange.hh:42
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
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
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
create_update_issues.origin
str origin
Definition create_update_issues.py:4
mc_usr_keys::energy_lost_in_can
const char *const energy_lost_in_can
Definition io_ascii.hh:46
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
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::coord_origin
Coordinate origin.
Definition JHead.hh:732
JAANET::detector
Detector file.
Definition JHead.hh:227
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::JTrk_t
Auxiliary class to set-up Trk.
Definition JSirene.hh:192
JSIRENE::JVertex
Vertex of energy loss of muon.
Definition JSireneToolkit.hh:337
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
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

Variable Documentation

◆ debug

int debug

debug level

Definition at line 74 of file JSirene.cc.

◆ numberOfBins

int numberOfBins = 200

number of bins for average CDF integral of optical module

Definition at line 75 of file JSirene.cc.

◆ safetyFactor

double safetyFactor = 1.7

safety factor for average CDF integral of optical module

Definition at line 76 of file JSirene.cc.

Generated by doxygen 1.12.0