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 "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.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...
 

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 173 of file JSirene.cc.

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

Variable Documentation

int debug

debug level

Definition at line 61 of file JSirene.cc.

Generated by   doxygen 1.8.5