223 string fileDescriptor;
225 JFileRecorder <JTYPELIST<JAAnetTypes_t, JMetaTypes_t, JRootTypes_t>::typelist>
outputFile;
226 JLimit_t& numberOfEvents = inputFile.getLimit();
250 JParser<> zap(
"Main program to simulate detector response to muons and showers.");
253 zap[
'F'] =
make_field(fileDescriptor,
"file name descriptor for CDF tables");
258 zap[
's'] =
make_field(writeEMShowers,
"store generated EM showers in event");
259 zap[
'N'] =
make_field(numberOfHits,
"minimum number of hits to output event") = 1;
260 zap[
'U'] =
make_field(factor,
"scaling factor applied to light yields") = 1.0;
266 catch(
const exception &error) {
267 FATAL(error.what() << endl);
274 const JMeta meta(argc, argv);
276 const double Zbed = 0.0;
281 if (fileDescriptor !=
"") {
282 CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_MUON));
283 CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_MUON));
284 CDF.push_back(JCDF_t(fileDescriptor, DIRECT_LIGHT_FROM_DELTARAYS));
285 CDF.push_back(JCDF_t(fileDescriptor, SCATTERED_LIGHT_FROM_DELTARAYS));
287 CDG.push_back(JCDG_t(fileDescriptor, DIRECT_LIGHT_FROM_EMSHOWER));
288 CDG.push_back(JCDG_t(fileDescriptor, SCATTERED_LIGHT_FROM_EMSHOWER));
291 double maximal_road_width = 0.0;
292 double maximal_distance = 0.0;
294 for (
size_t i = 0; i != CDF.size(); ++i) {
296 DEBUG(
"Range CDF["<< CDF[i].type <<
"] " << CDF[i].function.intensity.getXmax() <<
" m" << endl);
298 maximal_road_width = max(maximal_road_width, CDF[i].function.intensity.getXmax());
301 for (
size_t i = 0; i != CDG.size(); ++i) {
303 DEBUG(
"Range CDG["<< CDG[i].type <<
"] " << CDG[i].function.intensity.getXmax() <<
" m" << endl);
306 maximal_road_width = max(maximal_road_width, CDG[i].function.intensity.getXmax());
309 maximal_distance = max(maximal_distance, CDG[i].function.intensity.getXmax());
312 NOTICE(
"Maximal road width [m] " << maximal_road_width << endl);
313 NOTICE(
"Maximal distance [m] " << maximal_distance << endl);
316 if (detectorFile ==
"" || inputFile.empty()) {
317 STATUS(
"Nothing to be done." << endl);
325 STATUS(
"Load detector... " << flush);
337 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ) {
338 if (!module->empty())
341 module = detector.erase(module);
349 STATUS(
"Setting up radiation tables... " << flush);
351 const JRadiation hydrogen ( 1.0, 1.0, 40, 0.01, 0.1, 0.1);
352 const JRadiation oxygen ( 8.0, 16.0, 40, 0.01, 0.1, 0.1);
353 const JRadiation chlorine (17.0, 35.5, 40, 0.01, 0.1, 0.1);
354 const JRadiation sodium (11.0, 23.0, 40, 0.01, 0.1, 0.1);
356 const JRadiation calcium (20.0, 40.0, 40, 0.01, 0.1, 0.1);
357 const JRadiation magnesium(12.0, 24.3, 40, 0.01, 0.1, 0.1);
358 const JRadiation potassium(19.0, 39.0, 40, 0.01, 0.1, 0.1);
359 const JRadiation sulphur (16.0, 32.0, 40, 0.01, 0.1, 0.1);
427 if (cylinder.
getZmin() < Zbed) {
431 NOTICE(
"Light generation volume: " << cylinder << endl);
439 header = inputFile.getHeader();
441 JHead buffer(header);
455 buffer.
detector.rbegin()->filename = detectorFile;
459 offset +=
Vec(cylinder.
getX(), cylinder.
getY(), 0.0);
496 copy(buffer, header);
502 NOTICE(
"Offset applied to true tracks is: " << offset << endl);
504 TH1D job(
"job", NULL, 400, 0.5, 400.5);
505 TProfile cpu(
"cpu", NULL, 16, 0.0, 8.0);
506 TProfile2D rms(
"rms", NULL, 16, 0.0, 8.0, 201, -0.5, 200.5);
507 TProfile2D rad(
"rad", NULL, 16, 0.0, 8.0, 201, -0.5, 200.5);
520 const double epsilon = 1.0e-6;
527 STATUS(
"event: " << setw(10) << in.getCounter() <<
'\r');
DEBUG(endl);
531 Evt* evt = in.next();
534 track->pos += offset;
539 event.mc_hits.clear();
546 for (vector<Trk>::const_iterator track = evt->
mc_trks.begin(); track != evt->
mc_trks.end(); ++track) {
548 if (!track->is_finalstate()) {
564 double Zmin = intersection.first;
565 double Zmax = intersection.second;
567 if (Zmax - Zmin <= parameters.Dmin_m) {
571 JVertex vertex(0.0, track->t, track->E);
573 if (track->pos.z < Zbed) {
575 if (track->dir.z > 0.0)
576 vertex.
step(
gRock, (Zbed - track->pos.z) / track->dir.z);
581 if (vertex.
getZ() < Zmin) {
582 vertex.
step(gWater, Zmin - vertex.
getZ());
585 if (vertex.
getRange() <= parameters.Dmin_m) {
593 if (subdetector.empty()) {
601 while (vertex.
getE() >= parameters.Emin_GeV && vertex.
getZ() < Zmax) {
603 const int N = radiation.size();
608 for (
int i = 0; i != N; ++i) {
609 ls += li[i] = radiation[i]->getInverseInteractionLength(vertex.
getE());
614 for (
size_t i = 0; i != ionization.size(); ++i) {
615 As += ionization[i]->getA(vertex.
getE());
618 double step = gRandom->Exp(1.0) /
ls;
621 if (vertex.
getE() < parameters.Emax_GeV) {
622 if (parameters.Dmax_m < range) {
623 range = parameters.Dmax_m;
630 rms.Fill(log10(vertex.
getE()), (Double_t) 0, ts*ts);
634 vertex.
step(As, min(step,range));
640 if (vertex.
getE() >= parameters.Emin_GeV) {
642 double y = gRandom->Uniform(
ls);
644 for (
int i = 0; i != N; ++i) {
650 Es = radiation[i]->getEnergyOfShower(vertex.
getE());
651 ts = radiation[i]->getThetaRMS(vertex.
getE(), Es);
655 rms.Fill(log10(vertex.
getE()), (Double_t) radiation[i]->getID(), ts*ts);
656 rad.Fill(log10(vertex.
getE()), (Double_t) radiation[i]->getID(), Es);
664 vertex.
applyEloss(getRandomDirection(T2), Es);
666 muon.push_back(vertex);
675 muon.push_back(vertex);
682 make_predicate(&
Trk::id, track->id));
684 if (trk != event.
mc_trks.end()) {
685 trk->len = (muon.rbegin()->getZ() < Zmax ? +1 : -1) * (muon.rbegin()->getZ() - muon.begin()->getZ());
689 for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
691 const double z0 = muon.begin()->getZ();
692 const double t0 = muon.begin()->getT();
693 const double Z =
module->getZ() - module->getX() / getTanThetaC();
695 if (Z >= muon.begin()->getZ() && Z <= muon.rbegin()->getZ()) {
697 const JVector2D pos = muon.getPosition(Z);
698 const double R = hypot(module->getX() - pos.
getX(),
699 module->getY() - pos.
getY());
701 for (
size_t i = 0; i != CDF.size(); ++i) {
703 if (R < CDF[i].integral.getXmax()) {
713 const double NPE = CDF[i].integral.getNPE(R) *
module->size() * factor * W;
714 const size_t N = getPoisson(NPE);
720 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
722 const double R = hypot(pmt->getX() - pos.
getX(),
723 pmt->getY() - pos.
getY());
724 const double theta = pi.
constrain(pmt->getTheta());
725 const double phi = pi.
constrain(fabs(pmt->getPhi()));
727 npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
732 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
734 const double R = hypot(pmt->getX() - pos.
getX(),
735 pmt->getY() - pos.
getY());
736 const double Z = pmt->getZ() - z0;
737 const double theta = pi.
constrain(pmt->getTheta());
738 const double phi = pi.
constrain(fabs(pmt->getPhi()));
740 size_t n0 = min(ns[
distance(module->begin(),pmt)], parameters.Nmax_PMT);
742 job.Fill((
double) (100 + CDF[i].type), (
double) n0);
746 const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
749 mc_hits.push_back(
JHit_t(mc_hits.size() + 1,
753 t0 + (R * getTanThetaC() + Z) / C + t1,
760 if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
761 job.Fill((
double) (300 + CDF[i].type));
765 catch(
const exception& error) {
766 job.Fill((
double) (200 + CDF[i].type));
773 for (JTrack::const_iterator vertex = muon.begin(); vertex != muon.end(); ++vertex) {
775 const double Es = vertex->
getEs();
777 if (Es >= parameters.Ecut_GeV) {
779 const double z0 = vertex->
getZ();
780 const double t0 = vertex->
getT();
783 int origin = track->id;
785 if (writeEMShowers) {
786 origin =
event.mc_trks.size() + 1;
789 int number_of_hits = 0;
791 JDetectorSubset_t::range_type range = subdetector.getRange(z0 - maximal_distance,
792 z0 + maximal_distance);
794 for (JDetector::const_iterator module = range.begin(); module != range.end(); ++module) {
796 const double R = hypot(module->getX() - vertex->
getX(),
797 module->getY() - vertex->
getY());
798 const double Z =
module->getZ() - z0 - DZ;
799 const double D = sqrt(R*R + Z*Z);
800 const double cd = Z / D;
802 for (
size_t i = 0; i != CDG.size(); ++i) {
804 if (D < CDG[i].integral.getXmax()) {
808 const double NPE = CDG[i].integral.getNPE(D, cd) * Es * module->size() * factor;
809 const size_t N = getPoisson(NPE);
815 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
817 const double R = hypot(pmt->getX() - vertex->
getX(),
818 pmt->getY() - vertex->
getY());
819 const double Z = pmt->getZ() - z0 - DZ;
820 const double D = sqrt(R*R + Z*Z);
821 const double cd = Z / D;
822 const double theta = pi.
constrain(pmt->getTheta());
823 const double phi = pi.
constrain(fabs(pmt->getPhi()));
825 npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * Es * factor);
830 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
832 const double R = hypot(pmt->getX() - vertex->
getX(),
833 pmt->getY() - vertex->
getY());
834 const double theta = pi.
constrain(pmt->getTheta());
835 const double phi = pi.
constrain(fabs(pmt->getPhi()));
837 size_t n0 = min(ns[
distance(module->begin(),pmt)], parameters.Nmax_PMT);
839 job.Fill((
double) (100 + CDG[i].type), (
double) n0);
844 const double Z = pmt->getZ() - z0 - dz;
845 const double D = sqrt(R*R + Z*Z);
846 const double cd = Z / D;
848 const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
851 mc_hits.push_back(
JHit_t(mc_hits.size() + 1,
860 number_of_hits += n1;
864 if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
865 job.Fill((
double) (300 + CDG[i].type));
869 catch(
const exception& error) {
870 job.Fill((
double) (200 + CDG[i].type));
876 if (writeEMShowers && number_of_hits != 0) {
878 event.mc_trks.push_back(
JTrk_t(origin,
881 track->pos + track->dir * vertex->
getZ(),
889 }
else if (track->len > 0.0) {
897 const double z0 = 0.0;
898 const double z1 = z0 + track->len;
899 const double t0 = track->t;
900 const double E = track->E;
906 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
910 const double R = pos.
getX();
915 R > maximal_road_width) {
919 for (
size_t i = 0; i != CDF.size(); ++i) {
931 if (R < CDF[i].integral.getXmax()) {
935 const double NPE = CDF[i].integral.getNPE(R) *
module->size() * factor * W;
936 const size_t N = getPoisson(NPE);
946 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
948 const double R = pmt->getX();
949 const double theta = pi.
constrain(pmt->getTheta());
950 const double phi = pi.
constrain(fabs(pmt->getPhi()));
952 npe.push_back(CDF[i].function.getNPE(R, theta, phi) * factor * W);
957 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
959 const double R = pmt->getX();
960 const double Z = pmt->getZ() - z0;
961 const double theta = pi.
constrain(pmt->getTheta());
962 const double phi = pi.
constrain(fabs(pmt->getPhi()));
964 size_t n0 = min(ns[
distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
966 job.Fill((
double) (120 + CDF[i].type), (
double) n0);
970 const double t1 = CDF[i].function.getTime(R, theta, phi, gRandom->Rndm());
973 mc_hits.push_back(
JHit_t(mc_hits.size() + 1,
977 t0 + (R * getTanThetaC() + Z) / C + t1,
984 if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
985 job.Fill((
double) (320 + CDF[i].type));
989 catch(
const exception& error) {
990 job.Fill((
double) (220 + CDF[i].type));
996 if (!buffer.empty()) {
1010 double E = track->E;
1015 catch(
const exception& error) {
1016 ERROR(error.what() << endl);
1019 E =
pythia(track->type, E);
1023 const double z0 = 0.0;
1024 const double t0 = track->t;
1031 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
1035 const double R = pos.
getX();
1036 const double Z = pos.
getZ() - z0 - DZ;
1037 const double D = sqrt(R*R + Z*Z);
1038 const double cd = Z / D;
1040 for (
size_t i = 0; i != CDG.size(); ++i) {
1042 if (D < CDG[i].integral.getXmax()) {
1046 const double NPE = CDG[i].integral.getNPE(D, cd) * E * module->size() * factor;
1047 const size_t N = getPoisson(NPE);
1057 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1059 const double R = pmt->getX();
1060 const double Z = pmt->getZ() - z0 - DZ;
1061 const double D = sqrt(R*R + Z*Z);
1062 const double cd = Z / D;
1063 const double theta = pi.
constrain(pmt->getTheta());
1064 const double phi = pi.
constrain(fabs(pmt->getPhi()));
1066 npe.push_back(CDG[i].function.getNPE(D, cd, theta, phi) * E * factor);
1071 for (JModule::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
1073 const double theta = pi.
constrain(pmt->getTheta());
1074 const double phi = pi.
constrain(fabs(pmt->getPhi()));
1076 size_t n0 = min(ns[
distance(buffer.cbegin(),pmt)], parameters.Nmax_PMT);
1078 job.Fill((
double) (140 + CDG[i].type), (
double) n0);
1083 const double Z = pmt->getZ() - z0 - dz;
1084 const double D = sqrt(R*R + Z*Z);
1085 const double cd = Z / D;
1087 const double t1 = CDG[i].function.getTime(D, cd, theta, phi, gRandom->Rndm());
1090 mc_hits.push_back(
JHit_t(mc_hits.size() + 1,
1094 t0 + (dz + D * getIndexOfRefraction()) / C + t1,
1101 if (std::accumulate(npe.begin(), npe.end(), 0.0) > NPE) {
1102 job.Fill((
double) (340 + CDG[i].type));
1106 catch(
const exception& error) {
1107 job.Fill((
double) (240 + CDG[i].type));
1113 if (!buffer.empty()) {
1124 if (!mc_hits.empty()) {
1126 mc_hits.
merge(parameters.Tmax_ns);
1128 event.mc_hits.resize(mc_hits.size());
1130 copy(mc_hits.begin(), mc_hits.end(), event.
mc_hits.begin());
1139 if (event.
mc_hits.size() >= numberOfHits) {