#include #include #include #include #include #include "TROOT.h" #include "TFile.h" #include "TH1D.h" #include "TH2D.h" #include "TProfile2D.h" #include "JTools/JHashCollection.hh" #include "JTools/JHashMap.hh" #include "JTools/JRange.hh" #include "JTools/JQuantile.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JTripod.hh" #include "JDetector/JTransmitter.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JTreeScanner.hh" #include "JROOT/JRootToolkit.hh" #include "JROOT/JManager.hh" #include "JLang/JComparator.hh" #include "JAcoustics/JEmitter.hh" #include "JAcoustics/JEvent.hh" #include "JAcoustics/JSupport.hh" #include "Jeep/JContainer.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { using JDETECTOR::JLocation; using JGEOMETRY3D::JPosition3D; /** * Auxilisry data structure for module location and position. */ struct JModule_t : public JLocation, public JPosition3D { /** * Default constructor. */ JModule_t() {} /** * Constructor. * * \param location location * \param position position */ JModule_t(const JLocation& location, const JPosition3D& position) : JLocation (location), JPosition3D(position) {} }; } /** * \file * * Example program to monitor acoustic events. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; typedef JContainer< vector > tripods_container; typedef JContainer< vector > transmitters_container; JMultipleFileScanner<> inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string outputFile; string detectorFile; tripods_container tripods; // tripods transmitters_container transmitters; // transmitters double Q; int debug; try { JParser<> zap("Example program to monitor acoustic events."); zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]"); zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['o'] = make_field(outputFile) = "monitor.root"; zap['a'] = make_field(detectorFile); zap['T'] = make_field(tripods, "tripod data") = JPARSER::initialised(); zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised(); zap['Q'] = make_field(Q) = 0.9; 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); } sort(detector.begin(), detector.end(), make_comparator(&JModule::getLocation)); JHashMap receivers; JHashMap emitters; for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) { receivers[i->getID()] = JModule_t(i->getLocation(), i->getPosition()); } for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) { emitters[i->getID()] = (i->getUTMPosition() - detector.getUTMPosition()); } for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) { try { emitters[i->getID()] = (i->getPosition() + detector.getModule(i->getLocation()).getPosition()); } catch(const exception&) { continue; // if no module available, discard transmitter } } const JHashCollection string(make_array(detector.begin(), detector.end(), &JModule::getString)); const JRange floor (make_array(detector.begin(), detector.end(), &JModule::getFloor)); JManager HA(new TH1D("H[%].size", NULL, detector.size() + 1, -0.5, detector.size() + 0.5)); JManager HB(new TH1D("H[%].overlays", NULL, 201, -0.5, 200.5)); JManager HC(new TH1D("H[%].time", NULL, 800, 0.0, 4.0)); JManager HD(new TH1D("H[%].rms", NULL, 500, 0.0, 1.0e-1)); JManager HE(new TH1D("H[%].quantile", NULL, 500, 0.0, 1.0e-1)); JManager HQ(new TH1D("H[%].quality", NULL, 200, 0.0, 8.0)); JManager HR(new TH2D("H[%].QR", NULL, 40, 1.5, 3.5, 40, 3.0, 7.0)); JManager H1(new TH1D("H[%].multiplicity", NULL, 101, -0.5, 100.5)); JManager H2(new TH2D("H[%].event", NULL, string.size(), -0.5, string.size() - 0.5, floor.getUpperLimit() + 1, - 0.5, floor.getUpperLimit() + 0.5)); JManager H4(new TProfile2D("H[%].log10(Q)", NULL, string.size(), -0.5, string.size() - 0.5, floor.getUpperLimit() + 1, - 0.5, floor.getUpperLimit() + 0.5)); for (Int_t i = 1; i <= H2->GetXaxis()->GetNbins(); ++i) { H2->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1))); } for (Int_t i = 1; i <= H2->GetYaxis()->GetNbins(); ++i) { H2->GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i-1)); } for (Int_t i = 1; i <= H4->GetXaxis()->GetNbins(); ++i) { H4->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1))); } for (Int_t i = 1; i <= H4->GetYaxis()->GetNbins(); ++i) { H4->GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i-1)); } JManager H3((TH2D*) H2->Clone("H[%].doubles")); for (JHashMap::const_iterator i = emitters.begin(); i != emitters.end(); ++i) { H2[i->first]; H3[i->first]; } typedef JTreeScanner JTreeScanner_t; counter_type counter = 0; for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) { JTreeScanner_t in(*i); JHashMap toe; for (JTreeScanner_t::iterator event = in.begin(); event != in.end() && counter != inputFile.getLimit(); ++event, ++counter) { if (counter%1000 == 0) { STATUS("event " << setw(8) << counter << '\r'); DEBUG(endl); } const JEvent evt = *event; JQuantile Q1("", true); for (JEvent::const_iterator i = evt.begin(); i != evt.end(); ++i) { Q1.put(i->getToE()); } HA[evt.getID()]->Fill((double) evt.size()); HB[evt.getID()]->Fill((double) evt.getOverlays()); if (toe.has(evt.getID())) { HC[evt.getID()]->Fill(log10(evt.begin()->getToE() - toe[evt.getID()])); } HD[evt.getID()]->Fill(Q1.getSTDev()); HE[evt.getID()]->Fill(Q1.getQuantile(Q, JQuantile::symmetric_t)); toe[evt.getID()] = evt.begin()->getToE(); TH1D* hq = HQ[evt.getID()]; TH2D* hr = HR[evt.getID()]; TH1D* h1 = H1[evt.getID()]; TH2D* h2 = H2[evt.getID()]; TH2D* h3 = H3[evt.getID()]; TProfile2D* h4 = H4[evt.getID()]; for (JEvent::const_iterator i = evt.begin(); i != evt.end(); ++i) { hq->Fill(log10(i->getQ())); } if (emitters.has(evt.getID())) { const JPosition3D& p0 = emitters[evt.getID()]; for (JEvent::const_iterator i = evt.begin(); i != evt.end(); ++i) { if (receivers.has(i->getID())) { HR->Fill(log10(p0.getDistance(receivers[i->getID()])), log10(i->getQ())); hr->Fill(log10(p0.getDistance(receivers[i->getID()])), log10(i->getQ())); } } } map > buffer; for (JEvent::const_iterator i = evt.begin(); i != evt.end(); ++i) { buffer[i->getID()].insert(i->getQ()); } for (map >::const_iterator i = buffer.begin(); i != buffer.end(); ++i) { h1->Fill((double) i->second.size()); if (receivers.has(i->first)) { const JLocation& location = receivers[i->first]; const double x = string.getIndex(location.getString()); const double y = location.getFloor(); h2->Fill(x, y, 1.0); if (i->second.size() >= 2u) { h3->Fill(x, y, 1.0); } h4->Fill(x, y, log10(*(i->second.rbegin()))); } } } } STATUS(endl); TFile out(outputFile.c_str(), "recreate"); out << HA << HB << HC << HD << HE << HQ << HR << *HR << H1 << H2 << H3 << H4; out.Write(); out.Close(); }