#include #include #include #include #include "km3net-dataformat/offline/Head.hh" #include "km3net-dataformat/offline/Evt.hh" #include "km3net-dataformat/offline/Hit.hh" #include "JAAnet/JHead.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JMonteCarloFileSupportkit.hh" #include "JSupport/JSupport.hh" #include "JPhysics/JPDFToolkit.hh" #include "JPhysics/JPDFTransformer.hh" #include "JPhysics/JDispersion.hh" #include "JPhysics/JConstants.hh" #include "JPhysics/Antares.hh" #include "JPhysics/KM3NeT.hh" #include "JIO/JFileStreamIO.hh" #include "JGeometry3D/JPosition3D.hh" #include "JGeometry3D/JDirection3D.hh" #include "JGeometry3D/JAxis3D.hh" #include "JGeometry3D/JGeometry3DToolkit.hh" #include "JTools/JHistogram1D_t.hh" #include "JTools/JHistogramMap_t.hh" #include "JTools/JTransformableMultiHistogram.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JFunctionalMap_t.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JPMTRouter.hh" #include "JSirene/JVisibleEnergyToolkit.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { using namespace JPP; /** * Function for classifying hits generated with KM3 * * \param track muon track * \param hit hit * \return 1 for muon; 2 for shower; -1 for anything else */ inline int classify_km3(const Trk& track, const Hit& hit) { if (track.id == hit.origin) { switch(hit.type) { case HIT_TYPE_MUON_DIRECT: return 1; case HIT_TYPE_MUON_SCATTERED: return 1; case HIT_TYPE_BREMSSTRAHLUNG_DIRECT: return 2; case HIT_TYPE_BREMSSTRAHLUNG_SCATTERED: return 2; default: return -1; } } return -1; } /** * Function for classifying hits generated with KM3Sim * * \param track muon track * \param hit hit * * \return 1 for muon; 2 for shower */ inline int classify_km3sim(const Trk& track, const Hit& hit) { if (from_muon(hit)) return 1; else return 2; } } /** * \file * * Program to histogram event-by-event data of muon light for making PDFs. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace KM3NET; JMultipleFileScanner inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string outputFile; string detectorFile; int debug; try { JParser<> zap("Program to histogram event-by-event data of muon light for making PDFs."); zap['f'] = make_field(inputFile); zap['o'] = make_field(outputFile); zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['a'] = make_field(detectorFile); zap['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } JHead head(JMultipleFileScanner(inputFile).getHeader()); // Monte Carlo header const bool KM3Sim = is_km3sim(head); if (KM3Sim) { detector -= JCenter3D(detector.begin(), detector.end()); // For KM3Sim, you get it from the detector file } else { detector -= get(head); // For KM3, the MC origin is stored in the header } const JCylinder3D can = get(head); const JPMTRouter pmtRouter(detector); const double P_atm = NAMESPACE::getAmbientPressure(); const double wmax = getMaximalWavelength(); const JDispersion dispersion(P_atm); typedef JHistogram1D_t::abscissa_type abscissa_type; typedef JTransformableMultiHistogram::maplist> JMultiHistogram_t; typedef JPDFTransformer<3, abscissa_type> JFunction3DTransformer_t; JMultiHistogram_t h0; // occurrence rate of PMT (used for normalisation) JMultiHistogram_t h1; // light from muon JMultiHistogram_t h2; // light from showers const double cmin = dispersion.getKmin (wmax); h1.transformer.reset(new JFunction3DTransformer_t(35.0, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10)); h2.transformer.reset(new JFunction3DTransformer_t(35.0, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10)); // configure bins const double R[] = { 0.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 70.0, 80.0, 90.0, 100.0, 170.0, 250.0 }; for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) { const double R_m = R[i]; const double grid = 10.0 + 0.0 * R_m/100.0; // [deg] const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0)); // azimuth angle unit step size const int number_of_theta_points = max(2, (int) (180.0/(1.4 * grid))); const double theta_step = PI / (number_of_theta_points + 1); for (double theta = -0.5*theta_step; theta < PI + theta_step; theta += theta_step) { const int number_of_phi_points = max(2, (int) (PI * sin(theta) / alpha)); const double phi_step = PI / (number_of_phi_points + 1); for (double phi = -0.5*phi_step; phi < PI + phi_step; phi += phi_step) { for (JMultiHistogram_t* p : { &h0, &h1, &h2 }) { (*p)[R_m][theta][phi]; } } } } for (JMultiHistogram_t::super_iterator i1 = h1.super_begin(), i2 = h2.super_begin(); i1 != h1.super_end(); ++i1, ++i2) { for (double x : { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 7.5, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 85.0, 100.0, 125.0, 150.0, 200.0, 500.0} ) { i1.getValue()[x]; i2.getValue()[x]; } } while (inputFile.hasNext()) { NOTICE("event: " << setw(10) << inputFile.getCounter() << '\r'); STATUS(endl); const Evt* event = inputFile.next(); for (vector::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) { if (is_muon(*track)) { // track parameters const Vec Evis = getVisibleEnergy(*event, can); JVertex3D vertex = getVertex(*track); const JRotation3D R(getDirection(Evis / Evis.len())); vertex.rotate(R); JRange Z(0.0, gWater(track->E)); // range of muon for (vector::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) { try { JAxis3D axis = pmtRouter.getPMT(hit->pmt_id); axis.transform(R, vertex); const double t1 = vertex.getT() + (axis.getZ() + axis.getX()*getTanThetaC()) / getSpeedOfLight(); const double R_m = axis.getX(); const double theta = axis.getTheta(); const double phi = fabs(axis.getPhi()); const double dt = getTime(*hit) - t1; const double npe = getNPE (*hit); const int classID = (KM3Sim ? classify_km3sim(*track, *hit) : classify_km3(*track, *hit)); switch (classID) { case 1: h1.fill(R_m, theta, phi, dt, npe); // light from muon break; case 2: h2.fill(R_m, theta, phi, dt, npe / max(1.0, Evis.len())); // light from shower break; } } catch(const exception& error) { FATAL(error.what() << endl); } } for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) { JPosition3D P = module->getPosition(); P.rotate(R); P.sub(vertex); if (sqrt(P.getX()*P.getX() + P.getY()*P.getY()) <= h0.getXmax()) { for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) { JAxis3D axis = *pmt; axis.transform(R, vertex); if (Z(axis.getZ() - axis.getX()/getTanThetaC())) { h0.fill(axis.getX(), axis.getTheta(), fabs(axis.getPhi()), 0.0, 1.0); } } } } } } } NOTICE(endl); // output JFileStreamWriter out(outputFile.c_str()); for (const JMultiHistogram_t* p : { &h0, &h1, &h2 }) { out.store(*p); } out.close(); }