Jpp master_rocky-44-g75b7c4f75
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 "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 215 of file JSirene.cc.

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

◆ numberOfBins

int numberOfBins = 200

number of bins for average CDF integral of optical module

Definition at line 70 of file JSirene.cc.

◆ safetyFactor

double safetyFactor = 1.7

safety factor for average CDF integral of optical module

Definition at line 71 of file JSirene.cc.

Generated by doxygen 1.12.0