Jpp  15.0.0-rc.2
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 "TProfile.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/JRadiationSource.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JAAnet/JHit_t.hh"
#include "JAAnet/JTrk_t.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JPDB.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 175 of file JSirene.cc.

176 {
177  using namespace std;
178  using namespace JPP;
179 
180  string fileDescriptor;
183  JLimit_t& numberOfEvents = inputFile.getLimit();
184  string detectorFile;
185  double Ecut_GeV = 0.1; // minimal energy for generation of light from shower [GeV]
186  double Emin_GeV = 1.0; // minimal energy of muon for shower generation [GeV]
187  double Dmin_m = 0.1; // minimal distance for positioning [m]
188  double Tmax_ns = 0.1; // maximal time between hits on same PMT to be merged
189  int Nmax_npe = numeric_limits<int>::max(); // maximal number of photo-electrons
190  bool writeEMShowers;
191  bool keep;
192  size_t numberOfHits;
193  double factor;
194  UInt_t seed;
195 
196  try {
197 
198  JProperties properties;
199 
200  properties.insert(gmake_property(Ecut_GeV));
201  properties.insert(gmake_property(Emin_GeV));
202  properties.insert(gmake_property(Dmin_m));
203  properties.insert(gmake_property(Tmax_ns));
204  properties.insert(gmake_property(Nmax_npe));
205  properties.insert(gmake_property(numberOfBins));
206  properties.insert(gmake_property(safetyFactor));
207 
208  JParser<> zap("Main program to simulate detector response to muons and showers.");
209 
210  zap['@'] = make_field(properties) = JPARSER::initialised();
211  zap['F'] = make_field(fileDescriptor, "file name descriptor for CDF tables");
212  zap['f'] = make_field(inputFile) = JPARSER::initialised();
213  zap['o'] = make_field(outputFile) = "sirene.root";
214  zap['n'] = make_field(numberOfEvents) = JLimit::max();
215  zap['a'] = make_field(detectorFile) = "";
216  zap['s'] = make_field(writeEMShowers, "store generated EM showers in event");
217  zap['k'] = make_field(keep, "keep position of tracks");
218  zap['N'] = make_field(numberOfHits, "minimum number of hits to output event") = 1;
219  zap['U'] = make_field(factor, "scaling factor applied to light yields") = 1.0;
220  zap['S'] = make_field(seed) = 0;
221  zap['d'] = make_field(debug) = 1;
222 
223  zap(argc, argv);
224  }
225  catch(const exception &error) {
226  FATAL(error.what() << endl);
227  }
228 
229 
230  gRandom->SetSeed(seed);
231 
232 
233  const JMeta meta(argc, argv);
234 
235  const double Zbed = 0.0; // level of seabed [m]
236 
238  vector<JCDG_t> CDG;
239 
240  CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_MUON));
241  CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_MUON));
242  CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_DELTARAYS));
243  CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_DELTARAYS));
244 
245  CDG.push_back(JCDG_t(fileDescriptor, DIRECT_LIGHT_FROM_EMSHOWER));
246  CDG.push_back(JCDG_t(fileDescriptor, SCATTERED_LIGHT_FROM_EMSHOWER));
247 
248 
249  double maximal_road_width = 0.0; // road width [m]
250  double maximal_distance = 0.0; // road width [m]
251 
252  for (size_t i = 0; i != CDF.size(); ++i) {
253 
254  DEBUG("Range CDF["<< CDF[i].type << "] " << CDF[i].function.intensity.getXmax() << " m" << endl);
255 
256  maximal_road_width = max(maximal_road_width, CDF[i].function.intensity.getXmax());
257  }
258 
259  for (size_t i = 0; i != CDG.size(); ++i) {
260 
261  DEBUG("Range CDG["<< CDG[i].type << "] " << CDG[i].function.intensity.getXmax() << " m" << endl);
262 
263  if (!is_scattered(CDF[i].type)) {
264  maximal_road_width = max(maximal_road_width, CDG[i].function.intensity.getXmax());
265  }
266 
267  maximal_distance = max(maximal_distance, CDG[i].function.intensity.getXmax());
268  }
269 
270  NOTICE("Maximal road width [m] " << maximal_road_width << endl);
271  NOTICE("Maximal distance [m] " << maximal_distance << endl);
272 
273 
274  if (detectorFile == "" || inputFile.empty()) {
275  STATUS("Nothing to be done." << endl);
276  return 0;
277  }
278 
280 
281  try {
282 
283  STATUS("Load detector... " << flush);
284 
285  load(detectorFile, detector);
286 
287  STATUS("OK" << endl);
288  }
289  catch(const JException& error) {
290  FATAL(error);
291  }
292 
293  // remove empty modules
294 
295  for (JDetector::iterator module = detector.begin(); module != detector.end(); ) {
296  if (!module->empty())
297  ++module;
298  else
299  module = detector.erase(module);
300  }
301 
302  vector< JSharedPointer<JRadiationInterface> > radiation, ionization;
303 
304  if (true) {
305 
306  STATUS("Setting up radiation tables... " << flush);
307 
308  //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).
309 
310  const JRadiation hydrogen ( 1.0, 1.0, 40, 0.01, 0.1, 0.1);
311  const JRadiation oxygen ( 8.0, 16.0, 40, 0.01, 0.1, 0.1);
312  const JRadiation chlorine (17.0, 35.5, 40, 0.01, 0.1, 0.1);
313  const JRadiation sodium (11.0, 23.0, 40, 0.01, 0.1, 0.1);
314 /* const JRadiation calcium (20.0, 40.0, 40, 0.01, 0.1, 0.1);
315  const JRadiation magnesium(12.0, 24.3, 40, 0.01, 0.1, 0.1);
316  const JRadiation potassium(19.0, 39.0, 40, 0.01, 0.1, 0.1);
317  const JRadiation sulphur (16.0, 32.0, 40, 0.01, 0.1, 0.1);
318 */
319 
320  JSharedPointer<JRadiation> Hydrogen (new JRadiationFunction(hydrogen, 300, 0.2, 1.0e11));
321  JSharedPointer<JRadiation> Oxygen (new JRadiationFunction(oxygen, 300, 0.2, 1.0e11));
322  JSharedPointer<JRadiation> Chlorine (new JRadiationFunction(chlorine, 300, 0.2, 1.0e11));
323  JSharedPointer<JRadiation> Sodium (new JRadiationFunction(sodium, 300, 0.2, 1.0e11));
324 /* JSharedPointer<JRadiation> Calcium (new JRadiationFunction(calcium, 300, 0.2, 1.0e11));
325  JSharedPointer<JRadiation> Magnesium(new JRadiationFunction(magnesium,300, 0.2, 1.0e11));
326  JSharedPointer<JRadiation> Potassium(new JRadiationFunction(potassium,300, 0.2, 1.0e11));
327  JSharedPointer<JRadiation> Sulphur (new JRadiationFunction(sulphur, 300, 0.2, 1.0e11));
328 */
329 
330  JRadiationSource::source_type EErad(&JRadiation::TotalCrossSectionEErad, &JRadiation::EfromEErad);
331  JRadiationSource::source_type Brems(&JRadiation::TotalCrossSectionBrems, &JRadiation::EfromBrems);
332  JRadiationSource::source_type GNrad(&JRadiation::TotalCrossSectionGNrad, &JRadiation::EfromGNrad);
333  JRadiationSource::source_type Ion (&JRadiation::CalculateACoeff , NULL);
334 
335  radiation.push_back(new JRadiationSource(Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), EErad));
336  radiation.push_back(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad));
337  radiation.push_back(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), EErad));
338  radiation.push_back(new JRadiationSource(Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), EErad));
339 /* radiation.push_back(new JRadiationSource(Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), EErad));
340  radiation.push_back(new JRadiationSource(Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), EErad));
341  radiation.push_back(new JRadiationSource(Potassium,DENSITY_SEA_WATER * JSeaWater::K(), EErad));
342  radiation.push_back(new JRadiationSource(Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), EErad));
343 */
344 
345  radiation.push_back(new JRadiationSource(Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), Brems));
346  radiation.push_back(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems));
347  radiation.push_back(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), Brems));
348  radiation.push_back(new JRadiationSource(Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), Brems));
349 /* radiation.push_back(new JRadiationSource(Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), Brems));
350  radiation.push_back(new JRadiationSource(Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), Brems));
351  radiation.push_back(new JRadiationSource(Potassium,DENSITY_SEA_WATER * JSeaWater::K(), Brems));
352  radiation.push_back(new JRadiationSource(Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), Brems));
353 */
354 
355  radiation.push_back(new JRadiationSource(Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), GNrad));
356  radiation.push_back(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad));
357  radiation.push_back(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), GNrad));
358  radiation.push_back(new JRadiationSource(Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), GNrad));
359 /* radiation.push_back(new JRadiationSource(Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(), GNrad));
360  radiation.push_back(new JRadiationSource(Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(), GNrad));
361  radiation.push_back(new JRadiationSource(Potassium,DENSITY_SEA_WATER * JSeaWater::K(), GNrad));
362  radiation.push_back(new JRadiationSource(Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), GNrad));
363 */
364  ionization.push_back(new JRadiationSource(Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), Ion));
365  ionization.push_back(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(),Ion));
366  ionization.push_back(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), Ion));
367  ionization.push_back(new JRadiationSource(Sodium, DENSITY_SEA_WATER * JSeaWater::Na(),Ion));
368 /* ionization.push_back(new JRadiationSource(Calcium, DENSITY_SEA_WATER * JSeaWater::Ca(),Ion));
369  ionization.push_back(new JRadiationSource(Magnesium,DENSITY_SEA_WATER * JSeaWater::Mg(),Ion));
370  ionization.push_back(new JRadiationSource(Potassium,DENSITY_SEA_WATER * JSeaWater::K(), Ion));
371  ionization.push_back(new JRadiationSource(Sulphur, DENSITY_SEA_WATER * JSeaWater::S(), Ion));
372 */
373 
374  STATUS("OK" << endl);
375  }
376 
377 
378  Vec center(0,0,0);
379  Head header;
380 
381  try {
382 
383  header = inputFile.getHeader();
384 
385  JHead buffer(header);
386 
387  center = get<Vec>(buffer);
388 
389  buffer.simul.push_back(JAANET::simul());
390 
391  buffer.simul.rbegin()->program = APPLICATION_JSIRENE;
392  buffer.simul.rbegin()->version = getGITVersion();
393  buffer.simul.rbegin()->date = getDate();
394  buffer.simul.rbegin()->time = getTime();
395 
396  buffer.detector.push_back(JAANET::detector());
397 
398  buffer.detector.rbegin()->program = APPLICATION_JSIRENE;
399  buffer.detector.rbegin()->filename = detectorFile;
400 
401  buffer.push(&JHead::simul);
402  buffer.push(&JHead::detector);
403 
404  if (!keep) {
405 
406  buffer.coord_origin = coord_origin(0.0, 0.0, 0.0);
407  buffer.can.zmin += center.z;
408  buffer.can.zmax += center.z;
409 
410  buffer.push(&JHead::coord_origin);
411  buffer.push(&JHead::can);
412  }
413 
414  const JCircle2D circle(detector.begin(), detector.end());
415 
416  center += Vec(circle.getX(), circle.getY(), 0.0);
417 
418  copy(buffer, header);
419  }
420  catch(const JException& error) {
421  FATAL(error);
422  }
423 
424  if (!keep)
425  NOTICE("Offset applied to true tracks is: " << center << endl);
426  else
427  NOTICE("No offset applied to true tracks." << endl);
428 
429  JCylinder3D cylinder(detector.begin(), detector.end());
430 
431  cylinder.addMargin(maximal_distance);
432 
433  if (cylinder.getZmin() < Zbed) {
434  cylinder.setZmin(Zbed);
435  }
436 
437  NOTICE("Light generation volume: " << cylinder << endl);
438 
439  TProfile cpu("cpu", NULL, 14, 1.0, 8.0);
440  TH1D job("job", NULL, 400, 0.5, 400.5);
441 
442 
443  outputFile.open();
444 
445  if (!outputFile.is_open()) {
446  FATAL("Error opening file " << outputFile << endl);
447  }
448 
449  outputFile.put(meta);
450  outputFile.put(header);
451  outputFile.put(*gRandom);
452 
453 
454  JTimer timer;
455 
456  for (JMultipleFileScanner<Evt>& in = inputFile; in.hasNext(); ) {
457 
458  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
459 
460  job.Fill(1.0);
461 
462  Evt* evt = in.next();
463 
464  if (!keep) {
465  for (vector<Trk>::iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
466  track->pos += center;
467  }
468  }
469 
470  Evt event(*evt); // output
471 
472  event.mc_hits.clear();
473 
474  JHits_t mc_hits; // temporary buffer
475 
476  timer.reset();
477  timer.start();
478 
479  for (vector<Trk>::const_iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
480 
481  if (!track->is_finalstate()) continue; //only final state particles produce light
482 
483  if (is_muon(*track)) {
484 
485  // -----------------------------------------------
486  // muon
487  // -----------------------------------------------
488 
489  job.Fill(2.0);
490 
492 
493  const JCylinder3D::intersection_type intersection = cylinder.getIntersection(getAxis(*track));
494 
495  double Zmin = intersection.first;
496  double Zmax = intersection.second;
497 
498  if (Zmax - Zmin <= Dmin_m) {
499  continue;
500  }
501 
502  JVertex vertex(0.0, track->t, track->E); // start of muon
503 
504  if (track->pos.z < Zbed) { // propagate muon through rock
505 
506  if (track->dir.z > 0.0)
507  vertex.step(gRock, (Zbed - track->pos.z) / track->dir.z);
508  else
509  continue;
510  }
511 
512  if (vertex.getZ() < Zmin) { // propagate muon through water
513  vertex.step(gWater, Zmin - vertex.getZ());
514  }
515 
516  if (vertex.getRange() <= Dmin_m) {
517  continue;
518  }
519 
520  job.Fill(3.0);
521 
522  const JDetectorSubset_t subdetector(detector, getAxis(*track), maximal_road_width);
523 
524  if (subdetector.empty()) {
525  continue;
526  }
527 
528  job.Fill(4.0);
529 
530  JTrack muon(vertex); // propagate muon trough detector
531 
532  while (vertex.getE() >= Emin_GeV && vertex.getZ() < Zmax) {
533 
534  const int N = radiation.size();
535 
536  double li[N]; // inverse interaction lengths
537  double ls = 1.0e-5; // minimal total inverse interaction length (100 km)^-1
538 
539  for (int i = 0; i != N; ++i) {
540  ls += li[i] = radiation[i]->getInverseInteractionLength(vertex.getE());
541  }
542 
543  double As = 0.0; // ionization energy loss
544 
545  for (size_t i = 0; i != ionization.size(); ++i) {
546  As += ionization[i]->getA(vertex.getE());
547  }
548 
549  const double ds = min(gRandom->Exp(1.0) / ls, vertex.getRange(As));
550 
551  vertex.step(As, ds); // ionization energy loss
552 
553  if (vertex.getE() >= Emin_GeV) {
554 
555  double Es = 0.0; // shower energy [GeV]
556  double y = gRandom->Uniform(ls);
557 
558  for (int i = 0; i != N; ++i) {
559 
560  y -= li[i];
561 
562  if (y < 0.0) {
563  Es = radiation[i]->getEnergyOfShower(vertex.getE()); // shower energy [GeV]
564  break;
565  }
566  }
567 
568  vertex.applyEloss(Es);
569 
570  if (Es >= Ecut_GeV) {
571  muon.push_back(vertex);
572  }
573  }
574  }
575 
576  // add muon end point
577 
578  if (vertex.getE() < Emin_GeV && vertex.getRange() > 0.0) {
579 
580  vertex.step(vertex.getRange()); // for energies below minimum, default ionisation energy loss is used.
581 
582  muon.push_back(vertex);
583  }
584 
585  // add information to output muon
586 
587  vector<Trk>::iterator trk = find_if(event.mc_trks.begin(),
588  event.mc_trks.end(),
589  make_predicate(&Trk::id, track->id));
590 
591  if (trk != event.mc_trks.end()) {
592  trk->len = (muon.rbegin()->getZ() < Zmax ? +1 : -1) * (muon.rbegin()->getZ() - muon.begin()->getZ());
593  trk->setusr(mc_usr_keys::energy_lost_in_can, muon.begin()->getE() - muon.rbegin()->getE());
594  }
595 
596  for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
597 
598  const double z0 = muon.begin()->getZ();
599  const double t0 = muon.begin()->getT();
600  const double R = module->getX();
601  const double Z = module->getZ() - R / getTanThetaC();
602 
603  if (Z >= muon.begin()->getZ() && Z <= muon.rbegin()->getZ()) {
604 
605  for (size_t i = 0; i != CDF.size(); ++i) {
606 
607  if (R < CDF[i].integral.getXmax()) {
608 
609  try {
610 
611  double W = 1.0; // mip
612 
613  if (is_deltarays(CDF[i].type)) {
614  W = getDeltaRaysFromMuon(muon.getE(Z)); // delta-rays
615  }
616 
617  const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
618  const int N = gRandom->Poisson(NPE);
619 
620  if (N != 0) {
621 
622  double ns = 0.0;
623 
624  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
625 
626  const double R = pmt->getX();
627  const double Z = pmt->getZ() - z0;
628  const double theta = pmt->getTheta();
629  const double phi = fabs(pmt->getPhi());
630 
631  const double npe = CDF[i].function.getNPE(R, theta, phi) * factor * W;
632 
633  ns += npe;
634 
635  int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), Nmax_npe);
636 
637  job.Fill((double) (100 + CDF[i].type), (double) n0);
638 
639  while (n0 != 0) {
640 
641  const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
642  const int n1 = getNumberOfPhotoElectrons(n0);
643 
644  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
645  pmt->getID(),
646  getHitType(CDF[i].type),
647  track->id,
648  t0 + (R * getTanThetaC() + Z) / C + t1,
649  n1));
650 
651  n0 -= n1;
652  }
653  }
654 
655  if (ns > NPE) {
656  job.Fill((double) (300 + CDF[i].type));
657  }
658  }
659  }
660  catch(const exception& error) {
661  job.Fill((double) (200 + CDF[i].type));
662  }
663  }
664  }
665  }
666  }
667 
668  for (JTrack::const_iterator vertex = muon.begin(); vertex != muon.end(); ++vertex) {
669 
670  const double z0 = vertex->getZ();
671  const double t0 = vertex->getT();
672  const double Es = vertex->getEs();
673 
674  int origin = track->id;
675 
676  if (writeEMShowers) {
677  origin = event.mc_trks.size() + 1;
678  }
679 
680  int number_of_hits = 0;
681 
682  JDetectorSubset_t::range_type range = subdetector.getRange(z0 - maximal_distance,
683  z0 + maximal_distance);
684 
685  for (JDetector::const_iterator module = range.begin(); module != range.end(); ++module) {
686 
687  const double dz = geanz.getMaximum(Es);
688  const double R = module->getX();
689  const double Z = module->getZ() - z0 - dz;
690  const double D = sqrt(R*R + Z*Z);
691  const double cd = Z / D;
692 
693  for (size_t i = 0; i != CDG.size(); ++i) {
694 
695  if (D < CDG[i].integral.getXmax()) {
696 
697  try {
698 
699  const double NPE = CDG[i].integral.getNPE(D, cd) * Es * module->size() * factor;
700  const int N = gRandom->Poisson(NPE);
701 
702  if (N != 0) {
703 
704  double ns = 0.0;
705 
706  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
707 
708  const double dz = geanz.getMaximum(Es);
709  const double R = pmt->getX();
710  const double Z = pmt->getZ() - z0 - dz;
711  const double theta = pmt->getTheta();
712  const double phi = fabs(pmt->getPhi());
713  const double D = sqrt(R*R + Z*Z);
714  const double cd = Z / D;
715 
716  const double npe = CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor;
717 
718  ns += npe;
719 
720  int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), Nmax_npe);
721 
722  job.Fill((double) (100 + CDG[i].type), (double) n0);
723 
724  while (n0 != 0) {
725 
726  const double dz = geanz.getLength(Es, gRandom->Rndm());
727  const double Z = pmt->getZ() - z0 - dz;
728  const double D = sqrt(R*R + Z*Z);
729  const double cd = Z / D;
730 
731  const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
732  const int n1 = getNumberOfPhotoElectrons(n0);
733 
734  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
735  pmt->getID(),
736  getHitType(CDG[i].type),
737  origin,
738  t0 + (dz + D * getIndexOfRefraction()) / C + t1,
739  n1));
740 
741  n0 -= n1;
742 
743  number_of_hits += n1;
744  }
745  }
746 
747  if (ns > NPE) {
748  job.Fill((double) (300 + CDG[i].type));
749  }
750  }
751  }
752  catch(const exception& error) {
753  job.Fill((double) (200 + CDG[i].type));
754  }
755  }
756  }
757  }
758 
759  if (writeEMShowers && number_of_hits != 0) {
760 
761  event.mc_trks.push_back(JTrk_t(origin,
763  track->id,
764  track->pos + track->dir * vertex->getZ(),
765  track->dir,
766  vertex->getT(),
767  Es));
768  }
769  }
770 
771  } else if (track->len>0) {
772 
773  // -----------------------------------------------
774  // decayed particles treated as mip (tau includes mip+deltaray)
775  // -----------------------------------------------
776 
777  job.Fill(6.0);
778 
779  const double z0 = 0.0;
780  const double z1 = z0 + track->len;
781  const double t0 = track->t;
782  const double E = track->E;
783 
784  const JTransformation3D transformation = getTransformation(*track);
785 
786  JModule buffer;
787 
788  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
789 
790  const JPosition3D pos = transformation.transform(module->getPosition());
791 
792  const double R = pos.getX();
793  const double Z = pos.getZ() - R / getTanThetaC();
794 
795  if (Z < z0 ||
796  Z > z1 ||
797  R > maximal_road_width) {
798  continue;
799  }
800 
801  for (size_t i = 0; i != CDF.size(); ++i) {
802 
803  double W = 1.0; // mip
804 
805  if (is_deltarays(CDF[i].type)) {
806  if (is_tau(*track)) W = getDeltaRaysFromTau(E); // delta-rays
807  else continue;
808  }
809 
810  if (R < CDF[i].integral.getXmax()) {
811 
812  try {
813 
814  const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
815  const int N = gRandom->Poisson(NPE);
816 
817  if (N != 0) {
818 
819  buffer = *module;
820 
821  buffer.transform(transformation);
822 
823  double ns = 0.0;
824 
825  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
826 
827  const double R = pmt->getX();
828  const double Z = pmt->getZ() - z0;
829  const double theta = pmt->getTheta();
830  const double phi = fabs(pmt->getPhi());
831 
832  const double npe = CDF[i].function.getNPE(R, theta, phi) * factor * W;
833 
834  ns += npe;
835 
836  int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), Nmax_npe);
837 
838  job.Fill((double) (120 + CDF[i].type), (double) n0);
839 
840  while (n0 != 0) {
841 
842  const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
843  const int n1 = getNumberOfPhotoElectrons(n0);
844 
845  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
846  pmt->getID(),
847  getHitType(CDF[i].type),
848  track->id,
849  t0 + (R * getTanThetaC() + Z) / C + t1,
850  n1));
851 
852  n0 -= n1;
853  }
854  }
855 
856  if (ns > NPE) {
857  job.Fill((double) (320 + CDF[i].type));
858  }
859  }
860  }
861  catch(const exception& error) {
862  job.Fill((double) (220 + CDF[i].type));
863  }
864 
865  }
866  }
867  }
868 
869  if (!buffer.empty()) {
870  job.Fill(7.0);
871  }
872 
873  } else if (!is_neutrino(*track)) {
874 
875  if (JPDB::getInstance().hasPDG(track->type)) {
876 
877  // -----------------------------------------------
878  // electron or hadron
879  // -----------------------------------------------
880 
881  job.Fill(8.0);
882 
883  double E = track->E;
884 
885  try {
886  E = getKineticEnergy(E, JPDB::getInstance().getPDG(track->type).mass);
887  }
888  catch(const exception& error) {
889  ERROR(error.what() << endl);
890  }
891 
892  E = pythia(track->type, E);
893 
894  if (E < Ecut_GeV || cylinder.getDistance(getPosition(*track)) > Dmin_m) {
895  continue;
896  }
897 
898  const double z0 = 0.0;
899  const double t0 = track->t;
900 
901  const JTransformation3D transformation = getTransformation(*track);
902 
903  JModule buffer;
904 
905  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
906 
907  const JPosition3D pos = transformation.transform(module->getPosition());
908 
909  const double dz = geanz.getMaximum(E);
910  const double R = pos.getX();
911  const double Z = pos.getZ() - z0 - dz;
912  const double D = sqrt(R*R + Z*Z);
913  const double cd = Z / D;
914 
915  if (D > maximal_distance) {
916  continue;
917  }
918 
919  for (size_t i = 0; i != CDG.size(); ++i) {
920 
921  if (D < CDG[i].integral.getXmax()) {
922 
923  try {
924 
925  const double NPE = CDG[i].integral.getNPE(D, cd) * E * module->size() * factor;
926  const int N = gRandom->Poisson(NPE);
927 
928  if (N != 0) {
929 
930  buffer = *module;
931 
932  buffer.transform(transformation);
933 
934  double ns = 0.0;
935 
936  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
937 
938  const double dz = geanz.getMaximum(E);
939  const double R = pmt->getX();
940  const double Z = pmt->getZ() - z0 - dz;
941  const double theta = pmt->getTheta();
942  const double phi = fabs(pmt->getPhi());
943  const double D = sqrt(R*R + Z*Z);
944  const double cd = Z / D;
945 
946  const double npe = CDG[i].function.getNPE(D, cd, theta, phi) * E * factor;
947 
948  ns += npe;
949 
950  int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), Nmax_npe);
951 
952  job.Fill((double) (140 + CDG[i].type), (double) n0);
953 
954  while (n0 != 0) {
955 
956  const double dz = geanz.getLength(E, gRandom->Rndm());
957  const double Z = pmt->getZ() - z0 - dz;
958  const double D = sqrt(R*R + Z*Z);
959  const double cd = Z / D;
960 
961  const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
962  const int n1 = getNumberOfPhotoElectrons(n0);
963 
964  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
965  pmt->getID(),
966  getHitType(CDG[i].type, true),
967  track->id,
968  t0 + (dz + D * getIndexOfRefraction()) / C + t1,
969  n1));
970 
971  n0 -= n1;
972  }
973  }
974 
975  if (ns > NPE) {
976  job.Fill((double) (340 + CDG[i].type));
977  }
978  }
979  }
980  catch(const exception& error) {
981  job.Fill((double) (240 + CDG[i].type));
982  }
983  }
984  }
985  }
986 
987  if (!buffer.empty()) {
988  job.Fill(9.0);
989  }
990 
991  } else {
992  job.Fill(21.0);
993  }
994  }
995  }
996 
997  if (!mc_hits.empty()) {
998 
999  mc_hits.merge(Tmax_ns);
1000 
1001  event.mc_hits.resize(mc_hits.size());
1002 
1003  copy(mc_hits.begin(), mc_hits.end(), event.mc_hits.begin());
1004  }
1005 
1006  timer.stop();
1007 
1008  if (has_neutrino(event)) {
1009  cpu.Fill(log10(get_neutrino(event).E), (double) timer.usec_ucpu * 1.0e-3);
1010  }
1011 
1012  if (event.mc_hits.size() >= numberOfHits) {
1013 
1014  outputFile.put(event);
1015 
1016  job.Fill(10.0);
1017  }
1018  }
1019  STATUS(endl);
1020 
1021  outputFile.put(cpu);
1022  outputFile.put(job);
1023  outputFile.put(*gRandom);
1024  outputFile.close();
1025 }
Auxiliary class to set-up Hit.
Definition: JHit_t.hh:25
const char *const energy_lost_in_can
Definition: io_ascii.hh:45
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
Object writing to file.
then usage $script[detector file(input file)+[output file[CDF file descriptor]]] nAuxiliary script for simulation of detector response nNote that if more than one input file is all other arguments must be provided fi case set_variable CDF
Definition: JSirene.sh:34
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
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
Data structure for a composite optical module.
Definition: JModule.hh:68
JTransformation3D getTransformation(const Trk &track)
Get transformation.
Data structure for circle in two dimensions.
Definition: JCircle2D.hh:33
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
double safetyFactor
safety factor for average CDF integral of optical module
Definition: JSirene.cc:65
Generator for simulation.
Definition: JHead.hh:460
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
#define STATUS(A)
Definition: JMessage.hh:63
scattered light from EM shower
Definition: JPDFTypes.hh:41
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Detector data structure.
Definition: JDetector.hh:89
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
JHitType_t getHitType(const JPDFType_t pdf, const bool shower=false)
Get hit type corresponding to given PDF type.
double getDeltaRaysFromTau(const double E)
Equivalent EM-shower energy due to delta-rays per unit tau track length.
Definition: JPDFToolkit.hh:95
Utility class to parse parameter values.
Definition: JProperties.hh:496
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
static const double H
Planck constant [eV s].
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
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.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:323
void merge(const double Tmax_ns)
Merge hits on same PMT that are within given time window.
Definition: JHit_t.hh:116
static const double DENSITY_SEA_WATER
Fixed environment values.
static const JPythia pythia
Function object for relative light yield as a function of GEANT particle code.
Definition: JPythia.hh:96
double getTime(const Hit &hit)
Get true time of hit.
string outputFile
then usage E
Definition: JMuonPostfit.sh:35
static const char *const APPLICATION_JSIRENE
detector simulation
Definition: applications.hh:19
direct light from muon
Definition: JPDFTypes.hh:29
Coordinate origin.
Definition: JHead.hh:663
Fast implementation of class JRadiation.
Source of radiation.
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
static const double C
Physics constants.
scattered light from muon
Definition: JPDFTypes.hh:30
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition: Vec.hh:12
Cylinder object.
Definition: JCylinder3D.hh:39
Detector file.
Definition: JHead.hh:196
std::string getGITVersion(const std::string &tag)
Get GIT version for given GIT tag.
The template JSharedPointer class can be used to share a pointer to an object.
Auxiliary class for the calculation of the muon radiative cross sections.
Definition: JRadiation.hh:31
T & getInstance(const T &object)
Get static instance from temporary object.
Definition: JObject.hh:75
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
JAxis3D getAxis(const Trk &track)
Get axis.
Auxiliary class for CPU timing and usage.
Definition: JTimer.hh:32
JPosition3D getPosition(const Vec &pos)
Get position.
scattered light from delta-rays
Definition: JPDFTypes.hh:36
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
direct light from EM shower
Definition: JPDFTypes.hh:40
Muon trajectory.
Auxiliary data structure for list of hits with hit merging capability.
Definition: JHit_t.hh:105
void addMargin(const double D)
Add (safety) margin.
Definition: JCylinder3D.hh:173
int debug
debug level
Definition: JSirene.cc:63
double getKineticEnergy(const double E, const double m)
Get kinetic energy of particle with given mass.
const char * getDate(const JDateAndTimeFormat option=ISO8601)
Get ASCII formatted date.
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:43
Detector subset with binary search functionality.
Detector subset without binary search functionality.
Monte Carlo run header.
Definition: JHead.hh:1113
Vertex of energy loss of muon.
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:67
#define FATAL(A)
Definition: JMessage.hh:67
int id
track identifier
Definition: Trk.hh:16
direct light from delta-rays
Definition: JPDFTypes.hh:35
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
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 ...
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getMaximum(const double E) const
Get depth of shower maximum.
Definition: JGeanz.hh:162
General purpose class for object reading from a list of file names.
bool is_scattered(const int pdf)
Test if given PDF type corresponds to scattered light.
Definition: JPDFTypes.hh:191
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:350
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:177
double getX() const
Get x position.
Definition: JVector3D.hh:94
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:67
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:139
Auxiliary data structure to list files in directory.
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
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
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
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
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
do set_variable DETECTOR_TXT $WORKDIR detector
bool is_tau(const Trk &track)
Test whether given track is a (anti-)tau.
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:41
JPosition3D transform(const JPosition3D &pos) const
Transform position.
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
int numberOfBins
number of bins for average CDF integral of optical module
Definition: JSirene.cc:64
double getZ() const
Get z position.
Definition: JVector3D.hh:115
std::vector< Trk > mc_trks
MC: list of MC truth tracks.
Definition: Evt.hh:46
Auxiliary class to set-up Trk.
Definition: JTrk_t.hh:19
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 63 of file JSirene.cc.

int numberOfBins = 200

number of bins for average CDF integral of optical module

Definition at line 64 of file JSirene.cc.

double safetyFactor = 1.7

safety factor for average CDF integral of optical module

Definition at line 65 of file JSirene.cc.