#include #include #include #include #include "TRandom3.h" #include "km3net-dataformat/offline/Head.hh" #include "km3net-dataformat/offline/MultiHead.hh" #include "km3net-dataformat/offline/Evt.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JMonteCarloFileSupportkit.hh" #include "JSupport/JFileRecorder.hh" #include "JSupport/JSupport.hh" #include "JSupport/JMeta.hh" #include "JPhysics/JGeane.hh" #include "JPhysics/JGeanz.hh" #include "JPhysics/JConstants.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * * Auxiliary program to manipulate MUPAGE data. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; JMultipleFileScanner inputFile; JFileRecorder::typelist> outputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); double P; double u; double E; UInt_t seed; int debug; try { JParser<> zap("Auxiliary program to manipulate MUPAGE data."); zap['f'] = make_field(inputFile); zap['o'] = make_field(outputFile) = "mupage.root"; zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['P'] = make_field(P, "Keep muons in event with given probability (<= 1, 1 == all)") = 1.0; zap['u'] = make_field(u, "Reject events as a function of multiplicity (>= 0, 0 == none)") = 0.0; zap['E'] = make_field(E, "Energy scaling factor") = 1.0; zap['S'] = make_field(seed) = 0; zap['d'] = make_field(debug) = 2; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } gRandom->SetSeed(seed); outputFile.open(); if (!outputFile.is_open()) { FATAL("Error opening file " << outputFile << endl); } outputFile.put(JMeta(argc, argv)); outputFile.put(*gRandom); JMultipleFileScanner io(inputFile); io >> outputFile; while (inputFile.hasNext()) { STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); Evt* event = inputFile.next(); // randomly reject muons in event for (vector::iterator i = event->mc_trks.begin(); i != event->mc_trks.end(); ) { if (gRandom->Rndm() < P) ++i; else i = event->mc_trks.erase(i); } // apply energy loss for (vector::iterator i = event->mc_trks.begin(); i != event->mc_trks.end(); ) { i->E = gWater(i->E, i->t * getSpeedOfLight()); i->E *= E; if (i->E > MASS_MUON) ++i; else i = event->mc_trks.erase(i); } // randonly reject events as a function of muon multiplicity if (gRandom->Rndm() < pow(event->mc_trks.size(),-u)) { outputFile.put(*event); } } STATUS(endl); outputFile.put(*gRandom); outputFile.close(); }