Jpp  test_elongated_shower_pde
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/km3-opa_efrac.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, JMeta, 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 
224  properties.insert(gmake_property(numberOfBins));
225  properties.insert(gmake_property(safetyFactor));
226 
227  JParser<> zap("Main program to simulate detector response to muons and showers.");
228 
229  zap['@'] = make_field(properties) = JPARSER::initialised();
230  zap['F'] = make_field(fileDescriptor, "file name descriptor for CDF tables");
231  zap['f'] = make_field(inputFile) = JPARSER::initialised();
232  zap['o'] = make_field(outputFile) = "sirene.root";
233  zap['n'] = make_field(numberOfEvents) = JLimit::max();
234  zap['a'] = make_field(detectorFile) = "";
235  zap['s'] = make_field(writeEMShowers, "store generated EM showers in event");
236  zap['k'] = make_field(keep, "keep position of tracks");
237  zap['N'] = make_field(numberOfHits, "minimum number of hits to output event") = 1;
238  zap['U'] = make_field(factor, "scaling factor applied to light yields") = 1.0;
239  zap['S'] = make_field(seed) = 0;
240  zap['d'] = make_field(debug) = 1;
241 
242  zap(argc, argv);
243  }
244  catch(const exception &error) {
245  FATAL(error.what() << endl);
246  }
247 
248 
249  gRandom->SetSeed(seed);
250 
251 
252  const JMeta meta(argc, argv);
253 
254  const double Zbed = 0.0; // level of seabed [m]
255 
256  vector<JCDF_t> CDF;
257  vector<JCDG_t> CDG;
258 
259  CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_MUON));
260  CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_MUON));
261  CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_DELTARAYS));
262  CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_DELTARAYS));
263 
264  CDG.push_back(JCDG_t(fileDescriptor, DIRECT_LIGHT_FROM_EMSHOWER));
265  CDG.push_back(JCDG_t(fileDescriptor, SCATTERED_LIGHT_FROM_EMSHOWER));
266 
267 
268  double maximal_road_width = 0.0; // road width [m]
269  double maximal_distance = 0.0; // road width [m]
270 
271  for (size_t i = 0; i != CDF.size(); ++i) {
272 
273  DEBUG("Range CDF["<< CDF[i].type << "] " << CDF[i].function.intensity.getXmax() << " m" << endl);
274 
275  maximal_road_width = max(maximal_road_width, CDF[i].function.intensity.getXmax());
276  }
277 
278  for (size_t i = 0; i != CDG.size(); ++i) {
279 
280  DEBUG("Range CDG["<< CDG[i].type << "] " << CDG[i].function.intensity.getXmax() << " m" << endl);
281 
282  if (!is_scattered(CDF[i].type)) {
283  maximal_road_width = max(maximal_road_width, CDG[i].function.intensity.getXmax());
284  }
285 
286  maximal_distance = max(maximal_distance, CDG[i].function.intensity.getXmax());
287  }
288 
289  NOTICE("Maximal road width [m] " << maximal_road_width << endl);
290  NOTICE("Maximal distance [m] " << maximal_distance << endl);
291 
292 
293  if (detectorFile == "" || inputFile.empty()) {
294  STATUS("Nothing to be done." << endl);
295  return 0;
296  }
297 
298  JDetector detector;
299 
300  try {
301 
302  STATUS("Load detector... " << flush);
303 
304  load(detectorFile, detector);
305 
306  STATUS("OK" << endl);
307  }
308  catch(const JException& error) {
309  FATAL(error);
310  }
311 
312  // remove empty modules
313 
314  for (JDetector::iterator module = detector.begin(); module != detector.end(); ) {
315  if (!module->empty())
316  ++module;
317  else
318  module = detector.erase(module);
319  }
320 
321  vector< JSharedPointer<JRadiationInterface> > radiation;
322  vector< JSharedPointer<JACoeffInterface> > ionization;
323 
324  if (true) {
325 
326  STATUS("Setting up radiation tables... " << flush);
327 
328  //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).
329 
330  const JRadiation hydrogen ( 1.0, 1.0, 40, 0.01, 0.1, 0.1);
331  const JRadiation oxygen ( 8.0, 16.0, 40, 0.01, 0.1, 0.1);
332  const JRadiation chlorine (17.0, 35.5, 40, 0.01, 0.1, 0.1);
333  const JRadiation sodium (11.0, 23.0, 40, 0.01, 0.1, 0.1);
334  /* const JRadiation calcium (20.0, 40.0, 40, 0.01, 0.1, 0.1);
335  const JRadiation magnesium(12.0, 24.3, 40, 0.01, 0.1, 0.1);
336  const JRadiation potassium(19.0, 39.0, 40, 0.01, 0.1, 0.1);
337  const JRadiation sulphur (16.0, 32.0, 40, 0.01, 0.1, 0.1);
338  */
339 
340  JSharedPointer<JRadiation> Hydrogen (new JRadiationFunction(hydrogen, 300, 0.2, 1.0e11));
341  JSharedPointer<JRadiation> Oxygen (new JRadiationFunction(oxygen, 300, 0.2, 1.0e11));
342  JSharedPointer<JRadiation> Chlorine (new JRadiationFunction(chlorine, 300, 0.2, 1.0e11));
343  JSharedPointer<JRadiation> Sodium (new JRadiationFunction(sodium, 300, 0.2, 1.0e11));
344  /* JSharedPointer<JRadiation> Calcium (new JRadiationFunction(calcium, 300, 0.2, 1.0e11));
345  JSharedPointer<JRadiation> Magnesium(new JRadiationFunction(magnesium,300, 0.2, 1.0e11));
346  JSharedPointer<JRadiation> Potassium(new JRadiationFunction(potassium,300, 0.2, 1.0e11));
347  JSharedPointer<JRadiation> Sulphur (new JRadiationFunction(sulphur, 300, 0.2, 1.0e11));
348  */
349 
350  radiation.push_back(new JRadiationSource(11, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), EErad_t));
351  radiation.push_back(new JRadiationSource(12, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad_t));
352  radiation.push_back(new JRadiationSource(13, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), EErad_t));
353  radiation.push_back(new JRadiationSource(14, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), EErad_t));
354  /* radiation.push_back(new JRadiationSource(15, Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), EErad_t));
355  radiation.push_back(new JRadiationSource(16, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), EErad_t));
356  radiation.push_back(new JRadiationSource(17, Potassium,DENSITY_SEA_WATER * JSeaWater::K(), EErad_t));
357  radiation.push_back(new JRadiationSource(18, Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), EErad_t));
358  */
359 
360  radiation.push_back(new JRadiationSource(21, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), Brems_t));
361  radiation.push_back(new JRadiationSource(22, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems_t));
362  radiation.push_back(new JRadiationSource(23, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), Brems_t));
363  radiation.push_back(new JRadiationSource(24, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), Brems_t));
364  /* radiation.push_back(new JRadiationSource(25, Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), Brems_t));
365  radiation.push_back(new JRadiationSource(26, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), Brems_t));
366  radiation.push_back(new JRadiationSource(27, Potassium,DENSITY_SEA_WATER * JSeaWater::K(), Brems_t));
367  radiation.push_back(new JRadiationSource(28, Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), Brems_t));
368  */
369 
370  radiation.push_back(new JRadiationSource(31, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), GNrad_t));
371  radiation.push_back(new JRadiationSource(32, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad_t));
372  radiation.push_back(new JRadiationSource(33, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), GNrad_t));
373  radiation.push_back(new JRadiationSource(34, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), GNrad_t));
374  /* radiation.push_back(new JRadiationSource(35, Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), GNrad_t));
375  radiation.push_back(new JRadiationSource(36, Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), GNrad_t));
376  radiation.push_back(new JRadiationSource(37, Potassium,DENSITY_SEA_WATER * JSeaWater::K(), GNrad_t));
377  radiation.push_back(new JRadiationSource(38, Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), GNrad_t));
378  */
379  radiation.push_back(new JDISSource(100, DENSITY_SEA_WATER));
380 
381  ionization.push_back(new JACoeffSource(Oxygen, DENSITY_SEA_WATER * JSeaWater::O()));
382  ionization.push_back(new JACoeffSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl()));
383  ionization.push_back(new JACoeffSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H()));
384  ionization.push_back(new JACoeffSource(Sodium, DENSITY_SEA_WATER * JSeaWater::Na()));
385  /* ionization.push_back(new JACoeffSource(Calcium, DENSITY_SEA_WATER * JSeaWater::Ca()));
386  ionization.push_back(new JACoeffSource(Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg()));
387  ionization.push_back(new JACoeffSource(Potassium,DENSITY_SEA_WATER * JSeaWater::K()));
388  ionization.push_back(new JACoeffSource(Sulphur, DENSITY_SEA_WATER * JSeaWater::S()));
389  */
390 
391  STATUS("OK" << endl);
392  }
393 
394 
395  Vec center(0,0,0);
396  Head header;
397 
398  try {
399 
400  header = inputFile.getHeader();
401 
402  JHead buffer(header);
403 
404  center = get<Vec>(buffer);
405 
406  buffer.simul.push_back(JAANET::simul());
407 
408  buffer.simul.rbegin()->program = APPLICATION_JSIRENE;
409  buffer.simul.rbegin()->version = getGITVersion();
410  buffer.simul.rbegin()->date = getDate();
411  buffer.simul.rbegin()->time = getTime();
412 
413  buffer.detector.push_back(JAANET::detector());
414 
415  buffer.detector.rbegin()->program = APPLICATION_JSIRENE;
416  buffer.detector.rbegin()->filename = detectorFile;
417 
418  buffer.push(&JHead::simul);
419  buffer.push(&JHead::detector);
420 
421  if (!keep) {
422 
423  buffer.coord_origin = coord_origin(0.0, 0.0, 0.0);
424  buffer.can.zmin += center.z;
425  buffer.can.zmax += center.z;
426 
427  buffer.push(&JHead::coord_origin);
428  buffer.push(&JHead::can);
429  }
430 
431  const JCircle2D circle(detector.begin(), detector.end());
432 
433  center += Vec(circle.getX(), circle.getY(), 0.0);
434 
435  copy(buffer, header);
436  }
437  catch(const JException& error) {
438  FATAL(error);
439  }
440 
441  if (!keep)
442  NOTICE("Offset applied to true tracks is: " << center << endl);
443  else
444  NOTICE("No offset applied to true tracks." << endl);
445 
446  JCylinder3D cylinder(detector.begin(), detector.end());
447 
448  cylinder.addMargin(maximal_distance);
449 
450  if (cylinder.getZmin() < Zbed) {
451  cylinder.setZmin(Zbed);
452  }
453 
454  NOTICE("Light generation volume: " << cylinder << endl);
455 
456  TH1D job("job", NULL, 400, 0.5, 400.5);
457  TProfile cpu("cpu", NULL, 16, 0.0, 8.0);
458  TProfile2D rms("rms", NULL, 16, 0.0, 8.0, 201, -0.5, 200.5);
459  TProfile2D rad("rad", NULL, 16, 0.0, 8.0, 201, -0.5, 200.5);
460 
461 
462  outputFile.open();
463 
464  if (!outputFile.is_open()) {
465  FATAL("Error opening file " << outputFile << endl);
466  }
467 
468  outputFile.put(meta);
469  outputFile.put(header);
470  outputFile.put(*gRandom);
471 
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 = pmt->getTheta();
678  const double phi = 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 = pmt->getTheta();
767  const double phi = 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 = pmt->getTheta();
887  const double phi = 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 = pmt->getTheta();
994  const double phi = 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 }
detector
JDetector detector
Definition: JRunAnalyzer.hh:23
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
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:274
JGEOMETRY2D::JCircle2D
Data structure for circle in two dimensions.
Definition: JCircle2D.hh:33
gmake_property
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
Definition: JProperties.hh:1291
safetyFactor
double safetyFactor
safety factor for average CDF integral of optical module
Definition: JSirene.cc:70
JAANET::simul
Generator for simulation.
Definition: JHead.hh: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:41
JAANET::has_neutrino
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Definition: JAAnetToolkit.hh:427
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:330
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:120
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:325
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:84
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:29
JAANET::coord_origin
Coordinate origin.
Definition: JHead.hh:714
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:312
DZ
do echo $TRIPODS[${key}] read X Y Z let DZ
Definition: JFootprint.sh:66
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:30
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
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:238
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:166
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition: JPDFTypes.hh:36
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:40
JSIRENE::JTrack
Muon trajectory.
Definition: JSireneToolkit.hh:216
JGEOMETRY3D::JCylinder3D::addMargin
void addMargin(const double D)
Add (safety) margin.
Definition: JCylinder3D.hh:173
debug
int debug
debug level
Definition: JSirene.cc:68
JPHYSICS::Brems_t
static const JRadiationSource_t Brems_t
Definition: JRadiation.hh:511
JSIRENE::getKineticEnergy
double getKineticEnergy(const double E, const double m)
Get kinetic energy of particle with given mass.
Definition: JSireneToolkit.hh:446
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:1164
JSIRENE::JVertex
Vertex of energy loss of muon.
Definition: JSireneToolkit.hh:303
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:35
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:62
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:193
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
DEBUG
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
Definition: JCalibratePMT.sh:21
JPHYSICS::is_deltarays
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:179
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:438
JAANET::is_tau
bool is_tau(const Trk &track)
Test whether given track is a (anti-)tau.
Definition: JAAnetToolkit.hh:338
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 source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:42
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:46
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19

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