#include #include #include #include "TROOT.h" #include "TFile.h" #include "TH1D.h" #include "TH2D.h" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JModuleRouter.hh" #include "JDetector/JPMTIdentifier.hh" #include "JCalibrate/JCalibrateK40.hh" #include "JROOT/JRootToolkit.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * * Auxiliary program to project single PMT data from 2D histogram. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; string detectorFile; string inputFile; string outputFile; JPMTIdentifier pmt; string extension; int debug; try { JParser<> zap("Auxiliary program to project single PMT data from 2D histogram."); zap['a'] = make_field(detectorFile, "detector file."); zap['f'] = make_field(inputFile, "input file."); zap['o'] = make_field(outputFile, "output file.") = "k40.root"; zap['P'] = make_field(pmt, "PMT identifier"); zap['e'] = make_field(extension, "histogram name extension") = _2R, _2S, _2F; zap['d'] = make_field(debug, "debug flag.") = 1; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } gErrorIgnoreLevel = kError; JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } const JModuleRouter router(detector); const JModule& module = router.getModule(pmt.getModuleID()); JCombinatorics combinatorics; combinatorics.configure(module.size()); combinatorics.sort(JPairwiseComparator(module)); if (debug >= debug_t) { for (size_t i = 0; i != combinatorics.getNumberOfPairs(); ++i) { const JCombinatorics::pair_type pair = combinatorics.getPair(i); cout << setw(3) << i << " -> (" << FILL(2,'0') << pair.first << "," << FILL(2,'0') << pair.second << FILL() << ")" << endl; } } TFile in(inputFile.c_str(), "read"); TH2D* h2 = (TH2D*) in.Get(MAKE_CSTRING(pmt.getModuleID() << extension)); if (h2 == NULL) { FATAL("Missing histogram for module " << pmt.getModuleID() << endl); } TH1D ha(MAKE_CSTRING(pmt.getModuleID() << '.' << FILL(2,'0') << pmt.getTDC() << ".1D"), NULL, h2->GetYaxis()->GetNbins(), h2->GetYaxis()->GetXmin(), h2->GetYaxis()->GetXmax()); TH2D hb(MAKE_CSTRING(pmt.getModuleID() << '.' << FILL(2,'0') << pmt.getTDC() << ".2D"), NULL, module.size(), -0.5, module.size() - 0.5, h2->GetYaxis()->GetNbins(), h2->GetYaxis()->GetXmin(), h2->GetYaxis()->GetXmax()); for (int i = 0; i != (int) module.size(); ++i) { if (i != pmt.getTDC()) { const Int_t ix = combinatorics.getIndex(pmt.getTDC(), i) + 1; for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) { ha.SetBinContent(iy, ha.GetBinContent(iy) + h2->GetBinContent(ix,iy)); ha.SetBinError (iy, hypot(ha.GetBinError (iy), h2->GetBinError (ix,iy))); hb.SetBinContent(i + 1, iy, h2->GetBinContent(ix,iy)); hb.SetBinError (i + 1, iy, h2->GetBinError (ix,iy)); } } } TFile out(outputFile.c_str(), "recreate"); out << ha << hb; out.Write(); out.Close(); }