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