Jpp  17.3.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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/JDate.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/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.

sirene.png
Picture by Claudine Colnard

Note that CDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD;
The file names are obtained by replacing JPHYSICS::WILD_CARD; with

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

int main ( int  argc,
char **  argv 
)

Definition at line 195 of file JSirene.cc.

196 {
197  using namespace std;
198  using namespace JPP;
199 
200  string fileDescriptor;
201  JMultipleFileScanner<JAAnetTypes_t> inputFile;
202  JFileRecorder <JTYPELIST<JAAnetTypes_t, JMetaTypes_t, JRootTypes_t>::typelist> outputFile;
203  JLimit_t& numberOfEvents = inputFile.getLimit();
204  string detectorFile;
205  JParameters parameters;
206  bool writeEMShowers;
207  bool keep;
208  size_t numberOfHits;
209  double factor;
210  UInt_t seed;
211 
212  try {
213 
214  JProperties properties;
215 
216  properties.insert(gmake_property(parameters.Ecut_GeV));
217  properties.insert(gmake_property(parameters.Emin_GeV));
218  properties.insert(gmake_property(parameters.Dmin_m));
219  properties.insert(gmake_property(parameters.Emax_GeV));
220  properties.insert(gmake_property(parameters.Dmax_m));
221  properties.insert(gmake_property(parameters.Tmax_ns));
222  properties.insert(gmake_property(parameters.Nmax_NPE));
223  properties.insert(gmake_property(parameters.Nmax_PMT));
224 
225  properties.insert(gmake_property(numberOfBins));
226  properties.insert(gmake_property(safetyFactor));
227 
228  JParser<> zap("Main program to simulate detector response to muons and showers.");
229 
230  zap['@'] = make_field(properties) = JPARSER::initialised();
231  zap['F'] = make_field(fileDescriptor, "file name descriptor for CDF tables");
232  zap['f'] = make_field(inputFile) = JPARSER::initialised();
233  zap['o'] = make_field(outputFile) = "sirene.root";
234  zap['n'] = make_field(numberOfEvents) = JLimit::max();
235  zap['a'] = make_field(detectorFile) = "";
236  zap['s'] = make_field(writeEMShowers, "store generated EM showers in event");
237  zap['k'] = make_field(keep, "keep position of tracks");
238  zap['N'] = make_field(numberOfHits, "minimum number of hits to output event") = 1;
239  zap['U'] = make_field(factor, "scaling factor applied to light yields") = 1.0;
240  zap['S'] = make_field(seed) = 0;
241  zap['d'] = make_field(debug) = 1;
242 
243  zap(argc, argv);
244  }
245  catch(const exception &error) {
246  FATAL(error.what() << endl);
247  }
248 
249 
250  gRandom->SetSeed(seed);
251 
252 
253  const JMeta meta(argc, argv);
254 
255  const double Zbed = 0.0; // level of seabed [m]
256 
257  vector<JCDF_t> CDF;
258  vector<JCDG_t> CDG;
259 
260  CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_MUON));
261  CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_MUON));
262  CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_DELTARAYS));
263  CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_DELTARAYS));
264 
265  CDG.push_back(JCDG_t(fileDescriptor, DIRECT_LIGHT_FROM_EMSHOWER));
266  CDG.push_back(JCDG_t(fileDescriptor, SCATTERED_LIGHT_FROM_EMSHOWER));
267 
268 
269  double maximal_road_width = 0.0; // road width [m]
270  double maximal_distance = 0.0; // road width [m]
271 
272  for (size_t i = 0; i != CDF.size(); ++i) {
273 
274  DEBUG("Range CDF["<< CDF[i].type << "] " << CDF[i].function.intensity.getXmax() << " m" << endl);
275 
276  maximal_road_width = max(maximal_road_width, CDF[i].function.intensity.getXmax());
277  }
278 
279  for (size_t i = 0; i != CDG.size(); ++i) {
280 
281  DEBUG("Range CDG["<< CDG[i].type << "] " << CDG[i].function.intensity.getXmax() << " m" << endl);
282 
283  if (!is_scattered(CDF[i].type)) {
284  maximal_road_width = max(maximal_road_width, CDG[i].function.intensity.getXmax());
285  }
286 
287  maximal_distance = max(maximal_distance, CDG[i].function.intensity.getXmax());
288  }
289 
290  NOTICE("Maximal road width [m] " << maximal_road_width << endl);
291  NOTICE("Maximal distance [m] " << maximal_distance << endl);
292 
293 
294  if (detectorFile == "" || inputFile.empty()) {
295  STATUS("Nothing to be done." << endl);
296  return 0;
297  }
298 
299  JDetector detector;
300 
301  try {
302 
303  STATUS("Load detector... " << flush);
304 
305  load(detectorFile, detector);
306 
307  STATUS("OK" << endl);
308  }
309  catch(const JException& error) {
310  FATAL(error);
311  }
312 
313  // remove empty modules
314 
315  for (JDetector::iterator module = detector.begin(); module != detector.end(); ) {
316  if (!module->empty())
317  ++module;
318  else
319  module = detector.erase(module);
320  }
321 
322  vector< JSharedPointer<JRadiationInterface> > radiation;
323  vector< JSharedPointer<JACoeffInterface> > ionization;
324 
325  if (true) {
326 
327  STATUS("Setting up radiation tables... " << flush);
328 
329  //More accuracy can be achieved uncommenting the commented elements and its initializations, remember to do the same in JSeawater.hh (The code will run slower).
330 
331  const JRadiation hydrogen ( 1.0, 1.0, 40, 0.01, 0.1, 0.1);
332  const JRadiation oxygen ( 8.0, 16.0, 40, 0.01, 0.1, 0.1);
333  const JRadiation chlorine (17.0, 35.5, 40, 0.01, 0.1, 0.1);
334  const JRadiation sodium (11.0, 23.0, 40, 0.01, 0.1, 0.1);
335 #ifdef RADIATION
336  const JRadiation calcium (20.0, 40.0, 40, 0.01, 0.1, 0.1);
337  const JRadiation magnesium(12.0, 24.3, 40, 0.01, 0.1, 0.1);
338  const JRadiation potassium(19.0, 39.0, 40, 0.01, 0.1, 0.1);
339  const JRadiation sulphur (16.0, 32.0, 40, 0.01, 0.1, 0.1);
340 #endif
341 
342  JSharedPointer<JRadiation> Hydrogen (new JRadiationFunction(hydrogen, 300, 0.2, 1.0e11));
343  JSharedPointer<JRadiation> Oxygen (new JRadiationFunction(oxygen, 300, 0.2, 1.0e11));
344  JSharedPointer<JRadiation> Chlorine (new JRadiationFunction(chlorine, 300, 0.2, 1.0e11));
345  JSharedPointer<JRadiation> Sodium (new JRadiationFunction(sodium, 300, 0.2, 1.0e11));
346 #ifdef RADIATION
347  JSharedPointer<JRadiation> Calcium (new JRadiationFunction(calcium, 300, 0.2, 1.0e11));
348  JSharedPointer<JRadiation> Magnesium(new JRadiationFunction(magnesium,300, 0.2, 1.0e11));
349  JSharedPointer<JRadiation> Potassium(new JRadiationFunction(potassium,300, 0.2, 1.0e11));
350  JSharedPointer<JRadiation> Sulphur (new JRadiationFunction(sulphur, 300, 0.2, 1.0e11));
351 #endif
352 
353  radiation.push_back(new JRadiationSource(11, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), EErad_t));
354  radiation.push_back(new JRadiationSource(12, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad_t));
355  radiation.push_back(new JRadiationSource(13, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), EErad_t));
356  radiation.push_back(new JRadiationSource(14, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), EErad_t));
357 #ifdef RADIATION
358  radiation.push_back(new JRadiationSource(15, Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), EErad_t));
359  radiation.push_back(new JRadiationSource(16, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), EErad_t));
360  radiation.push_back(new JRadiationSource(17, Potassium,DENSITY_SEA_WATER * JSeaWater::K(), EErad_t));
361  radiation.push_back(new JRadiationSource(18, Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), EErad_t));
362 #endif
363 
364  radiation.push_back(new JRadiationSource(21, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), Brems_t));
365  radiation.push_back(new JRadiationSource(22, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems_t));
366  radiation.push_back(new JRadiationSource(23, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), Brems_t));
367  radiation.push_back(new JRadiationSource(24, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), Brems_t));
368 #ifdef RADIATION
369  radiation.push_back(new JRadiationSource(25, Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), Brems_t));
370  radiation.push_back(new JRadiationSource(26, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), Brems_t));
371  radiation.push_back(new JRadiationSource(27, Potassium,DENSITY_SEA_WATER * JSeaWater::K(), Brems_t));
372  radiation.push_back(new JRadiationSource(28, Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), Brems_t));
373 #endif
374 
375  radiation.push_back(new JRadiationSource(31, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), GNrad_t));
376  radiation.push_back(new JRadiationSource(32, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad_t));
377  radiation.push_back(new JRadiationSource(33, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), GNrad_t));
378  radiation.push_back(new JRadiationSource(34, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), GNrad_t));
379 #ifdef RADIATION
380  radiation.push_back(new JRadiationSource(35, Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), GNrad_t));
381  radiation.push_back(new JRadiationSource(36, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), GNrad_t));
382  radiation.push_back(new JRadiationSource(37, Potassium,DENSITY_SEA_WATER * JSeaWater::K(), GNrad_t));
383  radiation.push_back(new JRadiationSource(38, Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), GNrad_t));
384 #endif
385 
386  radiation.push_back(new JDISSource(100, DENSITY_SEA_WATER));
387 
388  ionization.push_back(new JACoeffSource(Oxygen, DENSITY_SEA_WATER * JSeaWater::O()));
389  ionization.push_back(new JACoeffSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl()));
390  ionization.push_back(new JACoeffSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H()));
391  ionization.push_back(new JACoeffSource(Sodium, DENSITY_SEA_WATER * JSeaWater::Na()));
392 #ifdef RADIATION
393  ionization.push_back(new JACoeffSource(Calcium, DENSITY_SEA_WATER * JSeaWater::Ca()));
394  ionization.push_back(new JACoeffSource(Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg()));
395  ionization.push_back(new JACoeffSource(Potassium,DENSITY_SEA_WATER * JSeaWater::K()));
396  ionization.push_back(new JACoeffSource(Sulphur, DENSITY_SEA_WATER * JSeaWater::S()));
397 #endif
398 
399  STATUS("OK" << endl);
400  }
401 
402 
403  Vec center(0,0,0);
404  Head header;
405 
406  try {
407 
408  header = inputFile.getHeader();
409 
410  JHead buffer(header);
411 
412  center = get<Vec>(buffer);
413 
414  buffer.simul.push_back(JAANET::simul());
415 
416  buffer.simul.rbegin()->program = APPLICATION_JSIRENE;
417  buffer.simul.rbegin()->version = getGITVersion();
418  buffer.simul.rbegin()->date = getDate();
419  buffer.simul.rbegin()->time = getTime();
420 
421  buffer.detector.push_back(JAANET::detector());
422 
423  buffer.detector.rbegin()->program = APPLICATION_JSIRENE;
424  buffer.detector.rbegin()->filename = detectorFile;
425 
426  buffer.push(&JHead::simul);
427  buffer.push(&JHead::detector);
428 
429  if (!keep) {
430 
431  buffer.coord_origin = coord_origin(0.0, 0.0, 0.0);
432  buffer.can.zmin += center.z;
433  buffer.can.zmax += center.z;
434 
435  buffer.push(&JHead::coord_origin);
436  buffer.push(&JHead::can);
437  }
438 
439  const JCircle2D circle(detector.begin(), detector.end());
440 
441  center += Vec(circle.getX(), circle.getY(), 0.0);
442 
443  copy(buffer, header);
444  }
445  catch(const JException& error) {
446  FATAL(error);
447  }
448 
449  if (!keep)
450  NOTICE("Offset applied to true tracks is: " << center << endl);
451  else
452  NOTICE("No offset applied to true tracks." << endl);
453 
454  JCylinder3D cylinder(detector.begin(), detector.end());
455 
456  cylinder.addMargin(maximal_distance);
457 
458  if (cylinder.getZmin() < Zbed) {
459  cylinder.setZmin(Zbed);
460  }
461 
462  NOTICE("Light generation volume: " << cylinder << endl);
463 
464  TH1D job("job", NULL, 400, 0.5, 400.5);
465  TProfile cpu("cpu", NULL, 16, 0.0, 8.0);
466  TProfile2D rms("rms", NULL, 16, 0.0, 8.0, 201, -0.5, 200.5);
467  TProfile2D rad("rad", NULL, 16, 0.0, 8.0, 201, -0.5, 200.5);
468 
469 
470  outputFile.open();
471 
472  if (!outputFile.is_open()) {
473  FATAL("Error opening file " << outputFile << endl);
474  }
475 
476  outputFile.put(meta);
477  outputFile.put(header);
478  outputFile.put(*gRandom);
479 
480  const double epsilon = 1.0e-6; // precision angle [rad]
481  const JRange<double> pi(epsilon, PI - epsilon); // constrain angle
482 
483  JTimer timer;
484 
485  for (JMultipleFileScanner<Evt>& in = inputFile; in.hasNext(); ) {
486 
487  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
488 
489  job.Fill(1.0);
490 
491  Evt* evt = in.next();
492 
493  if (!keep) {
494  for (vector<Trk>::iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
495  track->pos += center;
496  }
497  }
498 
499  Evt event(*evt); // output
500 
501  event.mc_hits.clear();
502 
503  JHits_t mc_hits; // temporary buffer
504 
505  timer.reset();
506  timer.start();
507 
508  for (vector<Trk>::const_iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
509 
510  if (!track->is_finalstate()) {
511  continue; // only final state particles produce light
512  }
513 
514  if (is_muon(*track)) {
515 
516  // -----------------------------------------------
517  // muon
518  // -----------------------------------------------
519 
520  job.Fill(2.0);
521 
522  typedef JDetectorSubset<JModuleComparator3Z_t> JDetectorSubset_t;
523 
524  const JCylinder3D::intersection_type intersection = cylinder.getIntersection(getAxis(*track));
525 
526  double Zmin = intersection.first;
527  double Zmax = intersection.second;
528 
529  if (Zmax - Zmin <= parameters.Dmin_m) {
530  continue;
531  }
532 
533  JVertex vertex(0.0, track->t, track->E); // start of muon
534 
535  if (track->pos.z < Zbed) { // propagate muon through rock
536 
537  if (track->dir.z > 0.0)
538  vertex.step(gRock, (Zbed - track->pos.z) / track->dir.z);
539  else
540  continue;
541  }
542 
543  if (vertex.getZ() < Zmin) { // propagate muon through water
544  vertex.step(gWater, Zmin - vertex.getZ());
545  }
546 
547  if (vertex.getRange() <= parameters.Dmin_m) {
548  continue;
549  }
550 
551  job.Fill(3.0);
552 
553  const JDetectorSubset_t subdetector(detector, getAxis(*track), maximal_road_width);
554 
555  if (subdetector.empty()) {
556  continue;
557  }
558 
559  job.Fill(4.0);
560 
561  JTrack muon(vertex); // propagate muon trough detector
562 
563  while (vertex.getE() >= parameters.Emin_GeV && vertex.getZ() < Zmax) {
564 
565  const int N = radiation.size();
566 
567  double li[N]; // inverse interaction lengths
568  double ls = 1.0e-5; // minimal total inverse interaction length [m^-1]
569 
570  for (int i = 0; i != N; ++i) {
571  ls += li[i] = radiation[i]->getInverseInteractionLength(vertex.getE());
572  }
573 
574  double As = 0.0; // ionization energy loss
575 
576  for (size_t i = 0; i != ionization.size(); ++i) {
577  As += ionization[i]->getA(vertex.getE());
578  }
579 
580  double step = gRandom->Exp(1.0) / ls; // distance to next radiation process
581  double range = vertex.getRange(As); // range of muon
582 
583  if (vertex.getE() < parameters.Emax_GeV) { // limited step size
584  if (parameters.Dmax_m < range) {
585  range = parameters.Dmax_m;
586  }
587  }
588 
589  double ts = getThetaMCS(vertex.getE(), min(step,range)); // multiple Coulomb scattering angle [rad]
590  double T2 = ts*ts; //
591 
592  rms.Fill(log10(vertex.getE()), (Double_t) 0, ts*ts);
593 
594  vertex.getDirection() += getRandomDirection(T2/3.0); // multiple Coulomb planar scattering
595 
596  vertex.step(As, min(step,range)); // ionization energy loss
597 
598  double Es = 0.0; // shower energy [GeV]
599 
600  if (step < range) {
601 
602  if (vertex.getE() >= parameters.Emin_GeV) {
603 
604  double y = gRandom->Uniform(ls);
605 
606  for (int i = 0; i != N; ++i) {
607 
608  y -= li[i];
609 
610  if (y < 0.0) {
611 
612  Es = radiation[i]->getEnergyOfShower(vertex.getE()); // shower energy [GeV]
613  ts = radiation[i]->getThetaRMS(vertex.getE(), Es); // scattering angle [rad]
614 
615  T2 += ts*ts;
616 
617  rms.Fill(log10(vertex.getE()), (Double_t) radiation[i]->getID(), ts*ts);
618  rad.Fill(log10(vertex.getE()), (Double_t) radiation[i]->getID(), Es);
619 
620  break;
621  }
622  }
623  }
624  }
625 
626  vertex.applyEloss(getRandomDirection(T2), Es);
627 
628  muon.push_back(vertex);
629  }
630 
631  // add muon end point
632 
633  if (vertex.getZ() < Zmax && vertex.getRange() > 0.0) {
634 
635  vertex.step(vertex.getRange());
636 
637  muon.push_back(vertex);
638  }
639 
640  // add information to output muon
641 
642  vector<Trk>::iterator trk = find_if(event.mc_trks.begin(),
643  event.mc_trks.end(),
644  make_predicate(&Trk::id, track->id));
645 
646  if (trk != event.mc_trks.end()) {
647  trk->len = (muon.rbegin()->getZ() < Zmax ? +1 : -1) * (muon.rbegin()->getZ() - muon.begin()->getZ());
648  trk->setusr(mc_usr_keys::energy_lost_in_can, muon.begin()->getE() - muon.rbegin()->getE());
649  }
650 
651  for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
652 
653  const double z0 = muon.begin()->getZ();
654  const double t0 = muon.begin()->getT();
655  const double Z = module->getZ() - module->getX() / getTanThetaC();
656 
657  if (Z >= muon.begin()->getZ() && Z <= muon.rbegin()->getZ()) {
658 
659  const JVector2D pos = muon.getPosition(Z);
660  const double R = hypot(module->getX() - pos.getX(),
661  module->getY() - pos.getY());
662 
663  for (size_t i = 0; i != CDF.size(); ++i) {
664 
665  if (R < CDF[i].integral.getXmax()) {
666 
667  try {
668 
669  double W = 1.0; // mip
670 
671  if (is_deltarays(CDF[i].type)) {
672  W = getDeltaRaysFromMuon(muon.getE(Z)); // delta-rays
673  }
674 
675  const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
676  const int N = gRandom->Poisson(NPE);
677 
678  if (N != 0) {
679 
680  vector<double> npe;
681 
682  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
683 
684  const double R = hypot(pmt->getX() - pos.getX(),
685  pmt->getY() - pos.getY());
686  const double theta = pi.constrain(pmt->getTheta());
687  const double phi = pi.constrain(fabs(pmt->getPhi()));
688 
689  npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
690  }
691 
692  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
693 
694  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
695 
696  const double R = hypot(pmt->getX() - pos.getX(),
697  pmt->getY() - pos.getY());
698  const double Z = pmt->getZ() - z0;
699  const double theta = pi.constrain(pmt->getTheta());
700  const double phi = pi.constrain(fabs(pmt->getPhi()));
701 
702  size_t n0 = min(ns[distance(module->begin(),pmt)], parameters.Nmax_PMT);
703 
704  job.Fill((double) (100 + CDF[i].type), (double) n0);
705 
706  while (n0 != 0) {
707 
708  const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
709  const int n1 = getNumberOfPhotoElectrons(n0);
710 
711  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
712  pmt->getID(),
713  getHitType(CDF[i].type),
714  track->id,
715  t0 + (R * getTanThetaC() + Z) / C + t1,
716  n1));
717 
718  n0 -= n1;
719  }
720  }
721 
722  if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
723  job.Fill((double) (300 + CDF[i].type));
724  }
725  }
726  }
727  catch(const exception& error) {
728  job.Fill((double) (200 + CDF[i].type));
729  }
730  }
731  }
732  }
733  }
734 
735  for (JTrack::const_iterator vertex = muon.begin(); vertex != muon.end(); ++vertex) {
736 
737  const double Es = vertex->getEs();
738 
739  if (Es >= parameters.Ecut_GeV) {
740 
741  const double z0 = vertex->getZ();
742  const double t0 = vertex->getT();
743  const double DZ = geanz.getMaximum(Es);
744 
745  int origin = track->id;
746 
747  if (writeEMShowers) {
748  origin = event.mc_trks.size() + 1;
749  }
750 
751  int number_of_hits = 0;
752 
753  JDetectorSubset_t::range_type range = subdetector.getRange(z0 - maximal_distance,
754  z0 + maximal_distance);
755 
756  for (JDetector::const_iterator module = range.begin(); module != range.end(); ++module) {
757 
758  const double R = hypot(module->getX() - vertex->getX(),
759  module->getY() - vertex->getY());
760  const double Z = module->getZ() - z0 - DZ;
761  const double D = sqrt(R*R + Z*Z);
762  const double cd = Z / D;
763 
764  for (size_t i = 0; i != CDG.size(); ++i) {
765 
766  if (D < CDG[i].integral.getXmax()) {
767 
768  try {
769 
770  const double NPE = CDG[i].integral.getNPE(D, cd) * Es * module->size() * factor;
771  const int N = gRandom->Poisson(NPE);
772 
773  if (N != 0) {
774 
775  vector<double> npe;
776 
777  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
778 
779  const double R = hypot(pmt->getX() - vertex->getX(),
780  pmt->getY() - vertex->getY());
781  const double Z = pmt->getZ() - z0 - DZ;
782  const double D = sqrt(R*R + Z*Z);
783  const double cd = Z / D;
784  const double theta = pi.constrain(pmt->getTheta());
785  const double phi = pi.constrain(fabs(pmt->getPhi()));
786 
787  npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor);
788  }
789 
790  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
791 
792  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
793 
794  const double R = hypot(pmt->getX() - vertex->getX(),
795  pmt->getY() - vertex->getY());
796  const double theta = pi.constrain(pmt->getTheta());
797  const double phi = pi.constrain(fabs(pmt->getPhi()));
798 
799  size_t n0 = min(ns[distance(module->begin(),pmt)], parameters.Nmax_PMT);
800 
801  job.Fill((double) (100 + CDG[i].type), (double) n0);
802 
803  while (n0 != 0) {
804 
805  const double dz = geanz.getLength(Es, gRandom->Rndm());
806  const double Z = pmt->getZ() - z0 - dz;
807  const double D = sqrt(R*R + Z*Z);
808  const double cd = Z / D;
809 
810  const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
811  const int n1 = getNumberOfPhotoElectrons(n0);
812 
813  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
814  pmt->getID(),
815  getHitType(CDG[i].type),
816  origin,
817  t0 + (dz + D * getIndexOfRefraction()) / C + t1,
818  n1));
819 
820  n0 -= n1;
821 
822  number_of_hits += n1;
823  }
824  }
825 
826  if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
827  job.Fill((double) (300 + CDG[i].type));
828  }
829  }
830  }
831  catch(const exception& error) {
832  job.Fill((double) (200 + CDG[i].type));
833  }
834  }
835  }
836  }
837 
838  if (writeEMShowers && number_of_hits != 0) {
839 
840  event.mc_trks.push_back(JTrk_t(origin,
841  TRACK_TYPE_PHOTON,
842  track->id,
843  track->pos + track->dir * vertex->getZ(),
844  track->dir,
845  vertex->getT(),
846  Es));
847  }
848  }
849  }
850 
851  } else if (track->len > 0.0) {
852 
853  // -----------------------------------------------
854  // decayed particles treated as mip (tau includes mip+deltaray)
855  // -----------------------------------------------
856 
857  job.Fill(6.0);
858 
859  const double z0 = 0.0;
860  const double z1 = z0 + track->len;
861  const double t0 = track->t;
862  const double E = track->E;
863 
864  const JTransformation3D transformation = getTransformation(*track);
865 
866  JModule buffer;
867 
868  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
869 
870  const JPosition3D pos = transformation.transform(module->getPosition());
871 
872  const double R = pos.getX();
873  const double Z = pos.getZ() - R / getTanThetaC();
874 
875  if (Z < z0 ||
876  Z > z1 ||
877  R > maximal_road_width) {
878  continue;
879  }
880 
881  for (size_t i = 0; i != CDF.size(); ++i) {
882 
883  double W = 1.0; // mip
884 
885  if (is_deltarays(CDF[i].type)) {
886 
887  if (is_tau(*track))
888  W = getDeltaRaysFromTau(E); // delta-rays
889  else
890  continue;
891  }
892 
893  if (R < CDF[i].integral.getXmax()) {
894 
895  try {
896 
897  const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
898  const int N = gRandom->Poisson(NPE);
899 
900  if (N != 0) {
901 
902  buffer = *module;
903 
904  buffer.transform(transformation);
905 
906  vector<double> npe;
907 
908  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
909 
910  const double R = pmt->getX();
911  const double theta = pi.constrain(pmt->getTheta());
912  const double phi = pi.constrain(fabs(pmt->getPhi()));
913 
914  npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
915  }
916 
917  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
918 
919  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
920 
921  const double R = pmt->getX();
922  const double Z = pmt->getZ() - z0;
923  const double theta = pi.constrain(pmt->getTheta());
924  const double phi = pi.constrain(fabs(pmt->getPhi()));
925 
926  size_t n0 = min(ns[distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
927 
928  job.Fill((double) (120 + CDF[i].type), (double) n0);
929 
930  while (n0 != 0) {
931 
932  const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
933  const int n1 = getNumberOfPhotoElectrons(n0);
934 
935  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
936  pmt->getID(),
937  getHitType(CDF[i].type),
938  track->id,
939  t0 + (R * getTanThetaC() + Z) / C + t1,
940  n1));
941 
942  n0 -= n1;
943  }
944  }
945 
946  if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
947  job.Fill((double) (320 + CDF[i].type));
948  }
949  }
950  }
951  catch(const exception& error) {
952  job.Fill((double) (220 + CDF[i].type));
953  }
954 
955  }
956  }
957  }
958 
959  if (!buffer.empty()) {
960  job.Fill(7.0);
961  }
962 
963  } else if (!is_neutrino(*track)) {
964 
965  if (JPDB::getInstance().hasPDG(track->type)) {
966 
967  // -----------------------------------------------
968  // electron or hadron
969  // -----------------------------------------------
970 
971  job.Fill(8.0);
972 
973  double E = track->E;
974 
975  try {
976  E = getKineticEnergy(E, JPDB::getInstance().getPDG(track->type).mass);
977  }
978  catch(const exception& error) {
979  ERROR(error.what() << endl);
980  }
981 
982  E = pythia(track->type, E);
983 
984  if (E >= parameters.Ecut_GeV && cylinder.getDistance(getPosition(*track)) < parameters.Dmin_m) {
985 
986  const double z0 = 0.0;
987  const double t0 = track->t;
988  const double DZ = geanz.getMaximum(E);
989 
990  const JTransformation3D transformation = getTransformation(*track);
991 
992  JModule buffer;
993 
994  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
995 
996  const JPosition3D pos = transformation.transform(module->getPosition());
997 
998  const double R = pos.getX();
999  const double Z = pos.getZ() - z0 - DZ;
1000  const double D = sqrt(R*R + Z*Z);
1001  const double cd = Z / D;
1002 
1003  for (size_t i = 0; i != CDG.size(); ++i) {
1004 
1005  if (D < CDG[i].integral.getXmax()) {
1006 
1007  try {
1008 
1009  const double NPE = CDG[i].integral.getNPE(D, cd) * E * module->size() * factor;
1010  const int N = gRandom->Poisson(NPE);
1011 
1012  if (N != 0) {
1013 
1014  buffer = *module;
1015 
1016  buffer.transform(transformation);
1017 
1018  vector<double> npe;
1019 
1020  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1021 
1022  const double R = pmt->getX();
1023  const double Z = pmt->getZ() - z0 - DZ;
1024  const double D = sqrt(R*R + Z*Z);
1025  const double cd = Z / D;
1026  const double theta = pi.constrain(pmt->getTheta());
1027  const double phi = pi.constrain(fabs(pmt->getPhi()));
1028 
1029  npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * E * factor);
1030  }
1031 
1032  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
1033 
1034  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1035 
1036  const double theta = pi.constrain(pmt->getTheta());
1037  const double phi = pi.constrain(fabs(pmt->getPhi()));
1038 
1039  size_t n0 = min(ns[distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
1040 
1041  job.Fill((double) (140 + CDG[i].type), (double) n0);
1042 
1043  while (n0 != 0) {
1044 
1045  const double dz = geanz.getLength(E, gRandom->Rndm());
1046  const double Z = pmt->getZ() - z0 - dz;
1047  const double D = sqrt(R*R + Z*Z);
1048  const double cd = Z / D;
1049 
1050  const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
1051  const int n1 = getNumberOfPhotoElectrons(n0);
1052 
1053  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
1054  pmt->getID(),
1055  getHitType(CDG[i].type, true),
1056  track->id,
1057  t0 + (dz + D * getIndexOfRefraction()) / C + t1,
1058  n1));
1059 
1060  n0 -= n1;
1061  }
1062  }
1063 
1064  if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
1065  job.Fill((double) (340 + CDG[i].type));
1066  }
1067  }
1068  }
1069  catch(const exception& error) {
1070  job.Fill((double) (240 + CDG[i].type));
1071  }
1072  }
1073  }
1074  }
1075 
1076  if (!buffer.empty()) {
1077  job.Fill(9.0);
1078  }
1079 
1080  } else {
1081  job.Fill(21.0);
1082  }
1083  }
1084  }
1085  }
1086 
1087  if (!mc_hits.empty()) {
1088 
1089  mc_hits.merge(parameters.Tmax_ns);
1090 
1091  event.mc_hits.resize(mc_hits.size());
1092 
1093  copy(mc_hits.begin(), mc_hits.end(), event.mc_hits.begin());
1094  }
1095 
1096  timer.stop();
1097 
1098  if (has_neutrino(event)) {
1099  cpu.Fill(log10(get_neutrino(event).E), (double) timer.usec_ucpu * 1.0e-3);
1100  }
1101 
1102  if (event.mc_hits.size() >= numberOfHits) {
1103 
1104  outputFile.put(event);
1105 
1106  job.Fill(10.0);
1107  }
1108  }
1109  STATUS(endl);
1110 
1111  outputFile.put(job);
1112  outputFile.put(cpu);
1113  outputFile.put(rms);
1114  outputFile.put(rad);
1115  outputFile.put(*gRandom);
1116 
1117  JMultipleFileScanner<JTYPELIST<JMetaTypes_t, MultiHead>::typelist> io(inputFile);
1118 
1119  io >> outputFile;
1120 
1121  outputFile.close();
1122 }
mc_usr_keys::energy_lost_in_can
const char *const energy_lost_in_can
Definition: io_ascii.hh:46
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
JGEOMETRY2D::JVector2D
Data structure for vector in two dimensions.
Definition: JVector2D.hh:32
JSUPPORT::JFileRecorder
Object writing to file.
Definition: JFileRecorder.hh:41
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1517
JLANG::JException
General exception.
Definition: JException.hh:23
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
K
static const uint32_t K[64]
Definition: crypt.cc:77
JPHYSICS::EErad_t
static const JRadiationSource_t EErad_t
Definition: JRadiation.hh:510
JDETECTOR::JModule
Data structure for a composite optical module.
Definition: JModule.hh:68
E
then usage E
Definition: JMuonPostfit.sh:35
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
DZ
then let DZ
Definition: module-Z:fit.sh:71
JSIRENE::JTrk_t
Auxiliary class to set-up Trk.
Definition: JSirene.hh:188
JAANET::getTransformation
JTransformation3D getTransformation(const Trk &track)
Get transformation.
Definition: JAAnetToolkit.hh:283
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
JGEOMETRY2D::JCircle2D
Data structure for circle in two dimensions.
Definition: JCircle2D.hh:33
gmake_property
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
Definition: JProperties.hh:1291
safetyFactor
double safetyFactor
safety factor for average CDF integral of optical module
Definition: JSirene.cc:70
JAANET::simul
Generator for simulation.
Definition: JHead.hh:513
JPHYSICS::getIndexOfRefraction
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
Definition: JPhysics/JConstants.hh:134
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER
scattered light from EM shower
Definition: JPDFTypes.hh:38
JAANET::has_neutrino
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Definition: JAAnetToolkit.hh:477
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:89
JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition: JAAnetToolkit.hh:339
JPHYSICS::JACoeffSource
Implementation for calculation of ionization constant.
Definition: JACoeffSource.hh:48
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
JCALIBRATE::JParameters
Auxiliary class for PMT parameters including threshold.
Definition: JParameters.hh:21
JPHYSICS::getDeltaRaysFromTau
double getDeltaRaysFromTau(const double E)
Equivalent EM-shower energy due to delta-rays per unit tau track length.
Definition: JPDFToolkit.hh:95
JEEP::JProperties
Utility class to parse parameter values.
Definition: JProperties.hh:496
N
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
Definition: JShowerPostfit.sh:95
JPHYSICS::H
static const double H
Planck constant [eV s].
Definition: JPhysics/JConstants.hh:29
parameters
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
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::gWater
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
JPHYSICS::DENSITY_SEA_WATER
static const double DENSITY_SEA_WATER
Fixed environment values.
Definition: JPhysics/JConstants.hh:44
JSIRENE::pythia
static const JPythia pythia
Function object for relative light yield as a function of GEANT particle code.
Definition: JPythia.hh:96
JAANET::getTime
double getTime(const Hit &hit)
Get true time of hit.
Definition: JAAnetToolkit.hh:93
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JGEOMETRY2D::JVector2D::getY
double getY() const
Get y position.
Definition: JVector2D.hh:74
APPLICATION_JSIRENE
static const char *const APPLICATION_JSIRENE
detector simulation
Definition: applications.hh:19
JPHYSICS::DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition: JPDFTypes.hh:26
JAANET::coord_origin
Coordinate origin.
Definition: JHead.hh:717
JPHYSICS::JRadiationFunction
Fast implementation of class JRadiation.
Definition: JRadiationFunction.hh:32
JPHYSICS::JRadiationSource
Implementation for calculation of inverse interaction length and shower energy.
Definition: JRadiationSource.hh:89
JSIRENE::JHits_t
Auxiliary data structure for list of hits with hit merging capability.
Definition: JSirene.hh:137
JAANET::is_neutrino
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
Definition: JAAnetToolkit.hh:321
JPHYSICS::C
static const double C
Physics constants.
Definition: JPhysics/JConstants.hh:25
JSIRENE::getNumberOfPhotoElectrons
struct JSIRENE::@63 getNumberOfPhotoElectrons
Auxiliary data structure for determination of number of photo-electrons.
JPHYSICS::SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition: JPDFTypes.hh:27
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
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::JDISSource
Implementation for calculation of inverse interaction length and shower energy due to deep-inelastic ...
Definition: JRadiationSource.hh:160
Vec
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition: Vec.hh:12
JGEOMETRY3D::JCylinder3D
Cylinder object.
Definition: JCylinder3D.hh:39
JAANET::detector
Detector file.
Definition: JHead.hh:226
JEEP::getGITVersion
std::string getGITVersion(const std::string &tag)
Get GIT version for given GIT tag.
Definition: Jeep/JGITTags.hh:34
JAANET::getKineticEnergy
double getKineticEnergy(const double E, const double m)
Get kinetic energy of particle with given mass.
Definition: JAAnetToolkit.hh:713
JPHYSICS::JRadiation
Auxiliary class for the calculation of the muon radiative cross sections.
Definition: JRadiation.hh:34
JLANG::getInstance
T & getInstance(const T &object)
Get static instance from temporary object.
Definition: JObject.hh:75
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1993
log10
set_variable E_E log10(E_{fit}/E_{#mu})"
JPHYSICS::geanz
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
JGEOMETRY2D::JVector2D::getX
double getX() const
Get x position.
Definition: JVector2D.hh:63
JAANET::getAxis
JAxis3D getAxis(const Trk &track)
Get axis.
Definition: JAAnetToolkit.hh:247
JEEP::JTimer
Auxiliary class for CPU timing and usage.
Definition: JTimer.hh:32
JAANET::getPosition
JPosition3D getPosition(const Vec &pos)
Get position.
Definition: JAAnetToolkit.hh:175
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition: JPDFTypes.hh:33
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition: JPDFTypes.hh:37
JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20
JSIRENE::JTrack
Muon trajectory.
Definition: JSireneToolkit.hh:247
JGEOMETRY3D::JCylinder3D::addMargin
void addMargin(const double D)
Add (safety) margin.
Definition: JCylinder3D.hh:173
JPHYSICS::Brems_t
static const JRadiationSource_t Brems_t
Definition: JRadiation.hh:511
JDETECTOR::JDetectorSubset
Detector subset with binary search functionality.
Definition: JDetectorSubset.hh:119
JDETECTOR::JDetectorSubset_t
Detector subset without binary search functionality.
Definition: JDetectorSubset.hh:38
JAANET::JHead
Monte Carlo run header.
Definition: JHead.hh:1167
JSIRENE::JVertex
Vertex of energy loss of muon.
Definition: JSireneToolkit.hh:334
Head
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:65
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
JSYSTEM::getDate
const char * getDate()
Get current date conform ISO-8601 standard.
Definition: JSystem/JDate.hh:33
Trk::id
int id
track identifier
Definition: Trk.hh:16
JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition: JPDFTypes.hh:32
range
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
Definition: module-Z:fit.sh:118
R
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:486
JPHYSICS::JGeanz::getMaximum
double getMaximum(const double E) const
Get depth of shower maximum.
Definition: JGeanz.hh:162
JSUPPORT::JMultipleFileScanner
General purpose class for object reading from a list of file names.
Definition: JMultipleFileScanner.hh:199
JPHYSICS::is_scattered
bool is_scattered(const int pdf)
Test if given PDF type corresponds to scattered light.
Definition: JPDFTypes.hh:165
JDETECTOR::JModule::transform
void transform(const JRotation3D &R, const JVector3D &pos)
Transformation of geometry (see method JGEOMETRY3D::JPosition3D::transform(const JRotation3D&amp;, const JVector3D&amp;)).
Definition: JModule.hh:346
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
JSIRENE::JHit_t
Auxiliary class to set-up Hit.
Definition: JSirene.hh:57
JGEOMETRY3D::JVector3D::getX
double getX() const
Get x position.
Definition: JVector3D.hh:94
JPHYSICS::getDeltaRaysFromMuon
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:67
JAANET::copy
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:162
JSYSTEM::ls
Auxiliary data structure to list files in directory.
Definition: JSystemToolkit.hh:302
JPHYSICS::getTanThetaC
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Definition: JPhysics/JConstants.hh:156
JAANET::getNPE
double getNPE(const Hit &hit)
Get true charge of hit.
Definition: JAAnetToolkit.hh:105
Z
do set_variable MODULE getModule a $WORKDIR detector_a datx L $STRING JEditDetector a $WORKDIR detector_a datx M $MODULE setz o $WORKDIR detector_a datx JEditDetector a $WORKDIR detector_b datx M $MODULE setz o $WORKDIR detector_b datx done echo Output stored at $WORKDIR detector_a datx and $WORKDIR tripod_a txt JDrawDetector2D a $WORKDIR detector_a datx a $WORKDIR detector_b datx L BL o detector $FORMAT $BATCH JDrawDetector2D T $WORKDIR tripod_a txt T $WORKDIR tripod_b txt L BL o tripod $FORMAT $BATCH JCompareDetector a $WORKDIR detector_a datx b $WORKDIR detector_b datx o $WORKDIR abc root &dev null for KEY in X Y Z
Definition: JAlignDetector.sh:81
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
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
JAANET::get_neutrino
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Definition: JAAnetToolkit.hh:497
detector
do set_variable DETECTOR_TXT $WORKDIR detector
Definition: detector-XY:fit1d.sh:38
JAANET::is_tau
bool is_tau(const Trk &track)
Test whether given track is a (anti-)tau.
Definition: JAAnetToolkit.hh:347
JPHYSICS::GNrad_t
static const JRadiationSource_t GNrad_t
Definition: JRadiation.hh:512
in
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:46
JGEOMETRY3D::JTransformation3D::transform
JPosition3D transform(const JPosition3D &pos) const
Transform position.
Definition: JTransformation3D.hh:49
D
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
numberOfBins
int numberOfBins
number of bins for average CDF integral of optical module
Definition: JSirene.cc:69
JGEOMETRY3D::JVector3D::getZ
double getZ() const
Get z position.
Definition: JVector3D.hh:115
Evt::mc_trks
std::vector< Trk > mc_trks
MC: list of MC truth tracks.
Definition: Evt.hh:49
GAUSS_LEGENDRE::epsilon
const double epsilon
Definition: JQuadrature.cc:21
debug
int debug
debug level
Definition: archive-put-wiki-detectors.sh:92
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20
JTOOLS::accumulate
void accumulate(T &value, JBool< false > option)
Accumulation of value.
Definition: JToolsToolkit.hh:321
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62

Variable Documentation

int debug

debug level

Definition at line 68 of file JSirene.cc.

int numberOfBins = 200

number of bins for average CDF integral of optical module

Definition at line 69 of file JSirene.cc.

double safetyFactor = 1.7

safety factor for average CDF integral of optical module

Definition at line 70 of file JSirene.cc.

Generated by   doxygen 1.8.5