#include <string>
#include <iostream>
#include <iomanip>
#include <map>
#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<Evt> inputFile;
JFileRecorder<JTYPELIST<JAAnetTypes_t, JMeta, JRootTypes_t>::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<Head> 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<Trk>::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<Trk>::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();
}