Jpp 20.0.0-72-g597b30bc9
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 "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 218 of file JSirene.cc.

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

◆ numberOfBins

int numberOfBins = 200

number of bins for average CDF integral of optical module

Definition at line 73 of file JSirene.cc.

◆ safetyFactor

double safetyFactor = 1.7

safety factor for average CDF integral of optical module

Definition at line 74 of file JSirene.cc.

Generated by doxygen 1.12.0