#include #include #include #include #include "TROOT.h" #include "TFile.h" #include "TH1D.h" #include "TH2D.h" #include "TGraph.h" #include "JROOT/JGraph.hh" #include "JROOT/JRootToolkit.hh" #include "JROOT/JManager.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JMeta.hh" #include "JSupport/JSupport.hh" #include "JAcoustics/JEvt.hh" #include "JAcoustics/JSupport.hh" #include "JDynamics/JDynamics.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * Example program to check contents of acoustic events. * * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; JMultipleFileScanner inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string detectorFile; string outputFile; double Tmax_s; int debug; try { JParser<> zap("Example program to check contents of acoustic events."); zap['f'] = make_field(inputFile, "output of JKatoomba[.sh]"); zap['n'] = make_field(numberOfEvents) = JLimit_t::max(); zap['a'] = make_field(detectorFile); zap['o'] = make_field(outputFile) = "narrabri.root"; zap['T'] = make_field(Tmax_s); zap['d'] = make_field(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); } JDynamics dynamics(detector, Tmax_s); STATUS("loading input from file(s) " << inputFile << "... " << flush); dynamics.load(inputFile); STATUS("OK" << endl); map ZO; map ZA; Double_t xmin = numeric_limits::max(); Double_t xmax = numeric_limits::lowest(); for (JDynamics::JPosition::const_iterator string = dynamics.position.begin(); string != dynamics.position.end(); ++string) { if (!string->second.empty()) { xmin = min(xmin, string->second.getXmin()); xmax = max(xmax, string->second.getXmax()); } } const JFormat_t format(4, 0, std::ios_base::fmtflags(), '0'); JManager H2(new TH2D ("[%].tilt", NULL, 300, -3.0, +3.0, 300, -3.0, +3.0), '%', format); JManager HT(new TH1D ("[%].time", NULL, 200, 0.0, 4.0), '%', format); JManager HO(new TH1D ("H[%].orientation", NULL, 1000, xmin, xmax), '%', format); JManager HA(new TH1D ("H[%].amplitude", NULL, 1000, xmin, xmax), '%', format); for (JDynamics::JPosition::const_iterator string = dynamics.position.begin(); string != dynamics.position.end(); ++string) { TH1D* ho = HO[string->first]; TH1D* ha = HA[string->first]; for (Int_t i = 1; i <= HO->GetXaxis()->GetNbins(); ++i) { const Double_t x = HO->GetXaxis()->GetBinCenter(i); try { const JMODEL::JString tilt = string->second(x); ho->SetBinContent(i, tilt.getAngle()); ha->SetBinContent(i, tilt.getLength()); } catch(const exception& error) { //ERROR("Error string " << setw(4) << string->first << " at " << FIXED(20,5) << x << ' ' << error.what() << endl); } } } for (JDynamics::JPosition::const_iterator string = dynamics.position.begin(); string != dynamics.position.end(); ++string) { if (string->second.size() > 1) { TH2D* h2 = H2[string->first]; TH1D* ht = HT[string->first]; for (typename JDynamics::JPosition::function_type::const_iterator q = string->second.begin(), p = q++; q != string->second.end(); ++p, ++q) { const double t1 = q->getX() - p->getX(); HT->Fill(log10(t1)); ht->Fill(log10(t1)); const double tx = 600.0e3 * (q->getY().tx - p->getY().tx) / t1; // [mrad/10 min] const double ty = 600.0e3 * (q->getY().ty - p->getY().ty) / t1; // [mrad/10 min] if (t1 > 0 && t1 < Tmax_s) { H2->Fill(tx, ty); h2->Fill(tx, ty); } } for (typename JDynamics::JPosition::function_type::const_iterator i = string->second.begin(); i != string->second.end(); ++i) { ZO[string->first].put(i->getX(), i->getY().getAngle()); ZA[string->first].put(i->getX(), i->getY().getLength()); } } } TFile out(outputFile.c_str(), "recreate"); out << *H2 << H2 << *HT << HT << HO << HA; for (map::const_iterator i = ZO.begin(); i != ZO.end(); ++i) { out << JGraph(i->second, MAKE_CSTRING("G[" << FILL(4,'0') << i->first << "].orientation")); } for (map::const_iterator i = ZA.begin(); i != ZA.end(); ++i) { out << JGraph(i->second, MAKE_CSTRING("G[" << FILL(4,'0') << i->first << "].amplitude")); } out.Write(); out.Close(); }