#include "km3net-dataformat/offline/Head.hh" #include "km3net-dataformat/offline/Evt.hh" #include "km3net-dataformat/offline/Trk.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JAAnet/JParticleTypes.hh" #include "JSupport/JSupport.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JMonteCarloFileSupportkit.hh" #include "JGeometry2D/JVector2D.hh" #include "JGeometry2D/JCircle2D.hh" #include "JGeometry3D/JCylinder3D.hh" #include "JPhysics/JConstants.hh" #include "JPhysics/JGeane.hh" #include "JSirene/JVisibleEnergyToolkit.hh" /** * \file * * Program to test visible energy toolkit. * \author bjung */ int main(int argc, char** argv) { using namespace std; using namespace JPP; JMultipleFileScanner inputFile; double epsilon; int debug; try { JParser<> zap("Program to test visible energy toolkit"); zap['f'] = make_field(inputFile); zap['e'] = make_field(epsilon) = 1e-4; zap['d'] = make_field(debug) = 2; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } unsigned int n = 0; unsigned int N = 0; const Head& header = getHeader(inputFile); const JCylinder3D can = get(header); STATUS(LEFT(15) << "event number" << RIGHT(15) << "E0 [GeV]" << RIGHT(15) << "Evis [GeV]" << endl); for (JMultipleFileScanner& in = inputFile; in.hasNext(); ) { const Evt& evt = *(in.next()); const double E0 = getE0(evt); const double Evis = getVisibleEnergy(evt, can); STATUS(RIGHT(15) << in.getCounter() << FIXED(15,3) << E0 << FIXED(15,3) << Evis << '\r'); DEBUG(endl); if (Evis > 0.0 && Evis > E0) { STATUS(endl); WARNING("Evis > E0 for event " << in.getCounter() << endl); DEBUG(RIGHT(15) << "track PDG ID" << RIGHT(37) << "position" << RIGHT(35) << "direction" << RIGHT(15) << "E [GeV]" << RIGHT(15) << "Evis [GeV]" << endl); for (vector::const_iterator i = evt.mc_trks.cbegin(); i != evt.mc_trks.cend(); ++i) { DEBUG(RIGHT(15) << i->type << RIGHT(20) << i->pos << RIGHT(20) << i->dir << RIGHT(15) << i->E << RIGHT(15) << getVisibleEnergy(*i, can) << endl); } ++N; } ++n; } const double fraction = N / ((double) n); ASSERT(fraction < epsilon); }