#include #include #include "km3net-dataformat/offline/Head.hh" #include "km3net-dataformat/offline/Evt.hh" #include "km3net-dataformat/offline/Vec.hh" #include "JDAQ/JDAQEventIO.hh" #include "JGeometry3D/JCylinder3D.hh" #include "JGeometry3D/JPosition3D.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JAAnet/JHeadToolkit.hh" #include "JAAnet/JVolume.hh" #include "JSupport/JTriggeredFileScanner.hh" #include "JSupport/JEvtWeightFileScannerSet.hh" #include "JSupport/JMonteCarloFileSupportkit.hh" #include "JSupport/JSupport.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JPhysics/JConstants.hh" #include "JPhysics/JGeane.hh" #include "JReconstruction/JEvt.hh" #include "JReconstruction/JEvtToolkit.hh" #include "JROOT/JManager.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" #include "TROOT.h" #include "TFile.h" #include "TH1D.h" namespace { /** * Check if given event has a muon intersecting the given cylinder. * * \param event event * \param cylinder cylinder * \return true if event has intersecting muon; else false. */ bool hasIntersectingMuon(const Evt& event, const JGEOMETRY3D::JCylinder3D& cylinder) { using namespace std; using namespace JPP; for (vector::const_iterator i = event.mc_trks.cbegin(); i != event.mc_trks.cend(); ++i) { if (is_finalstate(*i) && is_muon(*i)) { const JCylinder3D::intersection_type& intersection = cylinder.getIntersection(getAxis(*i)); const double Lmuon = gWater.getX(i->E, MASS_MUON / getSinThetaC()); const double Leff = (intersection.first < 0.0 ? min(Lmuon, intersection.second) : min(Lmuon, intersection.second) - intersection.first); if (Leff > 0.0) { return true; } } } return false; } } /** * \file * * Example program to histogram neutrino effective mass for triggered events inside the can volume.\n * For genie/gSeaGen events the effective mass for triggered events inside the instrumented volume is also histogrammed.\n * The unit of the contents of the histogram is Mton. * \author mdejong, bstrandberg */ int main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace KM3NETDAQ; //------------------------------------------------------------------------------ // parse command line arguments, check that header can be found in the input file //------------------------------------------------------------------------------ JMultipleFileScanner_t inputFiles; string outputFile; double wall; bool logx; int numberOfBins; string detectorFile; int debug; try { JParser<> zap("Example program to histogram neutrino effective mass for triggered events."); zap['f'] = make_field(inputFiles); zap['o'] = make_field(outputFile) = "Meff.root"; zap['R'] = make_field(wall, "Addition margin around the volume in which the considered events must reside") = 0.0; zap['X'] = make_field(logx, "Use logarithm of energy"); zap['N'] = make_field(numberOfBins, "Number of bins in the energy range of the MC simulation") = 10; zap['d'] = make_field(debug) = 1; zap['a'] = make_field(detectorFile, "Detector file: if not provided, trigger fraction is not calculated") = ""; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } //------------------------------------------------------------------------- // Initialise the detector //------------------------------------------------------------------------- JDetector detector; if (detectorFile != "") { try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } } //------------------------------------------------------------------------- // Retrieve auxiliary header info //------------------------------------------------------------------------- JCylinder3D canvol; JCylinder3D instvol; vector volumes; double Xmin = numeric_limits::max(); double Xmax = numeric_limits::lowest(); JEvtWeightFileScannerSet<> scanners(inputFiles); for (JEvtWeightFileScannerSet<>::const_iterator scanner = scanners.cbegin(); scanner != scanners.cend(); ++scanner) { const JHead& header = scanner->getHeader(); const JVolume volume(header, logx); if (volume.getXmin() < Xmin) { Xmin = volume.getXmin(); } if (volume.getXmax() > Xmax) { Xmax = volume.getXmax(); } volumes.push_back(volume); // the can volume comes from the generator header and the coordinate system is consistent by construction JCylinder3D can = get(header); canvol.addMargin(wall); if (can.getVolume() > canvol.getVolume()) { canvol = can; } // as the detector is constructed from Jpp detector file, the coordinate system needs to be synced // to the one used in the event generator JPosition3D center = getPosition(get(header)); detector -= center; JCylinder3D inst( detector.begin(), detector.end() ); instvol.addMargin(wall); if (inst.getVolume() > instvol.getVolume()) { instvol = inst; } detector += center; } //------------------------------------------------------------------------------ // initialise output //------------------------------------------------------------------------------ TFile out(outputFile.c_str(), "recreate"); TH1::SetDefaultSumw2(); JManager hM(new TH1D("hM[%]", "effective mass (can) [kg]" , numberOfBins, Xmin, Xmax)); JManager hm(new TH1D("hm[%]", "effective mass (instrumented) [kg]", numberOfBins, Xmin, Xmax)); JManager hNM(new TH1D("hNM[%]", "normalization effective mass (can) [kg]" , numberOfBins, Xmin, Xmax)); JManager hNm(new TH1D("hNm[%]", "normalization effective mass (instrumented) [kg]", numberOfBins, Xmin, Xmax)); for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) { const size_t scannerIndex = distance(scanners.begin(), scanner); const JVolume& volume = volumes[scannerIndex]; NOTICE("JEffectiveMass1D: Instrumented volume dimensions (zmin, zmax, r): " << instvol.getZmin() << " " << instvol.getZmax() << " " << instvol.getRadius() << endl ); //------------------------------------------------------------------------------ // loop over generated events in the input files //------------------------------------------------------------------------------ NOTICE("Scanning generated events for file type " << scannerIndex << endl); while (scanner->hasNext()) { STATUS("event: " << setw(10) << scanner->getCounter() << '\r'); DEBUG(endl); const Evt* event = scanner->next(); const Trk& primary = get_primary(*event); const JVertex3D& vertex = getVertex(*event); const bool isValid1 = (canvol.is_inside(vertex) || hasIntersectingMuon(*event, canvol)); const bool isValid2 = (instvol.is_inside(vertex) || hasIntersectingMuon(*event, instvol)); // events are filled to hM if they are inside the can hNM[primary.type]->Fill(volume.getX(primary.E), isValid1 ? scanner->getWeight(*event) : 0.0); // events are filled to hm if they are inside the instrumented volume hNm[primary.type]->Fill(volume.getX(primary.E), isValid2 ? scanner->getWeight(*event) : 0.0); } // end loop over generated events STATUS(endl); if (scanner->getCounter() == 0) { FATAL("JEffectiveMass1D: generated events not stored in the input file, JTriggerEfficiency should be run without option -O"); } //------------------------------------------------------------------------------ // loop over triggered events in the input files //------------------------------------------------------------------------------ NOTICE("Scanning triggered events for file type " << scannerIndex << endl); JEvtWeightFileScanner > in(scanner->getFilelist()); while (in.hasNext()) { STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl); JTriggeredFileScanner<>::multi_pointer_type mp = in.next(); const Evt* event = mp; const Trk& primary = get_primary(*event); const JVertex3D& vertex = getVertex(*event); const bool isValid1 = (canvol.is_inside(vertex) || hasIntersectingMuon(*event, canvol)); const bool isValid2 = (instvol.is_inside(vertex) || hasIntersectingMuon(*event, instvol)); hM[primary.type]->Fill(volume.getX(primary.E), isValid1 ? in.getWeight(*event) : 0.0); hm[primary.type]->Fill(volume.getX(primary.E), isValid2 ? in.getWeight(*event) : 0.0); } // end loop over generated events STATUS(endl); } //------------------------------------------------------------------------------ // convert 'generated' and 'triggered' histograms to effective mass //------------------------------------------------------------------------------ vector pdgIDs = get_keys(hM); for (vector::const_iterator i = pdgIDs.begin(); i != pdgIDs.end(); ++i) { const int pdgID = *i; hM[pdgID]->Divide(hNM[pdgID]); hM[pdgID]->Scale(canvol.getVolume() * DENSITY_SEA_WATER * 1e-6); // Mton hm[pdgID]->Divide(hNm[pdgID]); hm[pdgID]->Scale(instvol.getVolume() * DENSITY_SEA_WATER * 1e-6); // Mton } out << hM << hm; out.Close(); }