#include #include #include #include #include #include "TROOT.h" #include "TFile.h" #include "TH1D.h" #include "TH2D.h" #include "JLang/JComparator.hh" #include "JLang/JComparison.hh" #include "JDB/JSupport.hh" #include "JDB/JAHRS.hh" #include "JDB/JAHRSCalibration_t.hh" #include "JDB/JAHRSToolkit.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSystem/JStat.hh" #include "JTools/JQuantile.hh" #include "JROOT/JRootToolkit.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JModuleRouter.hh" #include "JDetector/JStringRouter.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * * Program to monitor AHRS data. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; JMultipleFileScanner inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string outputFile; string detectorFile; string ahrsFile; double precision; vector Q; int qaqc; int debug; try { JParser<> zap("Program to monitor AHRS data."); zap['f'] = make_field(inputFile, "output of JConvertDB -q ahrs"); zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['o'] = make_field(outputFile, "histogram output") = ""; zap['a'] = make_field(detectorFile); zap['c'] = make_field(ahrsFile, "output of JAHRSCalibration"); zap['p'] = make_field(precision, "precision") = 1.0e-6; zap['q'] = make_field(Q, "quantile"); zap['Q'] = make_field(qaqc) = 0; zap['d'] = make_field(debug) = 2; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } const JModuleRouter router(detector); const JAHRSCalibration_t calibration(ahrsFile.c_str()); const JAHRSValidity is_valid; const floor_range range = getRangeOfFloors(detector); const JStringRouter string(detector); TH1D h0("h0", NULL, 100, 0.0, 5.0); TH1D h1("h1", NULL, 100, 0.0, 5.0); TH2D h2("h2", NULL, 100, 0.0, 5.0, 100, 0.0, 5.0); TH1D hn("hn", NULL, 100, 0.0, 1.0e4); TH2D* ha = new TH2D("ha", NULL, string.size(), -0.5, string.size() - 0.5, range.getUpperLimit(), 1 - 0.5, range.getUpperLimit() + 0.5); for (Int_t i = 1; i <= ha->GetXaxis()->GetNbins(); ++i) { ha->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1))); } for (Int_t i = 1; i <= ha->GetYaxis()->GetNbins(); ++i) { ha->GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i)); } TH2D* hb = (TH2D*) ha->Clone("hb"); TH2D* hc = (TH2D*) ha->Clone("hc"); TH2D* hd = (TH2D*) ha->Clone("hd"); typedef map map_type; map_type data; for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) { if (module->getFloor() != 0) { data[module->getID()] = 0; } } map > buffer; while (inputFile.hasNext()) { STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); const JAHRS* parameters = inputFile.next(); if (router.hasModule(parameters->DOMID)) { const JLocation& location = router.getModule(parameters->DOMID); if (location.getFloor() != 0) { const double A = sqrt(parameters->AHRS_A0 * parameters->AHRS_A0 + parameters->AHRS_A1 * parameters->AHRS_A1 + parameters->AHRS_A2 * parameters->AHRS_A2); const double H = sqrt(parameters->AHRS_H0 * parameters->AHRS_H0 + parameters->AHRS_H1 * parameters->AHRS_H1 + parameters->AHRS_H2 * parameters->AHRS_H2); h0.Fill(A); h1.Fill(H); h2.Fill(A,H); ha->Fill((double) string.getIndex(location.getString()), (double) location.getFloor()); if (is_valid(*parameters)) { hb->Fill((double) string.getIndex(location.getString()), (double) location.getFloor()); if (calibration.has(parameters->DOMID)) { hc->Fill((double) string.getIndex(location.getString()), (double) location.getFloor()); } } if (is_valid(*parameters) && calibration.has(parameters->DOMID)) buffer[parameters->DOMID].push_back(*parameters); else DEBUG("AHRS " << setw(16) << parameters->UNIXTIME << ' ' << setw(4) << parameters->DUID << ' ' << setw(2) << parameters->FLOORID << ' ' << FIXED(7,3) << parameters->AHRS_A0 << ' ' << FIXED(7,3) << parameters->AHRS_A1 << ' ' << FIXED(7,3) << parameters->AHRS_A2 << ' ' << FIXED(7,3) << parameters->AHRS_H0 << ' ' << FIXED(7,3) << parameters->AHRS_H1 << ' ' << FIXED(7,3) << parameters->AHRS_H2 << ' ' << is_valid(*parameters) << calibration.has(parameters->DOMID) << endl); } } } STATUS(endl); for (map >::iterator i = buffer.begin(); i != buffer.end(); ++i) { data[i->first] = 1; if (i->second.size() >= 2u) { sort(i->second.begin(), i->second.end(), make_comparator(&JAHRS::UNIXTIME)); for (vector::const_iterator q = i->second.begin(), p = q++; q != i->second.end(); ++p, ++q) { if ((fabs(p->AHRS_A0 - q->AHRS_A0) > precision || fabs(p->AHRS_A1 - q->AHRS_A1) > precision || fabs(p->AHRS_A2 - q->AHRS_A2) > precision) && (fabs(p->AHRS_H0 - q->AHRS_H0) > precision || fabs(p->AHRS_H1 - q->AHRS_H1) > precision || fabs(p->AHRS_H2 - q->AHRS_H2) > precision)) { data[i->first] += 1; } } } } for (map_type::const_iterator i = data.begin(); i != data.end(); ++i) { const JLocation& location = router.getModule(i->first); hn.Fill(i->second); hd->Fill((double) string.getIndex(location.getString()), (double) location.getFloor(), i->second); } if (debug >= status_t) { for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) { map_type::const_iterator p = data.find(module->getID()); if (p != data.end()) { cout << getLabel(module->getLocation()) << ' ' << setw(8) << module->getID() << ' ' << setw(6) << p->second << endl; } } } if (outputFile != "") { TFile out(outputFile.c_str(), "recreate"); out << h0 << h1 << h2 << hn; out << *ha << *hb << *hc << *hd; out.Write(); out.Close(); } JQuantile q1("", true); int n_no_data_compass = 0; for (map_type::const_iterator i = data.begin(); i != data.end(); ++i) { if (i->second > 0) { q1.put(i->second); } else { n_no_data_compass++; } } for (vector::const_iterator i = Q.begin(); i != Q.end(); ++i) { QAQC(' ' << FIXED(5,0) << n_no_data_compass << ' ' << FIXED(8,3) << q1.getQuantile(*i)); } QAQC(endl); return 0; }