Jpp 20.0.0-rc.9-4-g539b3ce77
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 double As = 0.0; // ionization energy loss
616
617 for (size_t i = 0; i != ionization.size(); ++i) {
618 As += ionization[i]->getA(vertex.getE());
619 }
620
621 double step = gRandom->Exp(1.0) / ls; // distance to next radiation process
622 double range = vertex.getRange(As); // range of muon
623
624 if (vertex.getE() < parameters.Emax_GeV) { // limited step size
625 if (parameters.Dmax_m < range) {
626 range = parameters.Dmax_m;
627 }
628 }
629
630 double ts = getThetaMCS(vertex.getE(), min(step,range)); // multiple Coulomb scattering angle [rad]
631 double T2 = ts*ts; //
632
633 rms.Fill(log10(vertex.getE()), (Double_t) 0, ts*ts);
634
635 vertex.getDirection() += getRandomDirection(T2/3.0); // multiple Coulomb planar scattering
636
637 vertex.step(As, min(step,range)); // ionization energy loss
638
639 double Es = 0.0; // shower energy [GeV]
640
641 if (step < range) {
642
643 if (vertex.getE() >= parameters.Emin_GeV) {
644
645 double y = gRandom->Uniform(ls);
646
647 for (int i = 0; i != N; ++i) {
648
649 y -= li[i];
650
651 if (y < 0.0) {
652
653 Es = radiation[i]->getEnergyOfShower(vertex.getE()); // shower energy [GeV]
654 ts = radiation[i]->getThetaRMS(vertex.getE(), Es); // scattering angle [rad]
655
656 T2 += ts*ts;
657
658 rms.Fill(log10(vertex.getE()), (Double_t) radiation[i]->getID(), ts*ts);
659 rad.Fill(log10(vertex.getE()), (Double_t) radiation[i]->getID(), Es);
660
661 break;
662 }
663 }
664 }
665 }
666
667 vertex.applyEloss(getRandomDirection(T2), Es);
668
669 muon.push_back(vertex);
670
671 vertex.reset(); // reset after being added to muon
672 }
673
674 // add muon end point
675
676 if (vertex.getZ() < Zmax && vertex.getRange() > 0.0) {
677
678 vertex.step(vertex.getRange());
679
680 muon.push_back(vertex);
681 }
682
683 // add information to output muon
684
685 vector<Trk>::iterator trk = find_if(event.mc_trks.begin(),
686 event.mc_trks.end(),
687 make_predicate(&Trk::id, track->id));
688
689 if (trk != event.mc_trks.end()) {
690 trk->len = (muon.rbegin()->getZ() < Zmax ? +1 : -1) * (muon.rbegin()->getZ() - muon.begin()->getZ());
691 trk->setusr(mc_usr_keys::energy_lost_in_can, muon.begin()->getE() - muon.rbegin()->getE());
692 }
693
694 for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
695
696 const double z0 = muon.begin()->getZ();
697 const double t0 = muon.begin()->getT();
698 const double Z = module->getZ() - module->getX() / getTanThetaC();
699
700 if (Z >= muon.begin()->getZ() && Z <= muon.rbegin()->getZ()) {
701
702 const JVector2D pos = muon.getPosition(Z);
703 const double R = hypot(module->getX() - pos.getX(),
704 module->getY() - pos.getY());
705
706 for (size_t i = 0; i != CDF.size(); ++i) {
707
708 if (R < CDF[i].integral.getXmax()) {
709
710 try {
711
712 double W = 1.0; // mip
713
714 if (is_deltarays(CDF[i].type)) {
715 W = JDeltaRays::getEnergyLossFromMuon(muon.getE(Z), // delta-rays
716 JEnergyRange(JDeltaRays::getTmin(), parameters.Tmin_GeV));
717 }
718
719 const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
720 const size_t N = getPoisson(NPE);
721
722 if (N != 0) {
723
724 vector<double> npe;
725
726 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
727
728 const double R = hypot(pmt->getX() - pos.getX(),
729 pmt->getY() - pos.getY());
730 const double theta = pi.constrain(pmt->getTheta());
731 const double phi = pi.constrain(fabs(pmt->getPhi()));
732
733 npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
734 }
735
736 const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
737
738 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
739
740 const double R = hypot(pmt->getX() - pos.getX(),
741 pmt->getY() - pos.getY());
742 const double Z = pmt->getZ() - z0;
743 const double theta = pi.constrain(pmt->getTheta());
744 const double phi = pi.constrain(fabs(pmt->getPhi()));
745
746 size_t n0 = min(ns[distance(module->begin(),pmt)], parameters.Nmax_PMT);
747
748 job.Fill((double) (100 + CDF[i].type), (double) n0);
749
750 while (n0 != 0) {
751
752 const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
753 const int n1 = getNumberOfPhotoElectrons(n0);
754
755 mc_hits.push_back(JHit_t(mc_hits.size() + 1,
756 pmt->getID(),
757 getHitType(CDF[i].type),
758 track->id,
759 t0 + (R * getTanThetaC() + Z) / C + t1,
760 n1));
761
762 n0 -= n1;
763 }
764 }
765
766 if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
767 job.Fill((double) (300 + CDF[i].type));
768 }
769 }
770 }
771 catch(const exception& error) {
772 job.Fill((double) (200 + CDF[i].type));
773 }
774 }
775 }
776 }
777 }
778
779 for (JTrack::const_iterator vertex = muon.begin(); vertex != muon.end(); ++vertex) {
780
781 const double Es = vertex->getEs();
782
783 if (Es >= parameters.Ecut_GeV) {
784
785 const double z0 = vertex->getZ();
786 const double t0 = vertex->getT();
787 const double DZ = geanz.getMaximum(Es);
788
789 int origin = track->id;
790
791 if (writeEMShowers) {
792 origin = event.mc_trks.size() + 1;
793 }
794
795 int number_of_hits = 0;
796
797 JDetectorSubset_t::range_type range = subdetector.getRange(z0 - maximal_distance,
798 z0 + maximal_distance);
799
800 for (JDetector::const_iterator module = range.begin(); module != range.end(); ++module) {
801
802 const double R = hypot(module->getX() - vertex->getX(),
803 module->getY() - vertex->getY());
804 const double Z = module->getZ() - z0 - DZ;
805 const double D = sqrt(R*R + Z*Z);
806 const double cd = Z / D;
807
808 for (size_t i = 0; i != CDG.size(); ++i) {
809
810 if (D < CDG[i].integral.getXmax()) {
811
812 try {
813
814 const double NPE = CDG[i].integral.getNPE(D, cd) * Es * module->size() * factor;
815 const size_t N = getPoisson(NPE);
816
817 if (N != 0) {
818
819 vector<double> npe;
820
821 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
822
823 const double R = hypot(pmt->getX() - vertex->getX(),
824 pmt->getY() - vertex->getY());
825 const double Z = pmt->getZ() - z0 - DZ;
826 const double D = sqrt(R*R + Z*Z);
827 const double cd = Z / D;
828 const double theta = pi.constrain(pmt->getTheta());
829 const double phi = pi.constrain(fabs(pmt->getPhi()));
830
831 npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor);
832 }
833
834 const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
835
836 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
837
838 const double R = hypot(pmt->getX() - vertex->getX(),
839 pmt->getY() - vertex->getY());
840 const double theta = pi.constrain(pmt->getTheta());
841 const double phi = pi.constrain(fabs(pmt->getPhi()));
842
843 size_t n0 = min(ns[distance(module->begin(),pmt)], parameters.Nmax_PMT);
844
845 job.Fill((double) (100 + CDG[i].type), (double) n0);
846
847 while (n0 != 0) {
848
849 const double dz = geanz.getLength(Es, gRandom->Rndm());
850 const double Z = pmt->getZ() - z0 - dz;
851 const double D = sqrt(R*R + Z*Z);
852 const double cd = Z / D;
853
854 const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
855 const int n1 = getNumberOfPhotoElectrons(n0);
856
857 mc_hits.push_back(JHit_t(mc_hits.size() + 1,
858 pmt->getID(),
859 getHitType(CDG[i].type),
860 origin,
861 t0 + (dz + D * getIndexOfRefraction()) / C + t1,
862 n1));
863
864 n0 -= n1;
865
866 number_of_hits += n1;
867 }
868 }
869
870 if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
871 job.Fill((double) (300 + CDG[i].type));
872 }
873 }
874 }
875 catch(const exception& error) {
876 job.Fill((double) (200 + CDG[i].type));
877 }
878 }
879 }
880 }
881
882 if (writeEMShowers && number_of_hits != 0) {
883
884 event.mc_trks.push_back(JTrk_t(origin,
885 TRACK_TYPE_PHOTON,
886 track->id,
887 track->pos + track->dir * vertex->getZ(),
888 track->dir,
889 vertex->getT(),
890 Es));
891 }
892 }
893 }
894
895 } else if (track->len > 0.0) {
896
897 // -----------------------------------------------
898 // decayed particles treated as mip (tau includes mip+deltaray)
899 // -----------------------------------------------
900
901 job.Fill(6.0);
902
903 const double z0 = 0.0;
904 const double z1 = z0 + track->len;
905 const double t0 = track->t;
906 const double E = track->E;
907
908 const JTransformation3D transformation = getTransformation(*track);
909
910 JModule buffer;
911
912 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
913
914 const JPosition3D pos = transformation.transform(module->getPosition());
915
916 const double R = pos.getX();
917 const double Z = pos.getZ() - R / getTanThetaC();
918
919 if (Z < z0 ||
920 Z > z1 ||
921 R > maximal_road_width) {
922 continue;
923 }
924
925 for (size_t i = 0; i != CDF.size(); ++i) {
926
927 double W = 1.0; // mip
928
929 if (is_deltarays(CDF[i].type)) {
930
931 if (is_tau(*track))
932 W = JDeltaRays::getEnergyLossFromTau(E); // delta-rays
933 else
934 continue;
935 }
936
937 if (R < CDF[i].integral.getXmax()) {
938
939 try {
940
941 const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
942 const size_t N = getPoisson(NPE);
943
944 if (N != 0) {
945
946 buffer = *module;
947
948 buffer.transform(transformation);
949
950 vector<double> npe;
951
952 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
953
954 const double R = pmt->getX();
955 const double theta = pi.constrain(pmt->getTheta());
956 const double phi = pi.constrain(fabs(pmt->getPhi()));
957
958 npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
959 }
960
961 const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
962
963 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
964
965 const double R = pmt->getX();
966 const double Z = pmt->getZ() - z0;
967 const double theta = pi.constrain(pmt->getTheta());
968 const double phi = pi.constrain(fabs(pmt->getPhi()));
969
970 size_t n0 = min(ns[distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
971
972 job.Fill((double) (120 + CDF[i].type), (double) n0);
973
974 while (n0 != 0) {
975
976 const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
977 const int n1 = getNumberOfPhotoElectrons(n0);
978
979 mc_hits.push_back(JHit_t(mc_hits.size() + 1,
980 pmt->getID(),
981 getHitType(CDF[i].type),
982 track->id,
983 t0 + (R * getTanThetaC() + Z) / C + t1,
984 n1));
985
986 n0 -= n1;
987 }
988 }
989
990 if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
991 job.Fill((double) (320 + CDF[i].type));
992 }
993 }
994 }
995 catch(const exception& error) {
996 job.Fill((double) (220 + CDF[i].type));
997 }
998 }
999 }
1000 }
1001
1002 if (!buffer.empty()) {
1003 job.Fill(7.0);
1004 }
1005
1006 } else if (!is_neutrino(*track)) {
1007
1008 if (JPDB::getInstance().hasPDG(track->type)) {
1009
1010 // -----------------------------------------------
1011 // electron or hadron
1012 // -----------------------------------------------
1013
1014 job.Fill(8.0);
1015
1016 double E = track->E;
1017
1018 try {
1019 E = getKineticEnergy(E, JPDB::getInstance().getPDG(track->type).mass);
1020 }
1021 catch(const exception& error) {
1022 ERROR(error.what() << endl);
1023 }
1024
1025 E = pythia(track->type, E);
1026
1027 if (E >= parameters.Ecut_GeV && cylinder.getDistance(getPosition(*track)) < parameters.Dmin_m) {
1028
1029 const double z0 = 0.0;
1030 const double t0 = track->t;
1031 const double DZ = geanz.getMaximum(E);
1032
1033 const JTransformation3D transformation = getTransformation(*track);
1034
1035 JModule buffer;
1036
1037 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
1038
1039 const JPosition3D pos = transformation.transform(module->getPosition());
1040
1041 const double R = pos.getX();
1042 const double Z = pos.getZ() - z0 - DZ;
1043 const double D = sqrt(R*R + Z*Z);
1044 const double cd = Z / D;
1045
1046 for (size_t i = 0; i != CDG.size(); ++i) {
1047
1048 if (D < CDG[i].integral.getXmax()) {
1049
1050 try {
1051
1052 const double NPE = CDG[i].integral.getNPE(D, cd) * E * module->size() * factor;
1053 const size_t N = getPoisson(NPE);
1054
1055 if (N != 0) {
1056
1057 buffer = *module;
1058
1059 buffer.transform(transformation);
1060
1061 vector<double> npe;
1062
1063 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1064
1065 const double R = pmt->getX();
1066 const double Z = pmt->getZ() - z0 - DZ;
1067 const double D = sqrt(R*R + Z*Z);
1068 const double cd = Z / D;
1069 const double theta = pi.constrain(pmt->getTheta());
1070 const double phi = pi.constrain(fabs(pmt->getPhi()));
1071
1072 npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * E * factor);
1073 }
1074
1075 const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
1076
1077 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1078
1079 const double theta = pi.constrain(pmt->getTheta());
1080 const double phi = pi.constrain(fabs(pmt->getPhi()));
1081
1082 size_t n0 = min(ns[distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
1083
1084 job.Fill((double) (140 + CDG[i].type), (double) n0);
1085
1086 while (n0 != 0) {
1087
1088 const double dz = geanz.getLength(E, gRandom->Rndm());
1089 const double Z = pmt->getZ() - z0 - dz;
1090 const double D = sqrt(R*R + Z*Z);
1091 const double cd = Z / D;
1092
1093 const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
1094 const int n1 = getNumberOfPhotoElectrons(n0);
1095
1096 mc_hits.push_back(JHit_t(mc_hits.size() + 1,
1097 pmt->getID(),
1098 getHitType(CDG[i].type, true),
1099 track->id,
1100 t0 + (dz + D * getIndexOfRefraction()) / C + t1,
1101 n1));
1102
1103 n0 -= n1;
1104 }
1105 }
1106
1107 if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
1108 job.Fill((double) (340 + CDG[i].type));
1109 }
1110 }
1111 }
1112 catch(const exception& error) {
1113 job.Fill((double) (240 + CDG[i].type));
1114 }
1115 }
1116 }
1117 }
1118
1119 if (!buffer.empty()) {
1120 job.Fill(9.0);
1121 }
1122
1123 } else {
1124 job.Fill(21.0);
1125 }
1126 }
1127 }
1128 }
1129
1130 if (!mc_hits.empty()) {
1131
1132 mc_hits.merge(parameters.Tmax_ns);
1133
1134 event.mc_hits.resize(mc_hits.size());
1135
1136 copy(mc_hits.begin(), mc_hits.end(), event.mc_hits.begin());
1137 }
1138
1139 timer.stop();
1140
1141 if (has_neutrino(event)) {
1142 cpu.Fill(log10(get_neutrino(event).E), (double) timer.usec_ucpu * 1.0e-3);
1143 }
1144
1145 if (event.mc_hits.size() >= numberOfHits) {
1146
1147 outputFile.put(event);
1148
1149 job.Fill(10.0);
1150 }
1151 }
1152 STATUS(endl);
1153
1154 outputFile.put(job);
1155 outputFile.put(cpu);
1156 outputFile.put(rms);
1157 outputFile.put(rad);
1158 outputFile.put(*gRandom);
1159
1160 JMultipleFileScanner<JTYPELIST<JMetaTypes_t, MultiHead>::typelist> io(inputFile);
1161
1162 io >> outputFile;
1163
1164 outputFile.close();
1165}
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: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: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: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
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:403
JPHYSICS::JDeltaRays::getTmin
static double getTmin()
Get minimum delta-ray kinetic energy.
Definition JDeltaRays.hh:59
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:420
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