#include #include #include "TROOT.h" #include "TFile.h" #include "TH2D.h" #include "km3net-dataformat/definitions/pmt_status.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JCalibrate/JCalibrateToT.hh" #include "JTools/JRange.hh" #include "JLang/JLangToolkit.hh" #include "JSupport/JMeta.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * Auxiliary program to update the PMT status in the detector file based on the fraction of short time-over-threshold values. * * \author acreusot, mdejong */ int main(const int argc, const char * const argv[]) { using namespace std; using namespace JPP; typedef JRange JRange_t; string inputFile; string detectorFile; JRange_t range; double fraction; bool overwriteDetector; int debug; try { JParser<> zap("Auxiliary program to update the PMT status in the detector file based on the fraction of short time-over-threshold values."); zap['f'] = make_field(inputFile, "input file (output from JCalibrateToT)."); zap['a'] = make_field(detectorFile, "detector file."); zap['x'] = make_field(range, "integration range for noise.") = JRange_t(0.0, 8.0); zap['t'] = make_field(fraction, "maximal fraction of signal allowed for noise.") = 0.5; zap['A'] = make_field(overwriteDetector, "overwrite detector file."); zap['d'] = make_field(debug, "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); } detector.comment.add(JMeta(argc,argv)); if (detector.setToLatestVersion()) { NOTICE("Set detector version to " << detector.getVersion() << endl); } gErrorIgnoreLevel = kError; TFile* in = TFile::Open(inputFile.c_str(), "exist"); if (in == NULL || !in->IsOpen()) { FATAL("File: " << inputFile << " not opened." << endl); } for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) { if (!module->empty()) { TH2D* h2s = (TH2D*) in->Get(MAKE_CSTRING(module->getID() << _2SToT)); if (h2s != NULL) { for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) { double noise = 0.0; double signal = 0.0; for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) { const double x = h2s->GetYaxis()->GetBinCenter(iy); const double y = h2s->GetBinContent(ix, iy); signal += y; if (range(x)) { noise += y; } } const bool disable = noise > fraction * signal; if (debug >= debug_t || disable) { cout << "PMT " << FILL(10,'0') << module->getID() << "." << FILL(2,'0') << (ix - 1) << FILL() << " noise/signal " << FIXED(7,0) << noise << "/" << FIXED(7,0) << signal << (disable ? " *" : "") << endl; } if (overwriteDetector) { if (disable) { module->getPMT(ix - 1).getStatus().set(PMT_DISABLE); } } } } else { WARNING("No histogram for module " << module->getID() << "; skip." << endl); } } } in->Close(); if (overwriteDetector) { NOTICE("Store PMT status on file " << detectorFile << endl); store(detectorFile, detector); } }