180 string fileDescriptor;
185 double Ecut_GeV = 0.1;
186 double Emin_GeV = 1.0;
188 double Tmax_ns = 0.1;
189 int Nmax_npe = numeric_limits<int>::max();
208 JParser<> zap(
"Main program to simulate detector response to muons and showers.");
211 zap[
'F'] =
make_field(fileDescriptor,
"file name descriptor for CDF tables");
214 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
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;
225 catch(
const exception &error) {
226 FATAL(error.what() << endl);
230 gRandom->SetSeed(seed);
233 const JMeta meta(argc, argv);
235 const double Zbed = 0.0;
249 double maximal_road_width = 0.0;
250 double maximal_distance = 0.0;
252 for (
size_t i = 0; i !=
CDF.size(); ++i) {
254 DEBUG(
"Range CDF["<<
CDF[i].type <<
"] " <<
CDF[i].
function.intensity.getXmax() <<
" m" << endl);
256 maximal_road_width = max(maximal_road_width,
CDF[i].
function.intensity.getXmax());
259 for (
size_t i = 0; i != CDG.size(); ++i) {
261 DEBUG(
"Range CDG["<< CDG[i].type <<
"] " << CDG[i].
function.intensity.getXmax() <<
" m" << endl);
264 maximal_road_width = max(maximal_road_width, CDG[i].
function.intensity.getXmax());
267 maximal_distance = max(maximal_distance, CDG[i].
function.intensity.getXmax());
270 NOTICE(
"Maximal road width [m] " << maximal_road_width << endl);
271 NOTICE(
"Maximal distance [m] " << maximal_distance << endl);
274 if (detectorFile ==
"" || inputFile.empty()) {
275 STATUS(
"Nothing to be done." << endl);
283 STATUS(
"Load detector... " << flush);
295 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ) {
296 if (!module->empty())
306 STATUS(
"Setting up radiation tables... " << flush);
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);
383 header = inputFile.getHeader();
385 JHead buffer(header);
387 center = get<Vec>(buffer);
393 buffer.simul.rbegin()->date =
getDate();
394 buffer.simul.rbegin()->time =
getTime();
399 buffer.detector.rbegin()->filename = detectorFile;
401 buffer.push(&JHead::simul);
407 buffer.can.zmin += center.z;
408 buffer.can.zmax += center.z;
410 buffer.push(&JHead::coord_origin);
411 buffer.push(&JHead::can);
416 center +=
Vec(circle.getX(), circle.getY(), 0.0);
418 copy(buffer, header);
425 NOTICE(
"Offset applied to true tracks is: " << center << endl);
427 NOTICE(
"No offset applied to true tracks." << endl);
433 if (cylinder.getZmin() < Zbed) {
434 cylinder.setZmin(Zbed);
437 NOTICE(
"Light generation volume: " << cylinder << endl);
439 TProfile cpu(
"cpu", NULL, 14, 1.0, 8.0);
440 TH1D job(
"job", NULL, 400, 0.5, 400.5);
458 STATUS(
"event: " << setw(10) <<
in.getCounter() <<
'\r');
DEBUG(endl);
462 Evt* evt =
in.next();
466 track->pos += center;
472 event.mc_hits.clear();
481 if (!track->is_finalstate())
continue;
495 double Zmin = intersection.first;
496 double Zmax = intersection.second;
498 if (Zmax - Zmin <= Dmin_m) {
502 JVertex vertex(0.0, track->t, track->E);
504 if (track->pos.z < Zbed) {
506 if (track->dir.z > 0.0)
507 vertex.step(
gRock, (Zbed - track->pos.z) / track->dir.z);
512 if (vertex.getZ() < Zmin) {
513 vertex.step(
gWater, Zmin - vertex.getZ());
516 if (vertex.getRange() <= Dmin_m) {
522 const JDetectorSubset_t subdetector(
detector,
getAxis(*track), maximal_road_width);
524 if (subdetector.empty()) {
532 while (vertex.getE() >= Emin_GeV && vertex.getZ() < Zmax) {
534 const int N = radiation.size();
539 for (
int i = 0; i !=
N; ++i) {
540 ls += li[i] = radiation[i]->getInverseInteractionLength(vertex.getE());
545 for (
size_t i = 0; i != ionization.size(); ++i) {
546 As += ionization[i]->getA(vertex.getE());
549 const double ds = min(gRandom->Exp(1.0) / ls, vertex.getRange(As));
553 if (vertex.getE() >= Emin_GeV) {
556 double y = gRandom->Uniform(ls);
558 for (
int i = 0; i !=
N; ++i) {
563 Es = radiation[i]->getEnergyOfShower(vertex.getE());
568 vertex.applyEloss(Es);
570 if (Es >= Ecut_GeV) {
571 muon.push_back(vertex);
578 if (vertex.getE() < Emin_GeV && vertex.getRange() > 0.0) {
580 vertex.step(vertex.getRange());
582 muon.push_back(vertex);
591 if (trk != event.mc_trks.end()) {
592 trk->len = (muon.rbegin()->getZ() < Zmax ? +1 : -1) * (muon.rbegin()->getZ() - muon.begin()->getZ());
596 for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
598 const double z0 = muon.begin()->getZ();
599 const double t0 = muon.begin()->getT();
600 const double R = module->getX();
603 if (Z >= muon.begin()->getZ() && Z <= muon.rbegin()->getZ()) {
605 for (
size_t i = 0; i !=
CDF.size(); ++i) {
607 if (R <
CDF[i].integral.getXmax()) {
617 const double NPE =
CDF[i].integral.getNPE(R) * module->size() * factor * W;
618 const int N = gRandom->Poisson(NPE);
624 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
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());
631 const double npe =
CDF[i].function.getNPE(R, theta, phi) * factor * W;
637 job.Fill((
double) (100 +
CDF[i].type), (
double) n0);
641 const double t1 =
CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
644 mc_hits.push_back(
JHit_t(mc_hits.size() + 1,
656 job.Fill((
double) (300 +
CDF[i].type));
660 catch(
const exception& error) {
661 job.Fill((
double) (200 +
CDF[i].type));
668 for (JTrack::const_iterator vertex = muon.begin(); vertex != muon.end(); ++vertex) {
670 const double z0 = vertex->getZ();
671 const double t0 = vertex->getT();
672 const double Es = vertex->getEs();
674 int origin = track->id;
676 if (writeEMShowers) {
677 origin =
event.mc_trks.size() + 1;
680 int number_of_hits = 0;
682 JDetectorSubset_t::range_type
range = subdetector.getRange(z0 - maximal_distance,
683 z0 + maximal_distance);
685 for (JDetector::const_iterator module = range.begin(); module != range.end(); ++module) {
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;
693 for (
size_t i = 0; i != CDG.size(); ++i) {
695 if (D < CDG[i].integral.getXmax()) {
699 const double NPE = CDG[i].integral.getNPE(D, cd) * Es * module->size() * factor;
700 const int N = gRandom->Poisson(NPE);
706 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
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;
716 const double npe = CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor;
722 job.Fill((
double) (100 + CDG[i].type), (
double) n0);
727 const double Z = pmt->getZ() - z0 - dz;
728 const double D = sqrt(R*R + Z*Z);
729 const double cd = Z /
D;
731 const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
734 mc_hits.push_back(
JHit_t(mc_hits.size() + 1,
743 number_of_hits += n1;
748 job.Fill((
double) (300 + CDG[i].type));
752 catch(
const exception& error) {
753 job.Fill((
double) (200 + CDG[i].type));
759 if (writeEMShowers && number_of_hits != 0) {
761 event.mc_trks.push_back(
JTrk_t(origin,
764 track->pos + track->dir * vertex->getZ(),
771 }
else if (track->len>0) {
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;
788 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
792 const double R = pos.
getX();
797 R > maximal_road_width) {
801 for (
size_t i = 0; i !=
CDF.size(); ++i) {
810 if (R <
CDF[i].integral.getXmax()) {
814 const double NPE =
CDF[i].integral.getNPE(R) * module->size() * factor * W;
815 const int N = gRandom->Poisson(NPE);
825 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
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());
832 const double npe =
CDF[i].function.getNPE(R, theta, phi) * factor * W;
838 job.Fill((
double) (120 +
CDF[i].type), (
double) n0);
842 const double t1 =
CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
845 mc_hits.push_back(
JHit_t(mc_hits.size() + 1,
857 job.Fill((
double) (320 +
CDF[i].type));
861 catch(
const exception& error) {
862 job.Fill((
double) (220 +
CDF[i].type));
869 if (!buffer.empty()) {
888 catch(
const exception& error) {
889 ERROR(error.what() << endl);
892 E =
pythia(track->type, E);
894 if (E < Ecut_GeV || cylinder.getDistance(
getPosition(*track)) > Dmin_m) {
898 const double z0 = 0.0;
899 const double t0 = track->t;
905 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
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;
915 if (D > maximal_distance) {
919 for (
size_t i = 0; i != CDG.size(); ++i) {
921 if (D < CDG[i].integral.getXmax()) {
925 const double NPE = CDG[i].integral.getNPE(D, cd) * E * module->size() * factor;
926 const int N = gRandom->Poisson(NPE);
936 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
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;
946 const double npe = CDG[i].function.getNPE(D, cd, theta, phi) * E * factor;
952 job.Fill((
double) (140 + CDG[i].type), (
double) n0);
957 const double Z = pmt->getZ() - z0 - dz;
958 const double D = sqrt(R*R + Z*Z);
959 const double cd = Z /
D;
961 const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
964 mc_hits.push_back(
JHit_t(mc_hits.size() + 1,
976 job.Fill((
double) (340 + CDG[i].type));
980 catch(
const exception& error) {
981 job.Fill((
double) (240 + CDG[i].type));
987 if (!buffer.empty()) {
997 if (!mc_hits.empty()) {
999 mc_hits.
merge(Tmax_ns);
1001 event.mc_hits.resize(mc_hits.size());
1003 copy(mc_hits.begin(), mc_hits.end(),
event.mc_hits.begin());
1009 cpu.Fill(log10(
get_neutrino(event).E), (
double) timer.usec_ucpu * 1.0e-3);
1012 if (event.mc_hits.size() >= numberOfHits) {
Auxiliary class to set-up Hit.
const char *const energy_lost_in_can
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
Utility class to parse command line options.
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Data structure for a composite optical module.
JTransformation3D getTransformation(const Trk &track)
Get transformation.
Data structure for circle in two dimensions.
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
double safetyFactor
safety factor for average CDF integral of optical module
Generator for simulation.
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
scattered light from EM shower
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
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.
Utility class to parse parameter values.
static const double H
Planck constant [eV s].
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
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.
void merge(const double Tmax_ns)
Merge hits on same PMT that are within given time window.
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.
double getTime(const Hit &hit)
Get true time of hit.
static const char *const APPLICATION_JSIRENE
detector simulation
Fast implementation of class JRadiation.
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
static const double C
Physics constants.
scattered light from muon
Auxiliary class for defining the range of iterations of objects.
The Vec class is a straightforward 3-d vector, which also works in pyroot.
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.
T & getInstance(const T &object)
Get static instance from temporary object.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
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.
JPosition3D getPosition(const Vec &pos)
Get position.
scattered light from delta-rays
direct light from EM shower
Auxiliary data structure for list of hits with hit merging capability.
void addMargin(const double D)
Add (safety) margin.
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
Detector subset with binary search functionality.
Detector subset without binary search functionality.
Vertex of energy loss of muon.
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
direct light from delta-rays
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.
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.
void transform(const JRotation3D &R, const JVector3D &pos)
Transformation of geometry (see method JGEOMETRY3D::JPosition3D::transform(const JRotation3D&, const JVector3D&)).
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.
double getX() const
Get x position.
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
void copy(const Head &from, JHead &to)
Copy header from from to to.
Auxiliary data structure to list files in directory.
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Data structure for position in three dimensions.
const JLimit & getLimit() const
Get limit.
double getLength(const double E, const double P, const double eps=1.0e-3) const
Get shower length for a given integrated probability.
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
then usage $script[input file[working directory[option]]] nWhere option can be N
do echo Generating $dir eval D
int numberOfBins
number of bins for average CDF integral of optical module
double getZ() const
Get z position.
std::vector< Trk > mc_trks
MC: list of MC truth tracks.
Auxiliary class to set-up Trk.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
then usage $script[input file[working directory[option]]] nWhere option can be E