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