83 const int NUMBER_OF_BINS = 200;
84 const double SAFETY_FACTOR = 1.7;
90 template<
class function_type,
99 JABC(
const std::string& file_descriptor,
const JPDFType_t type)
105 const string file_name =
getFilename(file_descriptor, this->type);
107 STATUS(
"loading input from file " << file_name <<
"... " << flush);
110 function.load(file_name.c_str());
116 integral = integral_type(
function, NUMBER_OF_BINS, SAFETY_FACTOR);
122 function_type
function;
123 integral_type integral;
127 typedef JABC<JCDF4D_t, JCDF1D_t> JCDF_t;
128 typedef JABC<JCDF5D_t, JCDF2D_t> JCDG_t;
173 int main(
int argc,
char **argv)
178 string fileDescriptor;
183 double Ecut_GeV = 0.1;
184 double Emin_GeV = 1.0;
186 double Tmax_ns = 0.1;
187 int Nmax_npe = numeric_limits<int>::max();
204 JParser<> zap(
"Main program to simulate detector response to muons and showers.");
207 zap[
'F'] =
make_field(fileDescriptor,
"file name descriptor for CDF tables");
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;
221 catch(
const exception &error) {
222 FATAL(error.what() << endl);
226 gRandom->SetSeed(seed);
229 const JMeta meta(argc, argv);
231 const double Zbed = 0.0;
245 double maximal_road_width = 0.0;
246 double maximal_distance = 0.0;
248 for (
size_t i = 0; i !=
CDF.size(); ++i) {
250 DEBUG(
"Range CDF["<<
CDF[i].type <<
"] " <<
CDF[i].
function.intensity.getXmax() <<
" m" << endl);
252 maximal_road_width = max(maximal_road_width,
CDF[i].
function.intensity.getXmax());
255 for (
size_t i = 0; i != CDG.size(); ++i) {
257 DEBUG(
"Range CDG["<< CDG[i].type <<
"] " << CDG[i].
function.intensity.getXmax() <<
" m" << endl);
260 maximal_road_width = max(maximal_road_width, CDG[i].
function.intensity.getXmax());
263 maximal_distance = max(maximal_distance, CDG[i].
function.intensity.getXmax());
266 NOTICE(
"Maximal road width [m] " << maximal_road_width << endl);
267 NOTICE(
"Maximal distance [m] " << maximal_distance << endl);
270 if (detectorFile ==
"" || inputFile.empty()) {
271 STATUS(
"Nothing to be done." << endl);
279 STATUS(
"Load detector... " << flush);
294 STATUS(
"Setting up radiation tables... " << flush);
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);
329 header = inputFile.getHeader();
331 JHead buffer(header);
333 center = get<Vec>(buffer);
337 buffer.
simul.rbegin()->program =
"JSirene";
356 center +=
Vec(circle.getX(), circle.getY(), 0.0);
358 copy(buffer, header);
365 NOTICE(
"Offset applied to true tracks is: " << center << endl);
367 NOTICE(
"No offset applied to true tracks." << endl);
373 if (cylinder.getZmin() < Zbed) {
374 cylinder.setZmin(Zbed);
377 NOTICE(
"Light generation volume: " << cylinder << endl);
379 TProfile cpu(
"cpu", NULL, 14, 1.0, 8.0);
380 TH1D job(
"job", NULL, 400, 0.5, 400.5);
398 STATUS(
"event: " << setw(10) <<
in.getCounter() <<
'\r');
DEBUG(endl);
402 Evt* evt =
in.next();
406 track->pos += center;
412 event.mc_hits.clear();
433 double Zmin = intersection.first;
434 double Zmax = intersection.second;
436 if (Zmax - Zmin <= Dmin_m) {
440 JVertex vertex(0.0, track->t, track->E);
442 if (track->pos.z < Zbed) {
444 if (track->dir.z > 0.0)
445 vertex.
step(
gRock, (Zbed - track->pos.z) / track->dir.z);
450 if (vertex.
getZ() < Zmin) {
460 const JDetectorSubset_t subdetector(
detector,
getAxis(*track), maximal_road_width);
462 if (subdetector.empty()) {
470 while (vertex.
getE() >= Emin_GeV && vertex.
getZ() < Zmax) {
472 const int N = radiation.size();
477 for (
int i = 0; i !=
N; ++i) {
478 ls += li[i] = radiation[i]->getInverseInteractionLength(vertex.
getE());
481 const double ds = min(gRandom->Exp(1.0) / ls, vertex.
getRange());
485 if (vertex.
getE() >= Emin_GeV) {
488 double y = gRandom->Uniform(ls);
490 for (
int i = 0; i !=
N; ++i) {
495 Es = radiation[i]->getEnergyOfShower(vertex.
getE());
502 if (Es >= Ecut_GeV) {
503 muon.push_back(vertex);
510 if (vertex.
getE() < Emin_GeV && vertex.
getRange() > 0.0) {
514 muon.push_back(vertex);
523 if (trk != event.
mc_trks.end()) {
524 trk->len = (muon.rbegin()->getZ() < Zmax ? +1 : -1) * (muon.rbegin()->getZ() - muon.begin()->getZ());
528 for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
530 const double z0 = muon.begin()->getZ();
531 const double t0 = muon.begin()->getT();
532 const double R = module->getX();
535 if (Z >= muon.begin()->getZ() && Z <= muon.rbegin()->getZ()) {
537 for (
size_t i = 0; i !=
CDF.size(); ++i) {
539 if (R <
CDF[i].integral.getXmax()) {
549 const double NPE =
CDF[i].integral.getNPE(R) * module->size() * factor * W;
550 const int N = gRandom->Poisson(NPE);
556 for (JModule::const_iterator
pmt = module->begin();
pmt != module->end(); ++
pmt) {
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());
563 const double npe =
CDF[i].function.getNPE(R, theta, phi) * factor * W;
569 job.Fill((
double) (100 +
CDF[i].type), (
double) n0);
573 const double t1 =
CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
576 mc_hits.push_back(
JHit_t(mc_hits.size() + 1,
588 job.Fill((
double) (300 +
CDF[i].type));
592 catch(
const exception& error) {
593 job.Fill((
double) (200 +
CDF[i].type));
600 for (JTrack::const_iterator vertex = muon.begin(); vertex != muon.end(); ++vertex) {
602 const double z0 = vertex->
getZ();
603 const double t0 = vertex->
getT();
604 const double Es = vertex->
getEs();
606 int origin = track->id;
608 if (writeEMShowers) {
609 origin =
event.mc_trks.size() + 1;
612 int number_of_hits = 0;
614 JDetectorSubset_t::range_type range = subdetector.getRange(z0 - maximal_distance,
615 z0 + maximal_distance);
617 for (JDetector::const_iterator module = range.begin(); module != range.end(); ++module) {
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;
625 for (
size_t i = 0; i != CDG.size(); ++i) {
627 if (D < CDG[i].integral.getXmax()) {
631 const double NPE = CDG[i].integral.getNPE(D, cd) * Es * module->size() * factor;
632 const int N = gRandom->Poisson(NPE);
638 for (JModule::const_iterator
pmt = module->begin();
pmt != module->end(); ++
pmt) {
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;
648 const double npe = CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor;
654 job.Fill((
double) (100 + CDG[i].type), (
double) n0);
659 const double Z =
pmt->getZ() - z0 - dz;
660 const double D = sqrt(R*R + Z*Z);
661 const double cd = Z /
D;
663 const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
666 mc_hits.push_back(
JHit_t(mc_hits.size() + 1,
675 number_of_hits += n1;
680 job.Fill((
double) (300 + CDG[i].type));
684 catch(
const exception& error) {
685 job.Fill((
double) (200 + CDG[i].type));
691 if (writeEMShowers && number_of_hits != 0) {
693 event.mc_trks.push_back(
JTrk_t(origin,
696 track->pos + track->dir * vertex->
getZ(),
703 }
else if (
is_tau(*track)) {
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;
720 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
724 const double R = pos.
getX();
729 R > maximal_road_width) {
733 for (
size_t i = 0; i !=
CDF.size(); ++i) {
735 if (R <
CDF[i].integral.getXmax()) {
745 const double NPE =
CDF[i].integral.getNPE(R) * module->size() * factor * W;
746 const int N = gRandom->Poisson(NPE);
756 for (JModule::const_iterator
pmt = buffer.begin();
pmt != buffer.end(); ++
pmt) {
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());
763 const double npe =
CDF[i].function.getNPE(R, theta, phi) * factor * W;
769 job.Fill((
double) (120 +
CDF[i].type), (
double) n0);
773 const double t1 =
CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
776 mc_hits.push_back(
JHit_t(mc_hits.size() + 1,
788 job.Fill((
double) (320 +
CDF[i].type));
792 catch(
const exception& error) {
793 job.Fill((
double) (220 +
CDF[i].type));
800 if (!buffer.empty()) {
819 catch(
const exception& error) {
820 ERROR(error.what() << endl);
823 E =
pythia(track->type, E);
825 if (E < Ecut_GeV || cylinder.getDistance(
getPosition(*track)) > Dmin_m) {
829 const double z0 = 0.0;
830 const double t0 = track->t;
836 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
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;
846 if (D > maximal_distance) {
850 for (
size_t i = 0; i != CDG.size(); ++i) {
852 if (D < CDG[i].integral.getXmax()) {
856 const double NPE = CDG[i].integral.getNPE(D, cd) * E * module->size() * factor;
857 const int N = gRandom->Poisson(NPE);
867 for (JModule::const_iterator
pmt = buffer.begin();
pmt != buffer.end(); ++
pmt) {
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;
877 const double npe = CDG[i].function.getNPE(D, cd, theta, phi) * E * factor;
883 job.Fill((
double) (140 + CDG[i].type), (
double) n0);
888 const double Z =
pmt->getZ() - z0 - dz;
889 const double D = sqrt(R*R + Z*Z);
890 const double cd = Z /
D;
892 const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
895 mc_hits.push_back(
JHit_t(mc_hits.size() + 1,
907 job.Fill((
double) (340 + CDG[i].type));
911 catch(
const exception& error) {
912 job.Fill((
double) (240 + CDG[i].type));
918 if (!buffer.empty()) {
929 if (!mc_hits.empty()) {
931 mc_hits.
merge(Tmax_ns);
933 event.mc_hits.resize(mc_hits.size());
935 copy(mc_hits.begin(), mc_hits.end(),
event.mc_hits.begin());
941 cpu.Fill(log10(
get_neutrino(event).
E), (
double) timer.usec_ucpu * 1.0e-3);
944 if (event.
mc_hits.size() >= numberOfHits) {
Auxiliary class to set-up Hit.
const char *const energy_lost_in_can
Utility class to parse command line options.
Custom class for CDF table in 2 dimensions.
double getEs() const
Get shower energy.
do echo Generating $dir eval D
Multi-dimensional CDF table for arrival time of Cherenkov light.
double getT() const
Get time.
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
ROOT TTree parameter settings.
JVertex & step(const double ds)
Step.
Data structure for a composite optical module.
void applyEloss(const double Es)
Apply energy loss energy.
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
Generator for simulation.
double EfromBrems(const double E) const
Bremsstrahlung shower energy.
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.
Recording of objects on file according a format that follows from the file name extension.
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.
esac print_variable DETECTOR INPUT_FILE OUTPUT_FILE CDF for TYPE in
Utility class to parse parameter values.
static double O()
Estimated mass fractions of chemical elements in sea water.
double getE() const
Get muon energy.
Jpp environment information.
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.
static counter_type max()
Get maximum counter value.
void merge(const double Tmax_ns)
Merge hits on same PMT that are within given time window.
std::vector< JAANET::simul > simul
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.
Data structure for detector geometry and calibration.
Utility class to parse parameter values.
Fast implementation of class JRadiation.
Various implementations of functional maps.
esac $JPP_DIR examples JDetector JTransitTime o $OUTPUT_FILE n N $NPE T $TTS_NS d $DEBUG for HISTOGRAM in tts tt2 pmt
static const JPDB & getInstance()
Get particle data book.
double TotalCrossSectionBrems(const double E) const
Bremsstrahlung cross section.
Numbering scheme for PDF types.
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
Auxiliary class to extract a subset of optical modules from a detector.
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.
void push(T JHead::*pd)
Push given data member to Head.
double getZ() const
Get position.
I/O formatting auxiliaries.
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.
#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.
double getE(const double z) const
Get muon energy at given position.
int getID() const
Get identifier.
Auxiliary class for CPU timing and usage.
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.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
double EfromGNrad(const double E) const
Photo-nuclear shower energy.
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
double getKineticEnergy(const double E, const double m)
Get kinetic energy of particle with given mass.
then usage $script[distance] fi case set_variable R
General purpose messaging.
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 ...
Scanning of objects from multiple files according a format that follows from the extension of each fi...
const char * getDate()
Get ASCII formatted date.
Custom class for CDF table in 1 dimension.
direct light from delta-rays
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 ...
double getMaximum(const double E) const
Get depth of shower maximum.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
std::vector< Hit > mc_hits
MC: list of MC truth hits.
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&)).
virtual double TotalCrossSectionEErad(const double E) const
Pair production cross section.
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 class to list files.
std::string getFilename(const std::string &file_name)
Get file name part, i.e. part after last JEEP::PATHNAME_SEPARATOR if any.
Data structure for position in three dimensions.
const JLimit & getLimit() const
Get limit.
Longitudinal emission profile EM-shower.
double getLength(const double E, const double P, const double eps=1.0e-3) const
Get shower length for a given integrated probability.
double EfromEErad(const double E) const
Pair production shower energy.
virtual const char * what() const
Get error message.
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
bool is_tau(const Trk &track)
Test whether given track is a (anti-)tau.
virtual double TotalCrossSectionGNrad(const double E) const
Photo-nuclear cross section.
double getRange() const
Get range of muon.
JAANET::coord_origin coord_origin
then usage $script[input file[working directory[option]]] nWhere option can be N
Vec coord_origin() const
Get coordinate origin.
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
#define DEBUG(A)
Message macros.
JPosition3D getPosition(const Vec &v)
Get position.
int main(int argc, char *argv[])