#include #include #include #include #include "TROOT.h" #include "TFile.h" #include "TH1D.h" #include "km3net-dataformat/offline/Head.hh" #include "km3net-dataformat/offline/Evt.hh" #include "km3net-dataformat/online/JDAQClock.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JTools/JRange.hh" #include "JMath/JMathToolkit.hh" #include "JTrigger/JHitL0.hh" #include "JPhysics/JConstants.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JSupport.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { using JTRIGGER::JHitL0; using JDETECTOR::JModuleRouter; /** * Get transformation. * * \param hit hit * \param router module router * \return hit */ inline JHitL0 getHit(const Hit& hit, const JModuleRouter& router) { using namespace JPP; const JPMT& pmt = router.getModule(hit.dom_id).getPMT(hit.channel_id); return JHitL0(JDAQPMTIdentifier(hit.dom_id, hit.channel_id), pmt.getAxis(), JHit(getTime(hit.tdc, pmt.getCalibration()), getToT(hit.tot, pmt.getCalibration()))); } } /** * \file * * Example program to analyse track fit results from Evt formatted data. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace KM3NETDAQ; typedef JTOOLS::JRange JRange_t; JMultipleFileScanner inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string detectorFile; string outputFile; JRange_t E; int debug; try { JParser<> zap("Example program to analyse track fit results from Evt formatted data."); zap['f'] = make_field(inputFile); zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['a'] = make_field(detectorFile) = ""; zap['o'] = make_field(outputFile) = "histogram.root"; zap['E'] = make_field(E) = JRange_t(); zap['d'] = make_field(debug) = 2; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } JDetector detector; if (detectorFile != "") { try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } } const JModuleRouter router(detector); TFile out(outputFile.c_str(), "recreate"); TH1D hx("hx", NULL, 100, -3.0, +2.3); // [log(deg)] TH1D hd("hd", NULL, 100, 0.0, 10.0); // [m] TH1D ht("ht", NULL, 100, -100.0, 100.0); // [ns] TH1D he("he", NULL, 100, -5.0, +5.0); // [log10(E)] TH1D h1("h1", NULL, 100, -50.0, +50.0); // [ns] while (inputFile.hasNext()) { STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); const Evt* evt = inputFile.next(); //if (has_muon(*evt) && has_reconstructed_track(*evt, rec_stages_range(JMUONBEGIN, JMUONEND))) { if (has_muon(*evt) && has_reconstructed_jppmuon(*evt)) { DEBUG(get_best_reconstructed_jppmuon(*evt).comment << endl); JTrack3E ta = getTrack(get_muon(*evt)); //JTrack3E tb = getTrack(get_best_reconstructed_track(*evt, rec_stages_range(JMUONBEGIN, JMUONEND))); JTrack3E tb = getTrack(get_best_reconstructed_jppmuon(*evt)); ta.add(getTimeSinceRTS(evt->mc_t)); ta.move(ta.getIntersection(tb), getSpeedOfLight(), gWater); tb.move(tb.getIntersection(ta), getSpeedOfLight(), gWater); if (E(ta.getE())) { hx.Fill(log10(getAngle(ta.getDirection(), tb.getDirection()))); hd.Fill((ta.getPosition() - tb.getPosition()).getLength()); ht.Fill( ta.getT() - tb.getT()); he.Fill(log10(tb.getE()/ta.getE())); } for (vector::const_iterator i = evt->hits.begin(); i != evt->hits.end(); ++i) { if (router.hasModule(i->dom_id)) { const JHitL0 hit = getHit(*i, router); h1.Fill(hit.getT() - tb.getT(hit.getPosition())); } } } } STATUS(endl); out.Write(); out.Close(); }