Jpp  17.0.0-rc.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/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 193 of file JSirene.cc.

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

Variable Documentation

int debug

debug level

Definition at line 66 of file JSirene.cc.

int numberOfBins = 200

number of bins for average CDF integral of optical module

Definition at line 67 of file JSirene.cc.

double safetyFactor = 1.7

safety factor for average CDF integral of optical module

Definition at line 68 of file JSirene.cc.

Generated by   doxygen 1.8.5