#include #include #include "TROOT.h" #include "TFile.h" #include "TH2D.h" #include "km3net-dataformat/offline/Head.hh" #include "km3net-dataformat/offline/Evt.hh" #include "JDAQ/JDAQEventIO.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JTools/JAbstractHistogram.hh" #include "JSupport/JTriggeredFileScanner.hh" #include "JSupport/JMonteCarloFileSupportkit.hh" #include "JSupport/JSupport.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * * Example program to histogram trigger efficiency. * \author mdejong, bjjung */ int main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace KM3NETDAQ; typedef JAbstractHistogram JAbstractHistogram_t; JTriggeredFileScanner<> inputFile; string outputFile; JAbstractHistogram_t X; JAbstractHistogram_t Y; int debug; try { JParser<> zap("Example program to histogram trigger efficiency."); zap['f'] = make_field(inputFile); zap['o'] = make_field(outputFile) = "efficiency.root"; zap['X'] = make_field(X, "Energy binning") = JAbstractHistogram_t(20, 1e-3, 100.0); zap['Y'] = make_field(Y, "Zenith-angle binning") = JAbstractHistogram_t(20, -1.0, 1.0); zap['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } // output TFile out(outputFile.c_str(), "recreate"); TH2D h0("h0", NULL, X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit(), Y.getNumberOfBins(), Y.getLowerLimit(), Y.getUpperLimit()); TH2D h1("h1", NULL, X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit(), Y.getNumberOfBins(), Y.getLowerLimit(), Y.getUpperLimit()); TH2D h2("h2", NULL, X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit(), Y.getNumberOfBins(), Y.getLowerLimit(), Y.getUpperLimit()); h2.Sumw2(); while (inputFile.hasNext()) { STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next(); const Evt* event = ps; vector::const_iterator hit = find_if_not(event->mc_hits.cbegin(), event->mc_hits.cend(), &is_noise); vector::const_iterator primary = find_if(event->mc_trks.cbegin(), event->mc_trks.cend(), &is_initialstate); if (hit != event->mc_hits.cend() && primary != event->mc_trks.cend()) { h1.Fill(primary->E, primary->dir.z, 1.0); } } STATUS(endl); for (JMultipleFileScanner in(inputFile); in.hasNext(); ) { STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl); const Evt* event = in.next(); vector::const_iterator hit = find_if_not(event->mc_hits.cbegin(), event->mc_hits.cend(), &is_noise); vector::const_iterator primary = find_if(event->mc_trks.cbegin(), event->mc_trks.cend(), &is_initialstate); if (hit != event->mc_hits.cend() && primary != event->mc_trks.cend()) { h0.Fill(primary->E, primary->dir.z, 1.0); } } STATUS(endl); for (Int_t i = 1; i <= h1.GetNbinsX(); ++i) { for (Int_t j = 1; j <= h1.GetNbinsY(); ++j) { const Double_t z1 = h1.GetBinContent(i,j); const Double_t z0 = h0.GetBinContent(i,j); if (z0 != 0.0) { const Double_t z3 = z1 / z0; const Double_t w3 = sqrt(z1 * (z0 - z1) / (z0*z0*z0)); h2.SetBinContent(i, j, z3); h2.SetBinError (i, j, w3); } else { ERROR("Bin " << h0.GetName() << "[" << i << "," << j << "] empty." << endl); } } } out.Write(); out.Close(); }