Jpp 20.0.0-195-g190c9e876
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 "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/JTimer.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 JPHYSICS::WILDCARD;
The file names are obtained by replacing JPHYSICS::WILDCARD; with

More accuracy can be achieved by setting compile option RADITION but it will run slower.

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 220 of file JSirene.cc.

221{
222 using namespace std;
223 using namespace JPP;
224
225 typedef JMultipleFileScanner<calibration_types> JCalibration_t;
226
227 string fileDescriptor;
228 JMultipleFileScanner<JAAnetTypes_t> inputFile;
229 JFileRecorder <JTYPELIST<JAAnetTypes_t, JMetaTypes_t, JRootTypes_t>::typelist> outputFile;
230 JLimit_t& numberOfEvents = inputFile.getLimit();
231 string detectorFile;
232 JCalibration_t calibrationFile;
233 double Tmax_s;
234 JParameters parameters;
235 bool writeEMShowers;
236 size_t numberOfHits;
237 double factor;
238 JRandom seed;
239
240 try {
241
242 JProperties properties;
243
244 properties.insert(gmake_property(parameters.Ecut_GeV));
245 properties.insert(gmake_property(parameters.Emin_GeV));
246 properties.insert(gmake_property(parameters.Dmin_m));
247 properties.insert(gmake_property(parameters.Emax_GeV));
248 properties.insert(gmake_property(parameters.Dmax_m));
249 properties.insert(gmake_property(parameters.Tmax_ns));
250 properties.insert(gmake_property(parameters.Nmax_NPE));
251 properties.insert(gmake_property(parameters.Nmax_PMT));
252 properties.insert(gmake_property(parameters.Tmin_GeV));
253
254 properties.insert(gmake_property(numberOfBins));
255 properties.insert(gmake_property(safetyFactor));
256
257 JParser<> zap("Main program to simulate detector response to muons and showers.");
258
259 zap['@'] = make_field(properties) = JPARSER::initialised();
260 zap['F'] = make_field(fileDescriptor, "file name descriptor for CDF tables");
261 zap['f'] = make_field(inputFile) = JPARSER::initialised();
262 zap['o'] = make_field(outputFile) = "sirene.root";
263 zap['n'] = make_field(numberOfEvents) = JLimit::max();
264 zap['a'] = make_field(detectorFile) = "";
265 zap['+'] = make_field(calibrationFile, "dynamical calibration data") = JPARSER::initialised();
266 zap['T'] = make_field(Tmax_s, "dynamical update time [s]") = 100.0;
267 zap['s'] = make_field(writeEMShowers, "store generated EM showers in event");
268 zap['N'] = make_field(numberOfHits, "minimum number of hits to output event") = 1;
269 zap['U'] = make_field(factor, "scaling factor applied to light yields") = 1.0;
270 zap['S'] = make_field(seed) = 0;
271 zap['d'] = make_field(debug) = 1;
272
273 zap(argc, argv);
274 }
275 catch(const exception &error) {
276 FATAL(error.what() << endl);
277 }
278
279
280 seed.set(gRandom);
281
282
283 const JMeta meta(argc, argv);
284
285 const double Zbed = 0.0; // level of seabed [m]
286
287 vector<JCDF_t> CDF;
288 vector<JCDG_t> CDG;
289
290 if (fileDescriptor != "") {
291 CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_MUON));
292 CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_MUON));
293 CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_DELTARAYS));
294 CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_DELTARAYS));
295
296 CDG.push_back(JCDG_t(fileDescriptor, DIRECT_LIGHT_FROM_EMSHOWER));
297 CDG.push_back(JCDG_t(fileDescriptor, SCATTERED_LIGHT_FROM_EMSHOWER));
298 }
299
300 double maximal_road_width = 0.0; // road width [m]
301 double maximal_distance = 0.0; // road width [m]
302
303 for (size_t i = 0; i != CDF.size(); ++i) {
304
305 DEBUG("Range CDF["<< CDF[i].type << "] " << CDF[i].function.intensity.getXmax() << " m" << endl);
306
307 maximal_road_width = max(maximal_road_width, CDF[i].function.intensity.getXmax());
308 }
309
310 for (size_t i = 0; i != CDG.size(); ++i) {
311
312 DEBUG("Range CDG["<< CDG[i].type << "] " << CDG[i].function.intensity.getXmax() << " m" << endl);
313
314 if (!is_scattered(CDF[i].type)) {
315 maximal_road_width = max(maximal_road_width, CDG[i].function.intensity.getXmax());
316 }
317
318 maximal_distance = max(maximal_distance, CDG[i].function.intensity.getXmax());
319 }
320
321 NOTICE("Maximal road width [m] " << maximal_road_width << endl);
322 NOTICE("Maximal distance [m] " << maximal_distance << endl);
323
324
325 if (detectorFile == "" || inputFile.empty()) {
326 STATUS("Nothing to be done." << endl);
327 return 0;
328 }
329
330 struct {
331 JDetector detector;
332 } master;
333
334 try {
335
336 STATUS("Load detector... " << flush);
337
338 load(detectorFile, master.detector);
339
340 STATUS("OK" << endl);
341 }
342 catch(const JException& error) {
343 FATAL(error);
344 }
345
346 unique_ptr<JDynamics> dynamics;
347
348 if (!calibrationFile.empty()) {
349
350 try {
351
352 dynamics.reset(new JDynamics(master.detector, Tmax_s));
353
354 dynamics->load(calibrationFile);
355 }
356 catch(const exception& error) {
357 FATAL(error.what());
358 }
359 }
360
361 const JDetector& detector = (dynamics ? dynamics->getDetector() : master.detector);
362
363 vector< shared_ptr<JRadiationInterface> > radiation;
364 vector< shared_ptr<JACoeffInterface> > ionization;
365
366 if (true) {
367
368 STATUS("Setting up radiation tables... " << flush);
369
370 const JRadiation hydrogen (JSeaWater::H .Z, JSeaWater::H .A, 40, 0.01, 0.1, 0.1);
371 const JRadiation oxygen (JSeaWater::O .Z, JSeaWater::O .A, 40, 0.01, 0.1, 0.1);
372 const JRadiation chlorine (JSeaWater::Cl.Z, JSeaWater::Cl.A, 40, 0.01, 0.1, 0.1);
373 const JRadiation sodium (JSeaWater::Na.Z, JSeaWater::Na.A, 40, 0.01, 0.1, 0.1);
374#ifdef RADIATION
375 const JRadiation calcium (JSeaWater::Ca.Z, JSeaWater::Ca.A, 40, 0.01, 0.1, 0.1);
376 const JRadiation magnesium(JSeaWater::Mg.Z, JSeaWater::Mg.A, 40, 0.01, 0.1, 0.1);
377 const JRadiation potassium(JSeaWater::K .Z, JSeaWater::K .A, 40, 0.01, 0.1, 0.1);
378 const JRadiation sulphur (JSeaWater::S .Z, JSeaWater::S .A, 40, 0.01, 0.1, 0.1);
379#endif
380
381 shared_ptr<JRadiation> Hydrogen (make_shared<JRadiationFunction>(hydrogen, 300, 0.2, 1.0e11));
382 shared_ptr<JRadiation> Oxygen (make_shared<JRadiationFunction>(oxygen, 300, 0.2, 1.0e11));
383 shared_ptr<JRadiation> Chlorine (make_shared<JRadiationFunction>(chlorine, 300, 0.2, 1.0e11));
384 shared_ptr<JRadiation> Sodium (make_shared<JRadiationFunction>(sodium, 300, 0.2, 1.0e11));
385#ifdef RADIATION
386 shared_ptr<JRadiation> Calcium (make_shared<JRadiationFunction>(calcium, 300, 0.2, 1.0e11));
387 shared_ptr<JRadiation> Magnesium(make_shared<JRadiationFunction>(magnesium,300, 0.2, 1.0e11));
388 shared_ptr<JRadiation> Potassium(make_shared<JRadiationFunction>(potassium,300, 0.2, 1.0e11));
389 shared_ptr<JRadiation> Sulphur (make_shared<JRadiationFunction>(sulphur, 300, 0.2, 1.0e11));
390#endif
391
392 radiation.push_back(make_shared<JRadiationSource>(11, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), JRadiation::EErad_t));
393 radiation.push_back(make_shared<JRadiationSource>(12, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), JRadiation::EErad_t));
394 radiation.push_back(make_shared<JRadiationSource>(13, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), JRadiation::EErad_t));
395 radiation.push_back(make_shared<JRadiationSource>(14, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), JRadiation::EErad_t));
396#ifdef RADIATION
397 radiation.push_back(make_shared<JRadiationSource>(15, Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), JRadiation::EErad_t));
398 radiation.push_back(make_shared<JRadiationSource>(16, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), JRadiation::EErad_t));
399 radiation.push_back(make_shared<JRadiationSource>(17, Potassium,DENSITY_SEA_WATER * JSeaWater::K(), JRadiation::EErad_t));
400 radiation.push_back(make_shared<JRadiationSource>(18, Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), JRadiation::EErad_t));
401#endif
402
403 radiation.push_back(make_shared<JRadiationSource>(21, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), JRadiation::Brems_t));
404 radiation.push_back(make_shared<JRadiationSource>(22, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), JRadiation::Brems_t));
405 radiation.push_back(make_shared<JRadiationSource>(23, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), JRadiation::Brems_t));
406 radiation.push_back(make_shared<JRadiationSource>(24, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), JRadiation::Brems_t));
407#ifdef RADIATION
408 radiation.push_back(make_shared<JRadiationSource>(25, Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), JRadiation::Brems_t));
409 radiation.push_back(make_shared<JRadiationSource>(26, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), JRadiation::Brems_t));
410 radiation.push_back(make_shared<JRadiationSource>(27, Potassium,DENSITY_SEA_WATER * JSeaWater::K(), JRadiation::Brems_t));
411 radiation.push_back(make_shared<JRadiationSource>(28, Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), JRadiation::Brems_t));
412#endif
413
414 radiation.push_back(make_shared<JRadiationSource>(31, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), JRadiation::GNrad_t));
415 radiation.push_back(make_shared<JRadiationSource>(32, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), JRadiation::GNrad_t));
416 radiation.push_back(make_shared<JRadiationSource>(33, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), JRadiation::GNrad_t));
417 radiation.push_back(make_shared<JRadiationSource>(34, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), JRadiation::GNrad_t));
418#ifdef RADIATION
419 radiation.push_back(make_shared<JRadiationSource>(35, Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), JRadiation::GNrad_t));
420 radiation.push_back(make_shared<JRadiationSource>(36, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), JRadiation::GNrad_t));
421 radiation.push_back(make_shared<JRadiationSource>(37, Potassium,DENSITY_SEA_WATER * JSeaWater::K(), JRadiation::GNrad_t));
422 radiation.push_back(make_shared<JRadiationSource>(38, Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), JRadiation::GNrad_t));
423#endif
424
425 radiation.push_back(make_shared<JDISSource>(100, DENSITY_SEA_WATER));
426
427 radiation.push_back(make_shared<JDeltaRaysSource>(200, DENSITY_SEA_WATER, parameters.Tmin_GeV));
428
429 ionization.push_back(make_shared<JACoeffSource>(Oxygen, DENSITY_SEA_WATER * JSeaWater::O()));
430 ionization.push_back(make_shared<JACoeffSource>(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl()));
431 ionization.push_back(make_shared<JACoeffSource>(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H()));
432 ionization.push_back(make_shared<JACoeffSource>(Sodium, DENSITY_SEA_WATER * JSeaWater::Na()));
433#ifdef RADIATION
434 ionization.push_back(make_shared<JACoeffSource>(Calcium, DENSITY_SEA_WATER * JSeaWater::Ca()));
435 ionization.push_back(make_shared<JACoeffSource>(Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg()));
436 ionization.push_back(make_shared<JACoeffSource>(Potassium,DENSITY_SEA_WATER * JSeaWater::K()));
437 ionization.push_back(make_shared<JACoeffSource>(Sulphur, DENSITY_SEA_WATER * JSeaWater::S()));
438#endif
439
440 STATUS("OK" << endl);
441 }
442
443
444 JCylinder3D cylinder(detector.begin(), detector.end());
445
446 cylinder.addMargin(maximal_distance);
447
448 if (cylinder.getZmin() < Zbed) {
449 cylinder.setZmin(Zbed);
450 }
451
452 NOTICE("Light generation volume: " << cylinder << endl);
453
454
455 Vec offset(0,0,0);
456 Head header;
457
458 try {
459
460 header = inputFile.getHeader();
461
462 JAANET::JHead buffer(header);
463
464 buffer.simul.push_back(JAANET::simul());
465
466 buffer.simul.rbegin()->program = APPLICATION_JSIRENE;
467 buffer.simul.rbegin()->version = getGITVersion();
468 buffer.simul.rbegin()->date = getDate();
469 buffer.simul.rbegin()->time = getTime();
470
471 buffer.push(&JAANET::JHead::simul);
472
473 buffer.detector.push_back(JAANET::detector());
474
475 buffer.detector.rbegin()->program = APPLICATION_JSIRENE;
476 buffer.detector.rbegin()->filename = detectorFile;
477
478 buffer.push(&JAANET::JHead::detector);
479
480 offset += Vec(cylinder.getX(), cylinder.getY(), 0.0);
481 offset -= getOrigin(buffer);
482
483 if (buffer.is_valid(&JAANET::JHead::fixedcan)) {
484
485 buffer.fixedcan.xcenter += offset.x;
486 buffer.fixedcan.ycenter += offset.y;
487 buffer.fixedcan.zmin += offset.z;
488 buffer.fixedcan.zmax += offset.z;
489
490 } else {
491
492 buffer.fixedcan.xcenter = cylinder.getX();
493 buffer.fixedcan.ycenter = cylinder.getY();
494
495 if (buffer.is_valid(&JAANET::JHead::can)) {
496
497 buffer.fixedcan.radius = buffer.can.r;
498 buffer.fixedcan.zmin = buffer.can.zmin + offset.z;
499 buffer.fixedcan.zmax = buffer.can.zmax + offset.z;
500 } else {
501
502 buffer.fixedcan.radius = cylinder.getRadius();
503 buffer.fixedcan.zmin = cylinder.getZmin();
504 buffer.fixedcan.zmax = cylinder.getZmax();
505 }
506 }
507
508 buffer.push(&JAANET::JHead::fixedcan);
509
510 if (buffer.is_valid(&JAANET::JHead::coord_origin)) {
511
512 buffer.coord_origin = coord_origin(0.0, 0.0, 0.0);
513
514 buffer.push(&JAANET::JHead::coord_origin);
515 }
516
517 copy(buffer, header);
518 }
519 catch(const JException& error) {
520 FATAL(error);
521 }
522
523 NOTICE("Offset applied to true tracks is: " << offset << endl);
524
525 TH1D job("job", NULL, 400, 0.5, 400.5);
526 TProfile cpu("cpu", NULL, 16, 0.0, 8.0);
527 TProfile2D rms("rms", NULL, 16, 0.0, 8.0, 251, -0.5, 250.5);
528 TProfile2D rad("rad", NULL, 16, 0.0, 8.0, 251, -0.5, 250.5);
529
530
531 outputFile.open();
532
533 if (!outputFile.is_open()) {
534 FATAL("Error opening file " << outputFile << endl);
535 }
536
537 outputFile.put(meta);
538 outputFile.put(header);
539 outputFile.put(*gRandom);
540
541 const double epsilon = 1.0e-6; // precision angle [rad]
542 const JRange<double> pi(epsilon, PI - epsilon); // constrain angle
543
544 JTimer timer;
545
546 for (JMultipleFileScanner<Evt>& in = inputFile; in.hasNext(); ) {
547
548 STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
549
550 job.Fill(1.0);
551
552 Evt* evt = in.next();
553
554 for (vector<Trk>::iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
555 track->pos += offset;
556 }
557
558 Evt event(*evt); // output
559
560 event.mc_hits.clear();
561
562 if (dynamics) {
563
564 if (event.mc_event_time == TTimeStamp(0)) {
565 FATAL("Monte Carlo event time undefined." << endl);
566 }
567
568 dynamics->update(event.mc_event_time.AsDouble());
569 }
570
571
572 JHits_t mc_hits; // temporary buffer
573
574 timer.reset();
575 timer.start();
576
577 for (vector<Trk>::const_iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
578
579 if (!track->is_finalstate()) {
580 continue; // only final state particles produce light
581 }
582
583 if (is_muon(*track)) {
584
585 // -----------------------------------------------
586 // muon
587 // -----------------------------------------------
588
589 job.Fill(2.0);
590
591 typedef JDetectorSubset<JModuleComparator3Z_t> JDetectorSubset_t;
592
593 const JCylinder3D::intersection_type intersection = cylinder.getIntersection(getAxis(*track));
594
595 double Zmin = intersection.first;
596 double Zmax = intersection.second;
597
598 if (Zmax - Zmin <= parameters.Dmin_m) {
599 continue;
600 }
601
602 JVertex vertex(0.0, track->t, track->E); // start of muon
603
604 if (track->pos.z < Zbed) { // propagate muon through rock
605
606 if (track->dir.z > 0.0)
607 vertex.step(gRock, (Zbed - track->pos.z) / track->dir.z);
608 else
609 continue;
610 }
611
612 if (vertex.getZ() < Zmin) { // propagate muon through water
613 vertex.step(gWater, Zmin - vertex.getZ());
614 }
615
616 if (vertex.getRange() <= parameters.Dmin_m) {
617 continue;
618 }
619
620 job.Fill(3.0);
621
622 const JDetectorSubset_t subdetector(detector, getAxis(*track), maximal_road_width);
623
624 if (subdetector.empty()) {
625 continue;
626 }
627
628 job.Fill(4.0);
629
630 JTrack muon(vertex); // propagate muon trough detector
631
632 while (vertex.getE() >= parameters.Emin_GeV && vertex.getZ() < Zmax) {
633
634 const int N = radiation.size();
635
636 double li[N]; // inverse interaction lengths
637 double ls = 1.0e-5; // minimal total inverse interaction length [m^-1]
638
639 for (int i = 0; i != N; ++i) {
640 ls += li[i] = radiation[i]->getInverseInteractionLength(vertex.getE());
641 }
642
643 // ionization energy loss without stochastic energy losses due to delta-rays
644
645 double As = -JDeltaRays::getEnergyLossFromMuon(vertex.getE(), JEnergyRange(parameters.Tmin_GeV, JDeltaRays::TMAX_GEV));
646
647 for (size_t i = 0; i != ionization.size(); ++i) {
648 As += ionization[i]->getA(vertex.getE());
649 }
650
651 double step = gRandom->Exp(1.0) / ls; // distance to next radiation process
652 double range = vertex.getRange(As); // range of muon
653
654 if (vertex.getE() < parameters.Emax_GeV) { // limited step size
655 if (parameters.Dmax_m < range) {
656 range = parameters.Dmax_m;
657 }
658 }
659
660 double ts = getThetaMCS(vertex.getE(), min(step,range)); // multiple Coulomb scattering angle [rad]
661 double T2 = ts*ts; //
662
663 rms.Fill(log10(vertex.getE()), (Double_t) 0, ts*ts);
664
665 vertex.getDirection() += getRandomDirection(T2/3.0); // multiple Coulomb planar scattering
666
667 vertex.step(As, min(step,range)); // ionization energy loss
668
669 double Es = 0.0; // shower energy [GeV]
670
671 if (step < range) {
672
673 if (vertex.getE() >= parameters.Emin_GeV) {
674
675 double y = gRandom->Uniform(ls);
676
677 for (int i = 0; i != N; ++i) {
678
679 y -= li[i];
680
681 if (y < 0.0) {
682
683 Es = radiation[i]->getEnergyOfShower(vertex.getE()); // shower energy [GeV]
684 ts = radiation[i]->getThetaRMS(vertex.getE(), Es); // scattering angle [rad]
685
686 T2 += ts*ts;
687
688 rms.Fill(log10(vertex.getE()), (Double_t) radiation[i]->getID(), ts*ts);
689 rad.Fill(log10(vertex.getE()), (Double_t) radiation[i]->getID(), Es);
690
691 break;
692 }
693 }
694 }
695 }
696
697 vertex.applyEloss(getRandomDirection(T2), Es);
698
699 muon.push_back(vertex);
700
701 vertex.reset(); // reset after being added to muon
702 }
703
704 // add muon end point
705
706 if (vertex.getZ() < Zmax && vertex.getRange() > 0.0) {
707
708 vertex.step(vertex.getRange());
709
710 muon.push_back(vertex);
711 }
712
713 // add information to output muon
714
715 vector<Trk>::iterator trk = find_if(event.mc_trks.begin(),
716 event.mc_trks.end(),
717 make_predicate(&Trk::id, track->id));
718
719 if (trk != event.mc_trks.end()) {
720 trk->len = (muon.rbegin()->getZ() < Zmax ? +1 : -1) * (muon.rbegin()->getZ() - muon.begin()->getZ());
721 trk->setusr(mc_usr_keys::energy_lost_in_can, muon.begin()->getE() - muon.rbegin()->getE());
722 }
723
724 for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
725
726 if (module->empty()) {
727 continue;
728 }
729
730 const double z0 = muon.begin()->getZ();
731 const double t0 = muon.begin()->getT();
732 const double Z = module->getZ() - module->getX() / getTanThetaC();
733
734 if (Z >= muon.begin()->getZ() && Z <= muon.rbegin()->getZ()) {
735
736 const JVector2D pos = muon.getPosition(Z);
737 const double R = hypot(module->getX() - pos.getX(),
738 module->getY() - pos.getY());
739
740 for (size_t i = 0; i != CDF.size(); ++i) {
741
742 if (R < CDF[i].integral.getXmax()) {
743
744 try {
745
746 double W = 1.0; // mip
747
748 if (is_deltarays(CDF[i].type)) { // delta-rays
749 W = JDeltaRays::getEnergyLossFromMuon(muon.getE(Z), JEnergyRange(JDeltaRays::getTmin(), parameters.Tmin_GeV));
750 }
751
752 const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
753 const size_t N = getPoisson(NPE);
754
755 if (N != 0) {
756
757 vector<double> npe;
758
759 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
760
761 const double R = hypot(pmt->getX() - pos.getX(),
762 pmt->getY() - pos.getY());
763 const double theta = pi.constrain(pmt->getTheta());
764 const double phi = pi.constrain(fabs(pmt->getPhi()));
765
766 npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
767 }
768
769 const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
770
771 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
772
773 const double R = hypot(pmt->getX() - pos.getX(),
774 pmt->getY() - pos.getY());
775 const double Z = pmt->getZ() - z0;
776 const double theta = pi.constrain(pmt->getTheta());
777 const double phi = pi.constrain(fabs(pmt->getPhi()));
778
779 size_t n0 = min(ns[distance(module->begin(),pmt)], parameters.Nmax_PMT);
780
781 job.Fill((double) (100 + CDF[i].type), (double) n0);
782
783 while (n0 != 0) {
784
785 const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
786 const int n1 = getNumberOfPhotoElectrons(n0);
787
788 mc_hits.push_back(JHit_t(mc_hits.size() + 1,
789 pmt->getID(),
790 getHitType(CDF[i].type),
791 track->id,
792 t0 + (R * getTanThetaC() + Z) / C + t1,
793 n1));
794
795 n0 -= n1;
796 }
797 }
798
799 if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
800 job.Fill((double) (300 + CDF[i].type));
801 }
802 }
803 }
804 catch(const exception& error) {
805 job.Fill((double) (200 + CDF[i].type));
806 }
807 }
808 }
809 }
810 }
811
812 for (JTrack::const_iterator vertex = muon.begin(); vertex != muon.end(); ++vertex) {
813
814 const double Es = vertex->getEs();
815
816 if (Es >= parameters.Ecut_GeV) {
817
818 const double z0 = vertex->getZ();
819 const double t0 = vertex->getT();
820 const double DZ = geanz.getMaximum(Es);
821
822 int origin = track->id;
823
824 if (writeEMShowers) {
825 origin = event.mc_trks.size() + 1;
826 }
827
828 int number_of_hits = 0;
829
830 JDetectorSubset_t::range_type range = subdetector.getRange(z0 - maximal_distance,
831 z0 + maximal_distance);
832
833 for (JDetector::const_iterator module = range.begin(); module != range.end(); ++module) {
834
835 if (module->empty()) {
836 continue;
837 }
838
839 const double R = hypot(module->getX() - vertex->getX(),
840 module->getY() - vertex->getY());
841 const double Z = module->getZ() - z0 - DZ;
842 const double D = sqrt(R*R + Z*Z);
843 const double cd = Z / D;
844
845 for (size_t i = 0; i != CDG.size(); ++i) {
846
847 if (D < CDG[i].integral.getXmax()) {
848
849 try {
850
851 const double NPE = CDG[i].integral.getNPE(D, cd) * Es * module->size() * factor;
852 const size_t N = getPoisson(NPE);
853
854 if (N != 0) {
855
856 vector<double> npe;
857
858 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
859
860 const double R = hypot(pmt->getX() - vertex->getX(),
861 pmt->getY() - vertex->getY());
862 const double Z = pmt->getZ() - z0 - DZ;
863 const double D = sqrt(R*R + Z*Z);
864 const double cd = Z / D;
865 const double theta = pi.constrain(pmt->getTheta());
866 const double phi = pi.constrain(fabs(pmt->getPhi()));
867
868 npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor);
869 }
870
871 const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
872
873 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
874
875 const double R = hypot(pmt->getX() - vertex->getX(),
876 pmt->getY() - vertex->getY());
877 const double theta = pi.constrain(pmt->getTheta());
878 const double phi = pi.constrain(fabs(pmt->getPhi()));
879
880 size_t n0 = min(ns[distance(module->begin(),pmt)], parameters.Nmax_PMT);
881
882 job.Fill((double) (100 + CDG[i].type), (double) n0);
883
884 while (n0 != 0) {
885
886 const double dz = geanz.getLength(Es, gRandom->Rndm());
887 const double Z = pmt->getZ() - z0 - dz;
888 const double D = sqrt(R*R + Z*Z);
889 const double cd = Z / D;
890
891 const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
892 const int n1 = getNumberOfPhotoElectrons(n0);
893
894 mc_hits.push_back(JHit_t(mc_hits.size() + 1,
895 pmt->getID(),
896 getHitType(CDG[i].type),
897 origin,
898 t0 + (dz + D * getIndexOfRefraction()) / C + t1,
899 n1));
900
901 n0 -= n1;
902
903 number_of_hits += n1;
904 }
905 }
906
907 if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
908 job.Fill((double) (300 + CDG[i].type));
909 }
910 }
911 }
912 catch(const exception& error) {
913 job.Fill((double) (200 + CDG[i].type));
914 }
915 }
916 }
917 }
918
919 if (writeEMShowers && number_of_hits != 0) {
920
921 event.mc_trks.push_back(JTrk_t(origin,
922 TRACK_TYPE_PHOTON,
923 track->id,
924 track->pos + track->dir * vertex->getZ(),
925 track->dir,
926 vertex->getT(),
927 Es));
928 }
929 }
930 }
931
932 } else if (track->len > 0.0) {
933
934 // -----------------------------------------------
935 // decayed particles treated as mip (tau includes mip+deltaray)
936 // -----------------------------------------------
937
938 job.Fill(6.0);
939
940 const double z0 = 0.0;
941 const double z1 = z0 + track->len;
942 const double t0 = track->t;
943 const double E = track->E;
944
945 const JTransformation3D transformation = getTransformation(*track);
946
947 JModule buffer;
948
949 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
950
951 if (module->empty()) {
952 continue;
953 }
954
955 const JPosition3D pos = transformation.transform(module->getPosition());
956
957 const double R = pos.getX();
958 const double Z = pos.getZ() - R / getTanThetaC();
959
960 if (Z < z0 ||
961 Z > z1 ||
962 R > maximal_road_width) {
963 continue;
964 }
965
966 for (size_t i = 0; i != CDF.size(); ++i) {
967
968 double W = 1.0; // mip
969
970 if (is_deltarays(CDF[i].type)) {
971
972 if (is_tau(*track))
973 W = JDeltaRays::getEnergyLossFromTau(E); // delta-rays
974 else
975 continue;
976 }
977
978 if (R < CDF[i].integral.getXmax()) {
979
980 try {
981
982 const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
983 const size_t N = getPoisson(NPE);
984
985 if (N != 0) {
986
987 buffer = *module;
988
989 buffer.transform(transformation);
990
991 vector<double> npe;
992
993 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
994
995 const double R = pmt->getX();
996 const double theta = pi.constrain(pmt->getTheta());
997 const double phi = pi.constrain(fabs(pmt->getPhi()));
998
999 npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
1000 }
1001
1002 const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
1003
1004 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1005
1006 const double R = pmt->getX();
1007 const double Z = pmt->getZ() - z0;
1008 const double theta = pi.constrain(pmt->getTheta());
1009 const double phi = pi.constrain(fabs(pmt->getPhi()));
1010
1011 size_t n0 = min(ns[distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
1012
1013 job.Fill((double) (120 + CDF[i].type), (double) n0);
1014
1015 while (n0 != 0) {
1016
1017 const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
1018 const int n1 = getNumberOfPhotoElectrons(n0);
1019
1020 mc_hits.push_back(JHit_t(mc_hits.size() + 1,
1021 pmt->getID(),
1022 getHitType(CDF[i].type),
1023 track->id,
1024 t0 + (R * getTanThetaC() + Z) / C + t1,
1025 n1));
1026
1027 n0 -= n1;
1028 }
1029 }
1030
1031 if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
1032 job.Fill((double) (320 + CDF[i].type));
1033 }
1034 }
1035 }
1036 catch(const exception& error) {
1037 job.Fill((double) (220 + CDF[i].type));
1038 }
1039 }
1040 }
1041 }
1042
1043 if (!buffer.empty()) {
1044 job.Fill(7.0);
1045 }
1046
1047 } else if (!is_neutrino(*track)) {
1048
1049 if (JPDB::getInstance().hasPDG(track->type)) {
1050
1051 // -----------------------------------------------
1052 // electron or hadron
1053 // -----------------------------------------------
1054
1055 job.Fill(8.0);
1056
1057 double E = track->E;
1058
1059 try {
1060 E = getKineticEnergy(E, JPDB::getInstance().getPDG(track->type).mass);
1061 }
1062 catch(const exception& error) {
1063 ERROR(error.what() << endl);
1064 }
1065
1066 E = pythia(track->type, E);
1067
1068 if (E >= parameters.Ecut_GeV && cylinder.getDistance(getPosition(*track)) < parameters.Dmin_m) {
1069
1070 const double z0 = 0.0;
1071 const double t0 = track->t;
1072 const double DZ = geanz.getMaximum(E);
1073
1074 const JTransformation3D transformation = getTransformation(*track);
1075
1076 JModule buffer;
1077
1078 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
1079
1080 if (module->empty()) {
1081 continue;
1082 }
1083
1084 const JPosition3D pos = transformation.transform(module->getPosition());
1085
1086 const double R = pos.getX();
1087 const double Z = pos.getZ() - z0 - DZ;
1088 const double D = sqrt(R*R + Z*Z);
1089 const double cd = Z / D;
1090
1091 for (size_t i = 0; i != CDG.size(); ++i) {
1092
1093 if (D < CDG[i].integral.getXmax()) {
1094
1095 try {
1096
1097 const double NPE = CDG[i].integral.getNPE(D, cd) * E * module->size() * factor;
1098 const size_t N = getPoisson(NPE);
1099
1100 if (N != 0) {
1101
1102 buffer = *module;
1103
1104 buffer.transform(transformation);
1105
1106 vector<double> npe;
1107
1108 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1109
1110 const double R = pmt->getX();
1111 const double Z = pmt->getZ() - z0 - DZ;
1112 const double D = sqrt(R*R + Z*Z);
1113 const double cd = Z / D;
1114 const double theta = pi.constrain(pmt->getTheta());
1115 const double phi = pi.constrain(fabs(pmt->getPhi()));
1116
1117 npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * E * factor);
1118 }
1119
1120 const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
1121
1122 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1123
1124 const double theta = pi.constrain(pmt->getTheta());
1125 const double phi = pi.constrain(fabs(pmt->getPhi()));
1126
1127 size_t n0 = min(ns[distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
1128
1129 job.Fill((double) (140 + CDG[i].type), (double) n0);
1130
1131 while (n0 != 0) {
1132
1133 const double dz = geanz.getLength(E, gRandom->Rndm());
1134 const double Z = pmt->getZ() - z0 - dz;
1135 const double D = sqrt(R*R + Z*Z);
1136 const double cd = Z / D;
1137
1138 const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
1139 const int n1 = getNumberOfPhotoElectrons(n0);
1140
1141 mc_hits.push_back(JHit_t(mc_hits.size() + 1,
1142 pmt->getID(),
1143 getHitType(CDG[i].type, true),
1144 track->id,
1145 t0 + (dz + D * getIndexOfRefraction()) / C + t1,
1146 n1));
1147
1148 n0 -= n1;
1149 }
1150 }
1151
1152 if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
1153 job.Fill((double) (340 + CDG[i].type));
1154 }
1155 }
1156 }
1157 catch(const exception& error) {
1158 job.Fill((double) (240 + CDG[i].type));
1159 }
1160 }
1161 }
1162 }
1163
1164 if (!buffer.empty()) {
1165 job.Fill(9.0);
1166 }
1167
1168 } else {
1169 job.Fill(21.0);
1170 }
1171 }
1172 }
1173 }
1174
1175 if (!mc_hits.empty()) {
1176
1177 mc_hits.merge(parameters.Tmax_ns);
1178
1179 event.mc_hits.resize(mc_hits.size());
1180
1181 copy(mc_hits.begin(), mc_hits.end(), event.mc_hits.begin());
1182 }
1183
1184 timer.stop();
1185
1186 if (has_neutrino(event)) {
1187 cpu.Fill(log10(get_neutrino(event).E), (double) timer.usec_ucpu * 1.0e-3);
1188 }
1189
1190 if (event.mc_hits.size() >= numberOfHits) {
1191
1192 outputFile.put(event);
1193
1194 job.Fill(10.0);
1195 }
1196 }
1197 STATUS(endl);
1198
1199 outputFile.put(job);
1200 outputFile.put(cpu);
1201 outputFile.put(rms);
1202 outputFile.put(rad);
1203 outputFile.put(*gRandom);
1204
1205 JMultipleFileScanner<JTYPELIST<JMetaTypes_t, MultiHead>::typelist> io(inputFile);
1206
1207 io >> outputFile;
1208
1209 outputFile.close();
1210}
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:2142
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
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::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::JException
General exception.
Definition JException.hh:24
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
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:68
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