Jpp
 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 "evt/Head.hh"
#include "evt/Evt.hh"
#include "JLang/JSharedPointer.hh"
#include "JLang/Jpp.hh"
#include "JLang/JPredicate.hh"
#include "JSystem/JDate.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JCDFTable.hh"
#include "JPhysics/JPDFTypes.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...
 

Detailed Description

Main program to simulate detector response to muons and showers.

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 
)

< minimal energy for generation of light from shower [GeV]

< minimal energy of muon for shower generation [GeV]

< minimal distance for positioning [m]

< maximal time between hits on same PMT to be merged

< maximal number of photo-electrons

Definition at line 167 of file JSirene.cc.

168 {
169  using namespace std;
170  using namespace JPP;
171 
172  string fileDescriptor;
173  JMultipleFileScanner<JAAnetTypes_t> inputFile;
174  JFileRecorder <JTYPELIST<JAAnetTypes_t, JMeta, JRootTypes_t>::typelist> outputFile;
175  JLimit_t& numberOfEvents = inputFile.getLimit();
176  string detectorFile;
177  double Ecut_GeV = 0.1; //!< minimal energy for generation of light from shower [GeV]
178  double Emin_GeV = 1.0; //!< minimal energy of muon for shower generation [GeV]
179  double Dmin_m = 0.1; //!< minimal distance for positioning [m]
180  double Tmax_ns = 0.1; //!< maximal time between hits on same PMT to be merged
181  int Nmax_npe = numeric_limits<int>::max(); //!< maximal number of photo-electrons
182  bool writeEMShowers;
183  size_t numberOfHits;
184  double factor;
185  UInt_t seed;
186 
187  try {
188 
189  JProperties properties;
190 
191  properties.insert(gmake_property(Ecut_GeV));
192  properties.insert(gmake_property(Emin_GeV));
193  properties.insert(gmake_property(Dmin_m));
194  properties.insert(gmake_property(Tmax_ns));
195  properties.insert(gmake_property(Nmax_npe));
196 
197  JParser<> zap("Main program to simulate detector response to muons and showers.");
198 
199  zap['@'] = make_field(properties) = JPARSER::initialised();
200  zap['F'] = make_field(fileDescriptor, "file name descriptor for CDF tables");
201  zap['f'] = make_field(inputFile) = JPARSER::initialised();
202  zap['o'] = make_field(outputFile) = "sirene.root";
203  zap['n'] = make_field(numberOfEvents) = JLimit::max();
204  zap['a'] = make_field(detectorFile) = "";
205  zap['s'] = make_field(writeEMShowers, "store generated EM showers in event");
206  zap['N'] = make_field(numberOfHits, "minimum number of hits to output event") = 1;
207  zap['U'] = make_field(factor, "scaling factor applied to light yields") = 1.0;
208  zap['S'] = make_field(seed) = 0;
209  zap['d'] = make_field(debug) = 1;
210 
211  zap(argc, argv);
212  }
213  catch(const exception &error) {
214  FATAL(error.what() << endl);
215  }
216 
217 
218  gRandom->SetSeed(seed);
219 
220 
221  const JMeta meta(argc, argv);
222 
223  const double Zbed = 0.0; // level of seabed [m]
224 
225  vector<JCDF_t> CDF;
226  vector<JCDG_t> CDG;
227 
228  CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_MUON));
229  CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_MUON));
230  CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_DELTARAYS));
231  CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_DELTARAYS));
232 
233  CDG.push_back(JCDG_t(fileDescriptor, DIRECT_LIGHT_FROM_EMSHOWER));
234  CDG.push_back(JCDG_t(fileDescriptor, SCATTERED_LIGHT_FROM_EMSHOWER));
235 
236 
237  double maximal_road_width = 0.0; // road width [m]
238  double maximal_distance = 0.0; // road width [m]
239 
240  for (size_t i = 0; i != CDF.size(); ++i) {
241 
242  DEBUG("Range CDF["<< CDF[i].type << "] " << CDF[i].function.intensity.getXmax() << " m" << endl);
243 
244  maximal_road_width = max(maximal_road_width, CDF[i].function.intensity.getXmax());
245  }
246 
247  for (size_t i = 0; i != CDG.size(); ++i) {
248 
249  DEBUG("Range CDG["<< CDG[i].type << "] " << CDG[i].function.intensity.getXmax() << " m" << endl);
250 
251  if (!is_scattered(CDF[i].type)) {
252  maximal_road_width = max(maximal_road_width, CDG[i].function.intensity.getXmax());
253  }
254 
255  maximal_distance = max(maximal_distance, CDG[i].function.intensity.getXmax());
256  }
257 
258  NOTICE("Maximal road width [m] " << maximal_road_width << endl);
259  NOTICE("Maximal distance [m] " << maximal_distance << endl);
260 
261 
262  if (detectorFile == "" || inputFile.empty()) {
263  STATUS("Nothing to be done." << endl);
264  return 0;
265  }
266 
267  JDetector detector;
268 
269  try {
270 
271  STATUS("Load detector... " << flush);
272 
273  load(detectorFile, detector);
274 
275  STATUS("OK" << endl);
276  }
277  catch(const JException& error) {
278  FATAL(error);
279  }
280 
281 
282  vector< JSharedPointer<JRadiationInterface> > radiation;
283 
284  if (true) {
285 
286  STATUS("Setting up radiation tables... " << flush);
287 
288  const JRadiation hydrogen( 1.0, 1.0, 40, 0.01, 0.1, 0.1);
289  const JRadiation oxygen ( 8.0, 16.0, 40, 0.01, 0.1, 0.1);
290  const JRadiation chlorine(17.0, 35.0, 40, 0.01, 0.1, 0.1);
291 
292  JSharedPointer<JRadiation> Hydrogen(new JRadiationFunction(hydrogen, 300, 0.2, 1.0e11));
293  JSharedPointer<JRadiation> Oxygen (new JRadiationFunction(oxygen, 300, 0.2, 1.0e11));
294  JSharedPointer<JRadiation> Chlorine(new JRadiationFunction(chlorine, 300, 0.2, 1.0e11));
295 
296  JRadiationSource::source_type EErad(&JRadiation::TotalCrossSectionEErad, &JRadiation::EfromEErad);
297  JRadiationSource::source_type Brems(&JRadiation::TotalCrossSectionBrems, &JRadiation::EfromBrems);
298  JRadiationSource::source_type GNrad(&JRadiation::TotalCrossSectionGNrad, &JRadiation::EfromGNrad);
299 
300  radiation.push_back(new JRadiationSource(Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), EErad));
301  radiation.push_back(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad));
302  radiation.push_back(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), EErad));
303 
304  radiation.push_back(new JRadiationSource(Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), Brems));
305  radiation.push_back(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems));
306  radiation.push_back(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), Brems));
307 
308  radiation.push_back(new JRadiationSource(Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), GNrad));
309  radiation.push_back(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad));
310  radiation.push_back(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), GNrad));
311 
312  STATUS("OK" << endl);
313  }
314 
315 
316  Vec center(0,0,0);
317  Head header;
318 
319  try {
320 
321  header = inputFile.getHeader();
322 
323  JHead buffer(header);
324 
325  buffer.simul.program = "JSirene";
326  buffer.simul.version = getGITVersion();
327  buffer.simul.date = getDate();
328  buffer.simul.time = getTime();
329 
330  center = get<Vec>(buffer);
331 
332  buffer.coord_origin = coord_origin(0.0, 0.0, 0.0);
333  buffer.can.zmin += center.z;
334  buffer.can.zmax += center.z;
335 
336  const JCircle2D circle(detector.begin(), detector.end());
337 
338  center += Vec(circle.getX(), circle.getY(), 0.0);
339 
340  push(buffer, &JHead::simul);
341  push(buffer, &JHead::coord_origin);
342  push(buffer, &JHead::can);
343  copy(buffer, header);
344  }
345  catch(const JException& error) {
346  FATAL(error);
347  }
348 
349  NOTICE("Offset applied to true tracks is: " << center << endl);
350 
351  JCylinder3D cylinder(detector.begin(), detector.end());
352 
353  cylinder.addMargin(maximal_distance);
354 
355  if (cylinder.getZmin() < Zbed) {
356  cylinder.setZmin(Zbed);
357  }
358 
359 
360  TProfile cpu("cpu", NULL, 14, 1.0, 8.0);
361  TH1D job("job", NULL, 400, 0.5, 400.5);
362 
363 
364  outputFile.open();
365 
366  if (!outputFile.is_open()) {
367  FATAL("Error opening file " << outputFile << endl);
368  }
369 
370  outputFile.put(meta);
371  outputFile.put(header);
372  outputFile.put(*gRandom);
373 
374 
375  JTimer timer;
376 
377  for (JMultipleFileScanner<Evt>& in = inputFile; in.hasNext(); ) {
378 
379  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
380 
381  job.Fill(1.0);
382 
383  Evt* evt = in.next(); // input
384 
385  for (vector<Trk>::iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
386  track->pos += center;
387  }
388 
389  Evt event(*evt); // output
390 
391  event.mc_hits.clear();
392 
393  JHits_t mc_hits; // temporary buffer
394 
395  timer.reset();
396  timer.start();
397 
398  for (vector<Trk>::const_iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
399 
400  if (is_muon(*track)) {
401 
402  // -----------------------------------------------
403  // muon
404  // -----------------------------------------------
405 
406  job.Fill(2.0);
407 
408  typedef JDetectorSubset<JModuleComparator3Z_t> JDetectorSubset_t;
409 
410  const pair<double, double> intersection = cylinder.getIntersection(getAxis(*track));
411 
412  double Zmin = intersection.first;
413  double Zmax = intersection.second;
414 
415  if (Zmax - Zmin <= Dmin_m) {
416  continue;
417  }
418 
419  JVertex vertex(0.0, track->t, track->E); // start of muon
420 
421  if (track->pos.z < Zbed) { // propagate muon through rock
422 
423  if (track->dir.z > 0.0)
424  vertex.step(gRock, (Zbed - track->pos.z) / track->dir.z);
425  else
426  continue;
427  }
428 
429  if (vertex.getZ() < Zmin) { // propagate muon through water
430  vertex.step(gWater, Zmin - vertex.getZ());
431  }
432 
433  if (vertex.getRange() <= Dmin_m) {
434  continue;
435  }
436 
437  job.Fill(3.0);
438 
439  const JDetectorSubset_t subdetector(detector, getAxis(*track), maximal_road_width);
440 
441  if (subdetector.empty()) {
442  continue;
443  }
444 
445  job.Fill(4.0);
446 
447  JTrack muon(vertex); // propagate muon trough detector
448 
449  while (vertex.getE() >= Emin_GeV && vertex.getZ() < Zmax) {
450 
451  const int N = radiation.size();
452 
453  double li[N]; // inverse interaction lengths
454  double ls = 1.0e-5; // minimal total inverse interaction length (100 km)^-1
455 
456  for (int i = 0; i != N; ++i) {
457  ls += li[i] = radiation[i]->getInverseInteractionLength(vertex.getE());
458  }
459 
460  const double ds = min(gRandom->Exp(1.0) / ls, vertex.getRange());
461 
462  vertex.step(ds);
463 
464  if (vertex.getE() >= Emin_GeV) {
465 
466  double Es = 0.0; // shower energy [GeV]
467  double y = gRandom->Uniform(ls);
468 
469  for (int i = 0; i != N; ++i) {
470 
471  y -= li[i];
472 
473  if (y < 0.0) {
474  Es = radiation[i]->getEnergyOfShower(vertex.getE()); // shower energy [GeV]
475  break;
476  }
477  }
478 
479  vertex.applyEloss(Es);
480 
481  if (Es >= Ecut_GeV) {
482  muon.push_back(vertex);
483  }
484  }
485  }
486 
487  // add muon end point
488 
489  if (vertex.getE() < Emin_GeV && vertex.getRange() > 0.0) {
490 
491  vertex.step(vertex.getRange());
492 
493  muon.push_back(vertex);
494  }
495 
496  // add information to output muon
497 
498  vector<Trk>::iterator trk = find_if(event.mc_trks.begin(),
499  event.mc_trks.end(),
500  make_predicate(&Trk::id, track->id));
501 
502  if (trk != event.mc_trks.end()) {
503  trk->len = (muon.rbegin()->getZ() < Zmax ? +1 : -1) * (muon.rbegin()->getZ() - muon.begin()->getZ());
504  trk->setusr(mc_usr_keys::energy_lost_in_can, muon.begin()->getE() - muon.rbegin()->getE());
505  }
506 
507  for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
508 
509  const double z0 = muon.begin()->getZ();
510  const double t0 = muon.begin()->getT();
511  const double R = module->getX();
512  const double Z = module->getZ() - R / getTanThetaC();
513 
514  if (Z >= muon.begin()->getZ() && Z <= muon.rbegin()->getZ()) {
515 
516  for (size_t i = 0; i != CDF.size(); ++i) {
517 
518  if (R < CDF[i].integral.getXmax()) {
519 
520  try {
521 
522  double W = 1.0; // mip
523 
524  if (is_deltarays(CDF[i].type)) {
525  W = getDeltaRaysFromMuon(muon.getE(Z)); // delta-rays
526  }
527 
528  const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
529  const int N = gRandom->Poisson(NPE);
530 
531  if (N != 0) {
532 
533  double ns = 0.0;
534 
535  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
536 
537  const double R = pmt->getX();
538  const double Z = pmt->getZ() - z0;
539  const double theta = pmt->getTheta();
540  const double phi = fabs(pmt->getPhi());
541 
542  const double npe = CDF[i].function.getNPE(R, theta, phi) * factor * W;
543 
544  ns += npe;
545 
546  int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), Nmax_npe);
547 
548  job.Fill((double) (100 + CDF[i].type), (double) n0);
549 
550  while (n0 != 0) {
551 
552  const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
553  const int n1 = getNumberOfPhotoElectrons(n0);
554 
555  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
556  pmt->getID(),
557  getHitType(CDF[i].type),
558  track->id,
559  t0 + (R * getTanThetaC() + Z) / C + t1,
560  n1));
561 
562  n0 -= n1;
563  }
564  }
565 
566  if (ns > NPE) {
567  job.Fill((double) (300 + CDF[i].type));
568  }
569  }
570  }
571  catch(const exception& error) {
572  job.Fill((double) (200 + CDF[i].type));
573  }
574  }
575  }
576  }
577  }
578 
579  for (JTrack::const_iterator vertex = muon.begin(); vertex != muon.end(); ++vertex) {
580 
581  const double z0 = vertex->getZ();
582  const double t0 = vertex->getT();
583  const double Es = vertex->getEs();
584 
585  int origin = track->id;
586 
587  if (writeEMShowers) {
588  origin = event.mc_trks.size() + 1;
589  }
590 
591  int number_of_hits = 0;
592 
593  JDetectorSubset_t::range_type range = subdetector.getRange(z0 - maximal_distance,
594  z0 + maximal_distance);
595 
596  for (JDetector::const_iterator module = range.begin(); module != range.end(); ++module) {
597 
598  const double dz = geanz.getMaximum(Es);
599  const double R = module->getX();
600  const double Z = module->getZ() - z0 - dz;
601  const double D = sqrt(R*R + Z*Z);
602  const double cd = Z / D;
603 
604  for (size_t i = 0; i != CDG.size(); ++i) {
605 
606  if (D < CDG[i].integral.getXmax()) {
607 
608  try {
609 
610  const double NPE = CDG[i].integral.getNPE(D, cd) * Es * module->size() * factor;
611  const int N = gRandom->Poisson(NPE);
612 
613  if (N != 0) {
614 
615  double ns = 0.0;
616 
617  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
618 
619  const double dz = geanz.getMaximum(Es);
620  const double R = pmt->getX();
621  const double Z = pmt->getZ() - z0 - dz;
622  const double theta = pmt->getTheta();
623  const double phi = fabs(pmt->getPhi());
624  const double D = sqrt(R*R + Z*Z);
625  const double cd = Z / D;
626 
627  const double npe = CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor;
628 
629  ns += npe;
630 
631  int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), Nmax_npe);
632 
633  job.Fill((double) (100 + CDG[i].type), (double) n0);
634 
635  while (n0 != 0) {
636 
637  const double dz = geanz.getLength(Es, gRandom->Rndm());
638  const double Z = pmt->getZ() - z0 - dz;
639  const double D = sqrt(R*R + Z*Z);
640  const double cd = Z / D;
641 
642  const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
643  const int n1 = getNumberOfPhotoElectrons(n0);
644 
645  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
646  pmt->getID(),
647  getHitType(CDG[i].type),
648  origin,
649  t0 + (dz + D * getIndexOfRefraction()) / C + t1,
650  n1));
651 
652  n0 -= n1;
653 
654  number_of_hits += n1;
655  }
656  }
657 
658  if (ns > NPE) {
659  job.Fill((double) (300 + CDG[i].type));
660  }
661  }
662  }
663  catch(const exception& error) {
664  job.Fill((double) (200 + CDG[i].type));
665  }
666  }
667  }
668  }
669 
670  if (writeEMShowers && number_of_hits != 0) {
671 
672  event.mc_trks.push_back(JTrk_t(origin,
673  TRACK_TYPE_PHOTON,
674  track->id,
675  track->pos + track->dir * vertex->getZ(),
676  track->dir,
677  vertex->getT(),
678  Es));
679  }
680  }
681 
682  } else if (is_tau(*track)) {
683 
684  // -----------------------------------------------
685  // tau
686  // -----------------------------------------------
687 
688  job.Fill(6.0);
689 
690  const double z0 = 0.0;
691  const double z1 = z0 + track->len;
692  const double t0 = track->t;
693  const double E = track->E;
694 
695  const JTransformation3D transformation = getTransformation(*track);
696 
697  JModule buffer;
698 
699  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
700 
701  const JPosition3D pos = transformation.transform(module->getPosition());
702 
703  const double R = pos.getX();
704  const double Z = pos.getZ() - R / getTanThetaC();
705 
706  if (Z < z0 ||
707  Z > z1 ||
708  R > maximal_road_width) {
709  continue;
710  }
711 
712  for (size_t i = 0; i != CDF.size(); ++i) {
713 
714  if (R < CDF[i].integral.getXmax()) {
715 
716  try {
717 
718  double W = 1.0; // mip
719 
720  if (is_deltarays(CDF[i].type)) {
721  W = getDeltaRaysFromTau(E); // delta-rays
722  }
723 
724  const double NPE = CDF[i].integral.getNPE(R) * module->size() * factor * W;
725  const int N = gRandom->Poisson(NPE);
726 
727  if (N != 0) {
728 
729  buffer = *module;
730 
731  buffer.transform(transformation);
732 
733  double ns = 0.0;
734 
735  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
736 
737  const double R = pmt->getX();
738  const double Z = pmt->getZ() - z0;
739  const double theta = pmt->getTheta();
740  const double phi = fabs(pmt->getPhi());
741 
742  const double npe = CDF[i].function.getNPE(R, theta, phi) * factor * W;
743 
744  ns += npe;
745 
746  int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), Nmax_npe);
747 
748  job.Fill((double) (120 + CDF[i].type), (double) n0);
749 
750  while (n0 != 0) {
751 
752  const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
753  const int n1 = getNumberOfPhotoElectrons(n0);
754 
755  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
756  pmt->getID(),
757  getHitType(CDF[i].type),
758  track->id,
759  t0 + (R * getTanThetaC() + Z) / C + t1,
760  n1));
761 
762  n0 -= n1;
763  }
764  }
765 
766  if (ns > NPE) {
767  job.Fill((double) (320 + CDF[i].type));
768  }
769  }
770  }
771  catch(const exception& error) {
772  job.Fill((double) (220 + CDF[i].type));
773  }
774 
775  }
776  }
777  }
778 
779  if (!buffer.empty()) {
780  job.Fill(7.0);
781  }
782 
783  } else if (!is_neutrino(*track)) {
784 
785  if (JPDB::getInstance().hasPDG(track->type)) {
786 
787  // -----------------------------------------------
788  // electron or hadron
789  // -----------------------------------------------
790 
791  job.Fill(8.0);
792 
793  double E = track->E;
794 
795  try {
796  E = getKineticEnergy(E, JPDB::getInstance().getPDG(track->type).mass);
797  }
798  catch(const exception& error) {
799  ERROR(error.what() << endl);
800  }
801 
802  E = pythia(track->type, E);
803 
804  if (E < Ecut_GeV || cylinder.getDistance(getPosition(*track)) > Dmin_m) {
805  continue;
806  }
807 
808  const double z0 = 0.0;
809  const double t0 = track->t;
810 
811  const JTransformation3D transformation = getTransformation(*track);
812 
813  JModule buffer;
814 
815  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
816 
817  const JPosition3D pos = transformation.transform(module->getPosition());
818 
819  const double dz = geanz.getMaximum(E);
820  const double R = pos.getX();
821  const double Z = pos.getZ() - z0 - dz;
822  const double D = sqrt(R*R + Z*Z);
823  const double cd = Z / D;
824 
825  if (D > maximal_distance) {
826  continue;
827  }
828 
829  for (size_t i = 0; i != CDG.size(); ++i) {
830 
831  if (D < CDG[i].integral.getXmax()) {
832 
833  try {
834 
835  const double NPE = CDG[i].integral.getNPE(D, cd) * E * module->size() * factor;
836  const int N = gRandom->Poisson(NPE);
837 
838  if (N != 0) {
839 
840  buffer = *module;
841 
842  buffer.transform(transformation);
843 
844  double ns = 0.0;
845 
846  for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
847 
848  const double dz = geanz.getMaximum(E);
849  const double R = pmt->getX();
850  const double Z = pmt->getZ() - z0 - dz;
851  const double theta = pmt->getTheta();
852  const double phi = fabs(pmt->getPhi());
853  const double D = sqrt(R*R + Z*Z);
854  const double cd = Z / D;
855 
856  const double npe = CDG[i].function.getNPE(D, cd, theta, phi) * E * factor;
857 
858  ns += npe;
859 
860  int n0 = min(getNumberOfPhotoElectrons(NPE, N, npe), Nmax_npe);
861 
862  job.Fill((double) (140 + CDG[i].type), (double) n0);
863 
864  while (n0 != 0) {
865 
866  const double dz = geanz.getLength(E, gRandom->Rndm());
867  const double Z = pmt->getZ() - z0 - dz;
868  const double D = sqrt(R*R + Z*Z);
869  const double cd = Z / D;
870 
871  const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
872  const int n1 = getNumberOfPhotoElectrons(n0);
873 
874  mc_hits.push_back(JHit_t(mc_hits.size() + 1,
875  pmt->getID(),
876  getHitType(CDG[i].type, true),
877  track->id,
878  t0 + (dz + D * getIndexOfRefraction()) / C + t1,
879  n1));
880 
881  n0 -= n1;
882  }
883  }
884 
885  if (ns > NPE) {
886  job.Fill((double) (340 + CDG[i].type));
887  }
888  }
889  }
890  catch(const exception& error) {
891  job.Fill((double) (240 + CDG[i].type));
892  }
893  }
894  }
895  }
896 
897  if (!buffer.empty()) {
898  job.Fill(9.0);
899  }
900 
901  } else {
902 
903  job.Fill(21.0);
904  }
905  }
906  }
907 
908  if (!mc_hits.empty()) {
909 
910  mc_hits.merge(Tmax_ns);
911 
912  event.mc_hits.resize(mc_hits.size());
913 
914  copy(mc_hits.begin(), mc_hits.end(), event.mc_hits.begin());
915  }
916 
917  timer.stop();
918 
919  if (has_neutrino(event)) {
920  cpu.Fill(log10(get_neutrino(event).E), (double) timer.usec_ucpu * 1.0e-3);
921  }
922 
923  if (event.mc_hits.size() >= numberOfHits) {
924 
925  outputFile.put(event);
926 
927  job.Fill(10.0);
928  }
929  }
930  STATUS(endl);
931 
932  outputFile.put(cpu);
933  outputFile.put(job);
934  outputFile.put(*gRandom);
935  outputFile.close();
936 }
JTOOLS::H
static const double H
Planck constant [eV s].
Definition: JConstants.hh:25
JAANET::JHit_t
Auxiliary class to set-up Hit.
Definition: JHit_t.hh:22
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
JSUPPORT::JFileRecorder
Object writing to file.
Definition: JFileRecorder.hh:41
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1410
JLANG::JException
General exception.
Definition: JException.hh:40
JPHYSICS::gWater
static const JGeane gWater(2.67e-1 *JTOOLS::DENSITY_SEA_WATER, 3.4e-4 *JTOOLS::DENSITY_SEA_WATER)
Function object for Energy loss of muon in sea water.
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
JDETECTOR::JModule
Data structure for a composite optical module.
Definition: JModule.hh:47
JTOOLS::getIndexOfRefraction
double getIndexOfRefraction()
Get average index of refraction of water.
Definition: JConstants.hh:111
JAANET::getTransformation
JTransformation3D getTransformation(const Trk &track)
Get transformation.
Definition: JAAnetToolkit.hh:280
JGEOMETRY2D::JCircle2D
Data structure for circle in two dimensions.
Definition: JCircle2D.hh:29
gmake_property
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
Definition: JProperties.hh:1307
JTOOLS::DENSITY_SEA_WATER
static const double DENSITY_SEA_WATER
Fixed environment values.
Definition: JConstants.hh:34
STATUS
#define STATUS(A)
Definition: JMessage.hh:61
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:77
JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition: JAAnetToolkit.hh:367
JSIRENE::getHitType
JHitType_t getHitType(const JPDFType_t pdf, const bool shower=false)
Get hit type corresponding to given PDF type.
Definition: JSireneToolkit.hh:112
JPHYSICS::getDeltaRaysFromTau
double getDeltaRaysFromTau(const double E)
Equivalent EM-shower energy due to delta-rays per unit tau track length.
Definition: JPDFToolkit.hh:109
JEEP::JProperties
Utility class to parse parameter values.
Definition: JProperties.hh:484
event
Structure to store the ToT mean and standard deviation of the hits produced by a nanobeacon in a sour...
Definition: JVoltageOptimizer_utils.hh:26
JAANET::push
void push(JHead &header, T JHead::*p)
Push data of JHead for subsequent copy to Head.
Definition: JHeadToolkit.hh:102
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e.
Definition: JParser.hh:64
JAANET::JHits_t::merge
void merge(const double Tmax_ns)
Merge hits on same PMT that are within given time window.
Definition: JHit_t.hh:115
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:87
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JPHYSICS::gRock
static const JGeane gRock(2.67e-1 *0.9 *JTOOLS::DENSITY_ROCK, 3.4e-4 *1.2 *JTOOLS::DENSITY_ROCK)
Function object for Energy loss of muon in rock.
JPHYSICS::DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition: JPDFTypes.hh:29
JAANET::coord_origin
Coordinate origin.
Definition: JHead.hh:432
JPHYSICS::JRadiationFunction
Fast implementation of class JRadiation.
Definition: JRadiationSource.hh:35
JPHYSICS::JRadiationSource
Source of radiation.
Definition: JRadiationSource.hh:231
JAANET::is_neutrino
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
Definition: JAAnetToolkit.hh:349
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
JGEOMETRY3D::JCylinder3D
Cylinder object.
Definition: JCylinder3D.hh:37
JAANET::detector
Detector file.
Definition: JHead.hh:126
JLANG::JSharedPointer
The template JSharedPointer class can be used to share a pointer to an object.
Definition: JSharedPointer.hh:28
JPHYSICS::JRadiation
Auxiliary class for the calculation of the muon radiative cross sections.
Definition: JRadiation.hh:39
JLANG::getInstance
T & getInstance(const T &object)
Get static instance from temporary object.
Definition: JObject.hh:73
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
JPHYSICS::geanz
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
JAANET::getAxis
JAxis3D getAxis(const Trk &track)
Get axis.
Definition: JAAnetToolkit.hh:244
JEEP::JTimer
Auxiliary class for CPU timing and usage.
Definition: JTimer.hh:32
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition: JPDFTypes.hh:36
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:62
ERROR
#define ERROR(A)
Definition: JMessage.hh:64
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition: JPDFTypes.hh:40
JLANG::getGITVersion
const char * getGITVersion()
Get GIT version.
Definition: Jpp.hh:53
JSIRENE::JTrack
Muon trajectory.
Definition: JSireneToolkit.hh:303
JAANET::JHits_t
Auxiliary data structure for list of hits with hit merging capability.
Definition: JHit_t.hh:104
JGEOMETRY3D::JCylinder3D::addMargin
void addMargin(const double D)
Add (safety) margin.
Definition: JCylinder3D.hh:163
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:444
debug
int debug
debug level
Definition: JSirene.cc:59
JSIRENE::getKineticEnergy
double getKineticEnergy(const double E, const double m)
Get kinetic energy of particle with given mass.
Definition: JSireneToolkit.hh:348
JDETECTOR::JDetectorSubset
Detector subset with binary search functionality.
Definition: JDetectorSubset.hh:118
JDETECTOR::JDetectorSubset_t
Detector subset without binary search functionality.
Definition: JDetectorSubset.hh:37
JAANET::JHead
Monte Carlo run header.
Definition: JHead.hh:727
JSIRENE::JVertex
Vertex of energy loss of muon.
Definition: JSireneToolkit.hh:146
FATAL
#define FATAL(A)
Definition: JMessage.hh:65
JSYSTEM::getDate
const char * getDate()
Get ASCII formatted date.
Definition: JDate.hh:22
JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition: JPDFTypes.hh:35
JSUPPORT::JMultipleFileScanner::hasNext
virtual bool hasNext()
Check availability of next element.
Definition: JMultipleFileScanner.hh:289
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:54
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:98
JPHYSICS::is_scattered
bool is_scattered(const int pdf)
Test if given PDF type corresponds to scattered light.
Definition: JPDFTypes.hh:185
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:213
JPHYSICS::is_deltarays
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:171
JGEOMETRY3D::JVector3D::getX
double getX() const
Get x position.
Definition: JVector3D.hh:92
JPHYSICS::getDeltaRaysFromMuon
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:81
JAANET::copy
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:40
JSYSTEM::ls
Auxiliary class to list files.
Definition: JSystemToolkit.hh:296
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:72
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:375
JTOOLS::C
static const double C
Speed of light in vacuum [m/ns].
Definition: JConstants.hh:22
JTOOLS::getTanThetaC
double getTanThetaC()
Get average tangent of Cherenkov angle of water.
Definition: JConstants.hh:122
JGEOMETRY3D::JTransformation3D::transform
JPosition3D transform(const JPosition3D &pos) const
Transform position.
Definition: JTransformation3D.hh:49
JGEOMETRY3D::JVector3D::getZ
double getZ() const
Get z position.
Definition: JVector3D.hh:113
JAANET::JTrk_t
Auxiliary class to set-up Trk.
Definition: JTrk_t.hh:19
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
JAANET::getPosition
JPosition3D getPosition(const Vec &v)
Get position.
Definition: JAAnetToolkit.hh:197

Variable Documentation

int debug

debug level

Definition at line 59 of file JSirene.cc.

Generated by   doxygen 1.8.5