Jpp  17.1.1
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 "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "TRandom3.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/MultiHead.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/definitions/applications.hh"
#include "km3net-dataformat/definitions/trkmembers.hh"
#include "JLang/JSharedPointer.hh"
#include "JLang/Jpp.hh"
#include "JLang/JPredicate.hh"
#include "JSystem/JDate.hh"
#include "JPhysics/JCDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JGeane.hh"
#include "JPhysics/JGeanz.hh"
#include "JPhysics/JRadiation.hh"
#include "JPhysics/JRadiationFunction.hh"
#include "JPhysics/JRadiationSource.hh"
#include "JPhysics/JACoeffSource.hh"
#include "JPhysics/JPhysicsSupportkit.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JPDB.hh"
#include "JSirene/JSirene.hh"
#include "JSirene/pythia.hh"
#include "JSirene/JSeaWater.hh"
#include "JSirene/JCDFTable1D.hh"
#include "JSirene/JCDFTable2D.hh"
#include "JSirene/JSireneToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSubset.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Variables

int debug
 debug level More...
 
int numberOfBins = 200
 number of bins for average CDF integral of optical module More...
 
double safetyFactor = 1.7
 safety factor for average CDF integral of optical module More...
 

Detailed Description

Main program to simulate detector response to muons and showers.

sirene.png
Picture by Claudine Colnard

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

The CDF tables can be produced with the script JMakePDF.sh: 

JMakePDF.sh -P

(option -h will print all available options). Note that the script will launch a large number of processes (JMakePDF and JMakePDG)
which may take a considerable amount of time to completion.
On a standard desktop, all jobs should be finished within 1/2 a day or so.

The same script should then be run with option -M to merge the PDF files, i.e: 

JMakePDF.sh -M

CDF tables are obtained by running the same script with option -C, i.e: 

JMakePDF.sh -C

The various PDFs can be drawn using the JDrawPDF or JDrawPDG applications.
The tabulated PDFs can be plotted using the JPlotPDF or JPlotPDG applications.
The tabulated CDFs can be plotted using the JPlotCDF or JPlotCDG applications.

Author
mdejong

Definition in file JSirene.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 194 of file JSirene.cc.

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

Variable Documentation

int debug

debug level

Definition at line 67 of file JSirene.cc.

int numberOfBins = 200

number of bins for average CDF integral of optical module

Definition at line 68 of file JSirene.cc.

double safetyFactor = 1.7

safety factor for average CDF integral of optical module

Definition at line 69 of file JSirene.cc.

Generated by   doxygen 1.8.5