Jpp  18.0.0-rc.4
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::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  if (buffer.is_valid(&JHead::fixedcan)) {
432 
433  buffer.fixedcan.xcenter += center.x;
434  buffer.fixedcan.ycenter += center.y;
435  buffer.fixedcan.zmin += center.z;
436  buffer.fixedcan.zmax += center.z;
437 
438  buffer.push(&JHead::fixedcan);
439 
440  } else if (buffer.is_valid(&JHead::can)) {
441 
442  buffer.can.zmin += center.z;
443  buffer.can.zmax += center.z;
444 
445  buffer.push(&JHead::can);
446  }
447 
448  buffer.coord_origin = coord_origin(0.0, 0.0, 0.0);
449 
450  buffer.push(&JHead::coord_origin);
451  }
452 
453  const JCircle2D circle(detector.begin(), detector.end());
454 
455  center += Vec(circle.getX(), circle.getY(), 0.0);
456 
457  copy(buffer, header);
458  }
459  catch(const JException& error) {
460  FATAL(error);
461  }
462 
463  if (!keep)
464  NOTICE("Offset applied to true tracks is: " << center << endl);
465  else
466  NOTICE("No offset applied to true tracks." << endl);
467 
468  JCylinder3D cylinder(detector.begin(), detector.end());
469 
470  cylinder.addMargin(maximal_distance);
471 
472  if (cylinder.getZmin() < Zbed) {
473  cylinder.setZmin(Zbed);
474  }
475 
476  NOTICE("Light generation volume: " << cylinder << endl);
477 
478  TH1D job("job", NULL, 400, 0.5, 400.5);
479  TProfile cpu("cpu", NULL, 16, 0.0, 8.0);
480  TProfile2D rms("rms", NULL, 16, 0.0, 8.0, 201, -0.5, 200.5);
481  TProfile2D rad("rad", NULL, 16, 0.0, 8.0, 201, -0.5, 200.5);
482 
483 
484  outputFile.open();
485 
486  if (!outputFile.is_open()) {
487  FATAL("Error opening file " << outputFile << endl);
488  }
489 
490  outputFile.put(meta);
491  outputFile.put(header);
492  outputFile.put(*gRandom);
493 
494  const double epsilon = 1.0e-6; // precision angle [rad]
495  const JRange<double> pi(epsilon, PI - epsilon); // constrain angle
496 
497  JTimer timer;
498 
499  for (JMultipleFileScanner<Evt>& in = inputFile; in.hasNext(); ) {
500 
501  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
502 
503  job.Fill(1.0);
504 
505  Evt* evt = in.next();
506 
507  if (!keep) {
508  for (vector<Trk>::iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
509  track->pos += center;
510  }
511  }
512 
513  Evt event(*evt); // output
514 
515  event.mc_hits.clear();
516 
517  JHits_t mc_hits; // temporary buffer
518 
519  timer.reset();
520  timer.start();
521 
522  for (vector<Trk>::const_iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
523 
524  if (!track->is_finalstate()) {
525  continue; // only final state particles produce light
526  }
527 
528  if (is_muon(*track)) {
529 
530  // -----------------------------------------------
531  // muon
532  // -----------------------------------------------
533 
534  job.Fill(2.0);
535 
536  typedef JDetectorSubset<JModuleComparator3Z_t> JDetectorSubset_t;
537 
538  const JCylinder3D::intersection_type intersection = cylinder.getIntersection(getAxis(*track));
539 
540  double Zmin = intersection.first;
541  double Zmax = intersection.second;
542 
543  if (Zmax - Zmin <= parameters.Dmin_m) {
544  continue;
545  }
546 
547  JVertex vertex(0.0, track->t, track->E); // start of muon
548 
549  if (track->pos.z < Zbed) { // propagate muon through rock
550 
551  if (track->dir.z > 0.0)
552  vertex.step(gRock, (Zbed - track->pos.z) / track->dir.z);
553  else
554  continue;
555  }
556 
557  if (vertex.getZ() < Zmin) { // propagate muon through water
558  vertex.step(gWater, Zmin - vertex.getZ());
559  }
560 
561  if (vertex.getRange() <= parameters.Dmin_m) {
562  continue;
563  }
564 
565  job.Fill(3.0);
566 
567  const JDetectorSubset_t subdetector(detector, getAxis(*track), maximal_road_width);
568 
569  if (subdetector.empty()) {
570  continue;
571  }
572 
573  job.Fill(4.0);
574 
575  JTrack muon(vertex); // propagate muon trough detector
576 
577  while (vertex.getE() >= parameters.Emin_GeV && vertex.getZ() < Zmax) {
578 
579  const int N = radiation.size();
580 
581  double li[N]; // inverse interaction lengths
582  double ls = 1.0e-5; // minimal total inverse interaction length [m^-1]
583 
584  for (int i = 0; i != N; ++i) {
585  ls += li[i] = radiation[i]->getInverseInteractionLength(vertex.getE());
586  }
587 
588  double As = 0.0; // ionization energy loss
589 
590  for (size_t i = 0; i != ionization.size(); ++i) {
591  As += ionization[i]->getA(vertex.getE());
592  }
593 
594  double step = gRandom->Exp(1.0) / ls; // distance to next radiation process
595  double range = vertex.getRange(As); // range of muon
596 
597  if (vertex.getE() < parameters.Emax_GeV) { // limited step size
598  if (parameters.Dmax_m < range) {
599  range = parameters.Dmax_m;
600  }
601  }
602 
603  double ts = getThetaMCS(vertex.getE(), min(step,range)); // multiple Coulomb scattering angle [rad]
604  double T2 = ts*ts; //
605 
606  rms.Fill(log10(vertex.getE()), (Double_t) 0, ts*ts);
607 
608  vertex.getDirection() += getRandomDirection(T2/3.0); // multiple Coulomb planar scattering
609 
610  vertex.step(As, min(step,range)); // ionization energy loss
611 
612  double Es = 0.0; // shower energy [GeV]
613 
614  if (step < range) {
615 
616  if (vertex.getE() >= parameters.Emin_GeV) {
617 
618  double y = gRandom->Uniform(ls);
619 
620  for (int i = 0; i != N; ++i) {
621 
622  y -= li[i];
623 
624  if (y < 0.0) {
625 
626  Es = radiation[i]->getEnergyOfShower(vertex.getE()); // shower energy [GeV]
627  ts = radiation[i]->getThetaRMS(vertex.getE(), Es); // scattering angle [rad]
628 
629  T2 += ts*ts;
630 
631  rms.Fill(log10(vertex.getE()), (Double_t) radiation[i]->getID(), ts*ts);
632  rad.Fill(log10(vertex.getE()), (Double_t) radiation[i]->getID(), Es);
633 
634  break;
635  }
636  }
637  }
638  }
639 
640  vertex.applyEloss(getRandomDirection(T2), Es);
641 
642  muon.push_back(vertex);
643  }
644 
645  // add muon end point
646 
647  if (vertex.getZ() < Zmax && vertex.getRange() > 0.0) {
648 
649  vertex.step(vertex.getRange());
650 
651  muon.push_back(vertex);
652  }
653 
654  // add information to output muon
655 
656  vector<Trk>::iterator trk = find_if(event.mc_trks.begin(),
657  event.mc_trks.end(),
658  make_predicate(&Trk::id, track->id));
659 
660  if (trk != event.mc_trks.end()) {
661  trk->len = (muon.rbegin()->getZ() < Zmax ? +1 : -1) * (muon.rbegin()->getZ() - muon.begin()->getZ());
662  trk->setusr(mc_usr_keys::energy_lost_in_can, muon.begin()->getE() - muon.rbegin()->getE());
663  }
664 
665  for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
666 
667  const double z0 = muon.begin()->getZ();
668  const double t0 = muon.begin()->getT();
669  const double Z = module->getZ() - module->getX() / getTanThetaC();
670 
671  if (Z >= muon.begin()->getZ() && Z <= muon.rbegin()->getZ()) {
672 
673  const JVector2D pos = muon.getPosition(Z);
674  const double R = hypot(module->getX() - pos.getX(),
675  module->getY() - pos.getY());
676 
677  for (size_t i = 0; i != CDF.size(); ++i) {
678 
679  if (R < CDF[i].integral.getXmax()) {
680 
681  try {
682 
683  double W = 1.0; // mip
684 
685  if (is_deltarays(CDF[i].type)) {
686  W = getDeltaRaysFromMuon(muon.getE(Z)); // delta-rays
687  }
688 
689  const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
690  const int N = gRandom->Poisson(NPE);
691 
692  if (N != 0) {
693 
694  vector<double> npe;
695 
696  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
697 
698  const double R = hypot(pmt->getX() - pos.getX(),
699  pmt->getY() - pos.getY());
700  const double theta = pi.constrain(pmt->getTheta());
701  const double phi = pi.constrain(fabs(pmt->getPhi()));
702 
703  npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
704  }
705 
706  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
707 
708  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
709 
710  const double R = hypot(pmt->getX() - pos.getX(),
711  pmt->getY() - pos.getY());
712  const double Z = pmt->getZ() - z0;
713  const double theta = pi.constrain(pmt->getTheta());
714  const double phi = pi.constrain(fabs(pmt->getPhi()));
715 
716  size_t n0 = min(ns[distance(module->begin(),pmt)], parameters.Nmax_PMT);
717 
718  job.Fill((double) (100 + CDF[i].type), (double) n0);
719 
720  while (n0 != 0) {
721 
722  const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
723  const int n1 = getNumberOfPhotoElectrons(n0);
724 
725  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
726  pmt->getID(),
727  getHitType(CDF[i].type),
728  track->id,
729  t0 + (R * getTanThetaC() + Z) / C + t1,
730  n1));
731 
732  n0 -= n1;
733  }
734  }
735 
736  if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
737  job.Fill((double) (300 + CDF[i].type));
738  }
739  }
740  }
741  catch(const exception& error) {
742  job.Fill((double) (200 + CDF[i].type));
743  }
744  }
745  }
746  }
747  }
748 
749  for (JTrack::const_iterator vertex = muon.begin(); vertex != muon.end(); ++vertex) {
750 
751  const double Es = vertex->getEs();
752 
753  if (Es >= parameters.Ecut_GeV) {
754 
755  const double z0 = vertex->getZ();
756  const double t0 = vertex->getT();
757  const double DZ = geanz.getMaximum(Es);
758 
759  int origin = track->id;
760 
761  if (writeEMShowers) {
762  origin = event.mc_trks.size() + 1;
763  }
764 
765  int number_of_hits = 0;
766 
767  JDetectorSubset_t::range_type range = subdetector.getRange(z0 - maximal_distance,
768  z0 + maximal_distance);
769 
770  for (JDetector::const_iterator module = range.begin(); module != range.end(); ++module) {
771 
772  const double R = hypot(module->getX() - vertex->getX(),
773  module->getY() - vertex->getY());
774  const double Z = module->getZ() - z0 - DZ;
775  const double D = sqrt(R*R + Z*Z);
776  const double cd = Z / D;
777 
778  for (size_t i = 0; i != CDG.size(); ++i) {
779 
780  if (D < CDG[i].integral.getXmax()) {
781 
782  try {
783 
784  const double NPE = CDG[i].integral.getNPE(D, cd) * Es * module->size() * factor;
785  const int N = gRandom->Poisson(NPE);
786 
787  if (N != 0) {
788 
789  vector<double> npe;
790 
791  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
792 
793  const double R = hypot(pmt->getX() - vertex->getX(),
794  pmt->getY() - vertex->getY());
795  const double Z = pmt->getZ() - z0 - DZ;
796  const double D = sqrt(R*R + Z*Z);
797  const double cd = Z / D;
798  const double theta = pi.constrain(pmt->getTheta());
799  const double phi = pi.constrain(fabs(pmt->getPhi()));
800 
801  npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor);
802  }
803 
804  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
805 
806  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
807 
808  const double R = hypot(pmt->getX() - vertex->getX(),
809  pmt->getY() - vertex->getY());
810  const double theta = pi.constrain(pmt->getTheta());
811  const double phi = pi.constrain(fabs(pmt->getPhi()));
812 
813  size_t n0 = min(ns[distance(module->begin(),pmt)], parameters.Nmax_PMT);
814 
815  job.Fill((double) (100 + CDG[i].type), (double) n0);
816 
817  while (n0 != 0) {
818 
819  const double dz = geanz.getLength(Es, gRandom->Rndm());
820  const double Z = pmt->getZ() - z0 - dz;
821  const double D = sqrt(R*R + Z*Z);
822  const double cd = Z / D;
823 
824  const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
825  const int n1 = getNumberOfPhotoElectrons(n0);
826 
827  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
828  pmt->getID(),
829  getHitType(CDG[i].type),
830  origin,
831  t0 + (dz + D * getIndexOfRefraction()) / C + t1,
832  n1));
833 
834  n0 -= n1;
835 
836  number_of_hits += n1;
837  }
838  }
839 
840  if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
841  job.Fill((double) (300 + CDG[i].type));
842  }
843  }
844  }
845  catch(const exception& error) {
846  job.Fill((double) (200 + CDG[i].type));
847  }
848  }
849  }
850  }
851 
852  if (writeEMShowers && number_of_hits != 0) {
853 
854  event.mc_trks.push_back(JTrk_t(origin,
855  TRACK_TYPE_PHOTON,
856  track->id,
857  track->pos + track->dir * vertex->getZ(),
858  track->dir,
859  vertex->getT(),
860  Es));
861  }
862  }
863  }
864 
865  } else if (track->len > 0.0) {
866 
867  // -----------------------------------------------
868  // decayed particles treated as mip (tau includes mip+deltaray)
869  // -----------------------------------------------
870 
871  job.Fill(6.0);
872 
873  const double z0 = 0.0;
874  const double z1 = z0 + track->len;
875  const double t0 = track->t;
876  const double E = track->E;
877 
878  const JTransformation3D transformation = getTransformation(*track);
879 
880  JModule buffer;
881 
882  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
883 
884  const JPosition3D pos = transformation.transform(module->getPosition());
885 
886  const double R = pos.getX();
887  const double Z = pos.getZ() - R / getTanThetaC();
888 
889  if (Z < z0 ||
890  Z > z1 ||
891  R > maximal_road_width) {
892  continue;
893  }
894 
895  for (size_t i = 0; i != CDF.size(); ++i) {
896 
897  double W = 1.0; // mip
898 
899  if (is_deltarays(CDF[i].type)) {
900 
901  if (is_tau(*track))
902  W = getDeltaRaysFromTau(E); // delta-rays
903  else
904  continue;
905  }
906 
907  if (R < CDF[i].integral.getXmax()) {
908 
909  try {
910 
911  const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
912  const int N = gRandom->Poisson(NPE);
913 
914  if (N != 0) {
915 
916  buffer = *module;
917 
918  buffer.transform(transformation);
919 
920  vector<double> npe;
921 
922  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
923 
924  const double R = pmt->getX();
925  const double theta = pi.constrain(pmt->getTheta());
926  const double phi = pi.constrain(fabs(pmt->getPhi()));
927 
928  npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
929  }
930 
931  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
932 
933  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
934 
935  const double R = pmt->getX();
936  const double Z = pmt->getZ() - z0;
937  const double theta = pi.constrain(pmt->getTheta());
938  const double phi = pi.constrain(fabs(pmt->getPhi()));
939 
940  size_t n0 = min(ns[distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
941 
942  job.Fill((double) (120 + CDF[i].type), (double) n0);
943 
944  while (n0 != 0) {
945 
946  const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
947  const int n1 = getNumberOfPhotoElectrons(n0);
948 
949  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
950  pmt->getID(),
951  getHitType(CDF[i].type),
952  track->id,
953  t0 + (R * getTanThetaC() + Z) / C + t1,
954  n1));
955 
956  n0 -= n1;
957  }
958  }
959 
960  if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
961  job.Fill((double) (320 + CDF[i].type));
962  }
963  }
964  }
965  catch(const exception& error) {
966  job.Fill((double) (220 + CDF[i].type));
967  }
968 
969  }
970  }
971  }
972 
973  if (!buffer.empty()) {
974  job.Fill(7.0);
975  }
976 
977  } else if (!is_neutrino(*track)) {
978 
979  if (JPDB::getInstance().hasPDG(track->type)) {
980 
981  // -----------------------------------------------
982  // electron or hadron
983  // -----------------------------------------------
984 
985  job.Fill(8.0);
986 
987  double E = track->E;
988 
989  try {
990  E = getKineticEnergy(E, JPDB::getInstance().getPDG(track->type).mass);
991  }
992  catch(const exception& error) {
993  ERROR(error.what() << endl);
994  }
995 
996  E = pythia(track->type, E);
997 
998  if (E >= parameters.Ecut_GeV && cylinder.getDistance(getPosition(*track)) < parameters.Dmin_m) {
999 
1000  const double z0 = 0.0;
1001  const double t0 = track->t;
1002  const double DZ = geanz.getMaximum(E);
1003 
1004  const JTransformation3D transformation = getTransformation(*track);
1005 
1006  JModule buffer;
1007 
1008  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
1009 
1010  const JPosition3D pos = transformation.transform(module->getPosition());
1011 
1012  const double R = pos.getX();
1013  const double Z = pos.getZ() - z0 - DZ;
1014  const double D = sqrt(R*R + Z*Z);
1015  const double cd = Z / D;
1016 
1017  for (size_t i = 0; i != CDG.size(); ++i) {
1018 
1019  if (D < CDG[i].integral.getXmax()) {
1020 
1021  try {
1022 
1023  const double NPE = CDG[i].integral.getNPE(D, cd) * E * module->size() * factor;
1024  const int N = gRandom->Poisson(NPE);
1025 
1026  if (N != 0) {
1027 
1028  buffer = *module;
1029 
1030  buffer.transform(transformation);
1031 
1032  vector<double> npe;
1033 
1034  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1035 
1036  const double R = pmt->getX();
1037  const double Z = pmt->getZ() - z0 - DZ;
1038  const double D = sqrt(R*R + Z*Z);
1039  const double cd = Z / D;
1040  const double theta = pi.constrain(pmt->getTheta());
1041  const double phi = pi.constrain(fabs(pmt->getPhi()));
1042 
1043  npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * E * factor);
1044  }
1045 
1046  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < parameters.Nmax_NPE);
1047 
1048  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1049 
1050  const double theta = pi.constrain(pmt->getTheta());
1051  const double phi = pi.constrain(fabs(pmt->getPhi()));
1052 
1053  size_t n0 = min(ns[distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
1054 
1055  job.Fill((double) (140 + CDG[i].type), (double) n0);
1056 
1057  while (n0 != 0) {
1058 
1059  const double dz = geanz.getLength(E, gRandom->Rndm());
1060  const double Z = pmt->getZ() - z0 - dz;
1061  const double D = sqrt(R*R + Z*Z);
1062  const double cd = Z / D;
1063 
1064  const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
1065  const int n1 = getNumberOfPhotoElectrons(n0);
1066 
1067  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
1068  pmt->getID(),
1069  getHitType(CDG[i].type, true),
1070  track->id,
1071  t0 + (dz + D * getIndexOfRefraction()) / C + t1,
1072  n1));
1073 
1074  n0 -= n1;
1075  }
1076  }
1077 
1078  if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
1079  job.Fill((double) (340 + CDG[i].type));
1080  }
1081  }
1082  }
1083  catch(const exception& error) {
1084  job.Fill((double) (240 + CDG[i].type));
1085  }
1086  }
1087  }
1088  }
1089 
1090  if (!buffer.empty()) {
1091  job.Fill(9.0);
1092  }
1093 
1094  } else {
1095  job.Fill(21.0);
1096  }
1097  }
1098  }
1099  }
1100 
1101  if (!mc_hits.empty()) {
1102 
1103  mc_hits.merge(parameters.Tmax_ns);
1104 
1105  event.mc_hits.resize(mc_hits.size());
1106 
1107  copy(mc_hits.begin(), mc_hits.end(), event.mc_hits.begin());
1108  }
1109 
1110  timer.stop();
1111 
1112  if (has_neutrino(event)) {
1113  cpu.Fill(log10(get_neutrino(event).E), (double) timer.usec_ucpu * 1.0e-3);
1114  }
1115 
1116  if (event.mc_hits.size() >= numberOfHits) {
1117 
1118  outputFile.put(event);
1119 
1120  job.Fill(10.0);
1121  }
1122  }
1123  STATUS(endl);
1124 
1125  outputFile.put(job);
1126  outputFile.put(cpu);
1127  outputFile.put(rms);
1128  outputFile.put(rad);
1129  outputFile.put(*gRandom);
1130 
1131  JMultipleFileScanner<JTYPELIST<JMetaTypes_t, MultiHead>::typelist> io(inputFile);
1132 
1133  io >> outputFile;
1134 
1135  outputFile.close();
1136 }
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:1514
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 $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
JGEOMETRY2D::JCircle2D
Data structure for circle in two dimensions.
Definition: JCircle2D.hh:33
gmake_property
#define gmake_property(A)
macros 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
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
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:1989
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:172
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:1221
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 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:485
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
JSIRENE::getNumberOfPhotoElectrons
struct JSIRENE::@62 getNumberOfPhotoElectrons
Auxiliary data structure for determination of number of photo-electrons.
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: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
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