#include #include #include #include "TROOT.h" #include "TFile.h" #include "TH2D.h" #include "km3net-dataformat/online/JDAQHeader.hh" #include "JDB/JDB.hh" #include "JDB/JSelector.hh" #include "JDB/JSelectorSupportkit.hh" #include "JDB/JDBToolkit.hh" #include "JDB/JPersons.hh" #include "JDB/JPMTThreshold.hh" #include "JDB/JDetectorIntegration.hh" #include "JCalibrate/JCalibrateToT.hh" #include "JDetector/JDetectorCalibration.hh" #include "JTools/JRange.hh" #include "JLang/JLangToolkit.hh" #include "JLang/JComparator.hh" #include "JROOT/JRootFileReader.hh" #include "JSon/JSon.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { using KM3NETDAQ::JDAQHeader; using JLANG::JException; /** * Get DAQ header from file. * * \param file_name file name * \return DAQ header */ inline JDAQHeader getDAQHeader(const char* const file_name) { using namespace JPP; JRootFileReader in(file_name); if (in.hasNext()) { const JDAQHeader* p = in.next(); if (!in.hasNext()) { return *p; } } THROW(JException, "Error reading DAQ header from file " << file_name); } } /** * \file * Auxiliary program to compute the PMT thresholds according to the small time-over-threshold fraction. * * \author acreusot, mdejong */ int main(const int argc, const char * const argv[]) { using namespace std; using namespace JPP; typedef JRange JRange_t; JServer server; string usr; string pwd; string cookie; vector inputFile; string outputFile; JRange_t range; double fraction; string testType; int debug; try { JParser<> zap("Auxiliary program to compute the PMT thresholds according to the small time-over-threshold fraction."); zap['s'] = make_field(server) = getServernames(); zap['u'] = make_field(usr) = ""; zap['!'] = make_field(pwd) = ""; zap['C'] = make_field(cookie) = ""; zap['f'] = make_field(inputFile, "list of file names (output of JCalibrateToT)"); zap['o'] = make_field(outputFile, "output file for JSon") = ""; 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['T'] = make_field(testType, "test type") = "TH-TUNING-SEA-v1"; zap['d'] = make_field(debug, "debug.") = 2; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } if (inputFile.empty()) { FATAL("No input files."); } JDateAndTime timer; JPersons person; try { JDB::reset(usr, pwd, cookie); ResultSet& rs = getResultSet(getTable(), getSelector(JDB::get()->User())); rs >> person; rs.Close(); } catch(const exception& error) { FATAL(error.what() << endl); } const JDAQHeader header = getDAQHeader(inputFile[0].c_str()); const int ID = header.getDetectorID(); typedef map module_type; typedef map detector_type; detector_type detector; try { ResultSet& rs = getResultSet(getTable(), getSelector(ID)); for (JDetectorIntegration parameters; rs >> parameters; ) { detector[parameters.DOMID][parameters.CABLEPOS] = parameters; } rs.Close(); } catch(const exception& error) { FATAL(error.what() << endl); } struct parameters_type { int threshold; double signal; double noise; }; typedef vector data_type; typedef map map_type; map_type data; vector runs; for (vector::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) { const JDAQHeader header = getDAQHeader(i->c_str()); ASSERT(header.getDetectorID() == ID, "file: " << *i); runs.push_back(header.getRunNumber()); const JPMTThreshold getPMTThreshold(header.getDetectorID(), header.getRunNumber()); TFile* in = TFile::Open(i->c_str(), "exist"); if (in == NULL || !in->IsOpen()) { FATAL("File: " << *i << " not opened." << endl); } for (detector_type::const_iterator module = detector.begin(); module != detector.end(); ++module) { if (!module->second.empty()) { TH2D* h2s = (TH2D*) in->Get(MAKE_CSTRING(module->first << _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 JUPI_t upi = module->second.at(ix-1).PMTUPI; data[upi].push_back({ getPMTThreshold(upi).value, signal, noise}); } } else { WARNING("No histogram for module " << module->first << "; skip." << endl); } } } in->Close(); } JPMTThresholdCalibration calibration; for (map_type::iterator i = data.begin(); i != data.end(); ++i) { sort(i->second.begin(), i->second.end(), make_comparator(¶meters_type::threshold)); int threshold = i->second.begin()->threshold; for (data_type::const_iterator p = i->second.begin(); p != i->second.end(); ++p) { if (p->noise <= fraction * p->signal) { threshold = p->threshold; break; } } if (debug >= debug_t || threshold > i->second.begin()->threshold) { cout << "PMT " << left << setw(32) << i->first << " -> " << right << setw(3) << threshold << " (" << i->second.begin()->threshold << ")" << endl; for (data_type::const_iterator p = i->second.begin(); p != i->second.end(); ++p) { DEBUG(setw(3) << p->threshold << ' ' << FIXED(7,0) << p->noise << "/" << FIXED(7,0) << p->signal << ' ' << (p->noise <= fraction * p->signal) << endl); } } calibration.push_back(JPMTThresholdCalibration_t(i->first.toString(), OK_t, threshold, runs)); } if (outputFile != "") { json js; js[User_t] = person.LOGIN; js[Location_t] = person.LOCATIONID; js[Start_t + Time_t] = timer.toString(); js[End_t + Time_t] = timer().toString(); js[Test_t + Type_t] = testType; js[Tests_t] = json(calibration); ofstream out(outputFile.c_str()); out << setw(2) << setprecision(8); out << js; out.close(); } }