#include #include "JDAQ/JDAQEventIO.hh" #include "Jeep/JParser.hh" #include "JSupport/JSupport.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JSupportToolkit.hh" #include "JTrigger/JHitL0.hh" #include "JTrigger/JBuildL0.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JModuleRouter.hh" #include "JROOT/JManager.hh" #include "JPhysics/JDispersion.hh" #include "JTrigger/JPMTSelector.hh" #include "JTools/JCombinatorics.hh" #include "JMath/JMathToolkit.hh" #include "JNanobeacon.hh" #include "TH1.h" #include "TH2.h" #include "TFile.h" using namespace std; using namespace JPP; using namespace JDETECTOR; using namespace KM3NETDAQ; using namespace JNANOBEACON; int main(int argc , char** argv){ JMultipleFileScanner inputFiles; JLimit_t& numberOfEvents = inputFiles.getLimit(); string detectorFile; string outFile; JDetector detector; double depth; double wavelength; JPMTSelector selector; JPMTSelector default_pmts; default_pmts.push_back(JPMTIdentifier_t(-1,0)); default_pmts.push_back(JPMTIdentifier_t(-1,3)); default_pmts.push_back(JPMTIdentifier_t(-1,4)); try { JParser<> zap; zap['f'] = make_field(inputFiles ); zap['a'] = make_field(detectorFile ); zap['n'] = make_field(numberOfEvents ) = JLimit::max(); zap['o'] = make_field(outFile ) = "out.root"; zap['h'] = make_field(depth , "detector depth [m]" ) = 2500.0; zap['l'] = make_field(wavelength , "LED wavelength [nm]") = 400.0; zap['T'] = make_field(selector ) = default_pmts; // for example '-1 0 -1 3 -1 4', to select PMTs 0, 3 and 4 on every DOM zap(argc,argv); } catch(const exception &error) { ERROR(error.what() << endl); } typedef JHitL0 JHit_t; typedef vector JHitBuffer; JBuildL0 builder; JHitBuffer triggerBuffer; JHitBuffer snapshotBuffer; load(detectorFile , detector); JModuleRouter moduleRouter(detector); JDispersion jd (0.1 * depth); double cw = getSpeedOfLight() / jd.getIndexOfRefractionGroup(wavelength); int num_of_floors = getNumberOfFloors(detector); JCombinatorics c(num_of_floors); int npairs = c.getNumberOfPairs(); c.sort(JNANOBEACON::comparepair); JManager < int, TH2D >* timeDifferences; timeDifferences = new JManager < int, TH2D > ( new TH2D("%", "", npairs, 0.5, npairs+0.5, 501, -250.5, 250.5) ); while(inputFiles.hasNext()){ JDAQEvent* event = inputFiles.next(); triggerBuffer .clear(); snapshotBuffer.clear(); builder(*event, moduleRouter, true , back_inserter(snapshotBuffer)); builder(*event, moduleRouter, false , back_inserter( triggerBuffer)); JHitBuffer::iterator __end = partition(triggerBuffer.begin(), triggerBuffer.end(), selector); if(triggerBuffer.begin() != __end){ sort(triggerBuffer.begin(), __end, less()); JHitBuffer::iterator triggeredHit = triggerBuffer.begin(); if( moduleRouter.hasModule( triggeredHit -> getModuleID() ) ){ const JModule& triggeredModule = moduleRouter.getModule( triggeredHit -> getModuleID() ); for(JHitBuffer::const_iterator snapshotHit = snapshotBuffer.begin(); snapshotHit != snapshotBuffer.end(); ++snapshotHit){ if( moduleRouter.hasModule( snapshotHit -> getModuleID() ) ){ const JModule& snapshotModule = moduleRouter.getModule( snapshotHit -> getModuleID() ); if( triggeredModule.getString() == snapshotModule.getString() ){ //double ToF = triggeredHit->getDistance(*snapshotHit) / cw; double Dt = snapshotHit->getT() - triggeredHit->getT(); double ToF = getDistance( triggeredHit->getPosition(), snapshotHit->getPosition() ) / cw; //indexes are pairs running from 0-17 between 0-152, while floors are 1-18 and bins are 1-153 int xbin = c.getIndex( triggeredModule.getFloor()-1 , snapshotModule.getFloor()-1 ) + 1; if( triggeredModule.getFloor() < snapshotModule.getFloor() ){ (*timeDifferences)[triggeredModule.getString()]->Fill(xbin,Dt-ToF); } } } else { FATAL("JModuleRouter trying to access non existing identifier: " << snapshotHit->getModuleID() ); } } } else { FATAL("JModuleRouter trying to access non existing identifier: " << triggeredHit->getModuleID() ); } } } TFile output(outFile.c_str() , "recreate") ; timeDifferences -> Write(output); }