Jpp  master_rocky-43-ge265d140c
the software that should make you happy
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 More...
 
int numberOfBins = 200
 number of bins for average CDF integral of optical module More...
 
double safetyFactor = 1.7
 safety factor for average CDF integral of optical module More...
 

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
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
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
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::JCylinder3D
Cylinder object.
Definition: JCylinder3D.hh:41
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::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
JSUPPORT::JFileRecorder
Object writing to file.
Definition: JFileRecorder.hh:44
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::getRange
const range_type & getRange() const
Get range.
Definition: JRange.hh:135
K
static const uint32_t K[64]
Definition: crypt.cc:77
GAUSS_LEGENDRE::epsilon
const double epsilon
Definition: JQuadrature.cc:21
JAANET::getAxis
JAxis3D getAxis(const Trk &track)
Get axis.
Definition: JAAnetToolkit.hh:250
JAANET::get_neutrino
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Definition: JAAnetToolkit.hh:513
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::getTime
double getTime(const Hit &hit)
Get true time of hit.
Definition: JAAnetToolkit.hh:96
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::getNPE
double getNPE(const Hit &hit)
Get true charge of hit.
Definition: JAAnetToolkit.hh:108
JAANET::TRACK_TYPE_PHOTON
@ TRACK_TYPE_PHOTON
Definition: JParticleTypes.hh:102
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
JEEP::getGITVersion
std::string getGITVersion(const std::string &tag)
Get GIT version for given GIT tag.
Definition: Jeep/JGITTags.hh:34
JLANG::getInstance
T & getInstance(const T &object)
Get static instance from temporary object.
Definition: JObject.hh:75
JLANG::make_predicate
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Definition: JPredicate.hh:128
JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20
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::H
static const double H
Planck constant [eV s].
Definition: JPhysics/JConstants.hh:29
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::gWater
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
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::SCATTERED_LIGHT_FROM_EMSHOWER
@ SCATTERED_LIGHT_FROM_EMSHOWER
scattered light from EM shower
Definition: JPDFTypes.hh:38
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
@ SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition: JPDFTypes.hh:33
JPHYSICS::SCATTERED_LIGHT_FROM_MUON
@ SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition: JPDFTypes.hh:27
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER
@ DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition: JPDFTypes.hh:37
JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS
@ DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition: JPDFTypes.hh:32
JPHYSICS::DIRECT_LIGHT_FROM_MUON
@ DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition: JPDFTypes.hh:26
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
JPHYSICS::C
static const double C
Physics constants.
Definition: JPhysics/JConstants.hh:25
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::getDate
const char * getDate()
Get current local date conform ISO-8601 standard.
Definition: JDateAndTime.hh:36
JTOOLS::accumulate
void accumulate(T &value, JBool< false > option)
Accumulation of value.
Definition: JToolsToolkit.hh:321
KM3NETDAQ::range_type
JRange< int > range_type
Definition: JDQSimulator.cc:160
create_update_issues.origin
string 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
std
Definition: JSTDTypes.hh:14
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::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::JTrack
Muon trajectory.
Definition: JSireneToolkit.hh:249
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::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.9.1