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

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Variables

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

Detailed Description

Main program to simulate detector response to muons and showers.

sirene.png
Picture by Claudine Colnard

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

More accuracy can be achieved by setting compile option RADITION but it will run slower.

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

JMakePDF.sh -P

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

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

JMakePDF.sh -M

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

JMakePDF.sh -C

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

Author
mdejong

Definition in file JSirene.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 198 of file JSirene.cc.

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

Variable Documentation

int debug

debug level

Definition at line 69 of file JSirene.cc.

int numberOfBins = 200

number of bins for average CDF integral of optical module

Definition at line 70 of file JSirene.cc.

double safetyFactor = 1.7

safety factor for average CDF integral of optical module

Definition at line 71 of file JSirene.cc.

Generated by   doxygen 1.8.5