Jpp  18.3.0-209-g56ce19a
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/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/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::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

int main ( int  argc,
char **  argv 
)

Definition at line 197 of file JSirene.cc.

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