#include #include #include #include #include "km3net-dataformat/definitions/weightlist.hh" #include "km3net-dataformat/offline/MultiHead.hh" #include "km3net-dataformat/offline/Head.hh" #include "km3net-dataformat/offline/Evt.hh" #include "JDAQ/JDAQEventIO.hh" #include "JDAQ/JDAQTimesliceIO.hh" #include "JDAQ/JDAQSummarysliceIO.hh" #include "JTrigger/JTriggerParameters.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" #include "JOscProb/JOscParameters.hh" #include "JOscProb/JOscProbInterpolator.hh" #include "JOscProb/JOscProbHelper.hh" #include "JAAnet/JHead.hh" #include "JAAnet/JMultiHead.hh" #include "JAAnet/JHeadToolkit.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JAAnet/JParticleTypes.hh" #include "JAAnet/JEvtCategoryMap.hh" #include "JReconstruction/JEvt.hh" #include "JSupport/JMeta.hh" #include "JSupport/JSupport.hh" #include "JSupport/JFileRecorder.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JEvtWeightFileScannerSet.hh" #include "JSupport/JMonteCarloFileSupportkit.hh" #include "JGizmo/JGizmoToolkit.hh" /** * \file * Program for reweighting Monte-Carlo simulations. * * \author bjjung */ int main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace KM3NETDAQ; typedef JTYPELIST::typelist, JFIT::JEvt>::typelist typelist; JMultipleFileScanner_t inputFiles; JFileRecorder outputFile; string oscProbTable; JOscParameters oscParameters; JEvtWeightFactorMap factors{JOscProbInterpolator<>()}; int debug; try { JParser<> zap; zap['f'] = make_field(inputFiles); zap['o'] = make_field(outputFile); zap['P'] = make_field(oscProbTable) = JPARSER::initialised(); zap['@'] = make_field(oscParameters) = JPARSER::initialised(); zap['%'] = make_field(factors) = JPARSER::initialised(); zap['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } // Load oscillation probability table if (!oscProbTable.empty()) { JOscProbInterpolator<>& interpolator = static_cast&>(factors.oscProb.getOscProbInterface()); interpolator.load(oscProbTable.c_str()); interpolator.set (oscParameters); } // Reweight events JEvtWeightFileScannerSet<> scanners(inputFiles); if (scanners.setEvtWeightFactor(factors) == 0) { FATAL("No weight factors could be set." << endl); } JMultiHead eventHeaders = getMultiHeader(inputFiles); JHead commonHeader = scanners.begin()->getHeader(); outputFile.open(); outputFile.put(JMeta(argc, argv)); STATUS(RIGHT(15) << "event" << RIGHT(15) << "weight" << endl); for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) { JHead header = scanner->getHeader(); header.createUUID(); // Ensure UUID is set commonHeader = commonHeader.getMatch(header); eventHeaders.merge(header); while (scanner->hasNext()) { Evt* event = scanner->next(); if (event == NULL) { WARNING("Event " << scanner->getCounter() << " is empty; skip."); continue; } const double weight = scanner->getWeight(*event); STATUS(RIGHT (15) << scanner->getCounter() << SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl); if (!isfinite(weight)) { WARNING("Non-finite reweighting factor " << weight << " for event " << scanner->getCounter() << "!"); } if (event->w.size() <= WEIGHTLIST_RESCALED_EVENT_RATE) { event->w.resize(WEIGHTLIST_RESCALED_EVENT_RATE + 1, 0.0); } const JUUID headerUUID = (uuid_is_null(event->header_uuid) == 0 ? // is valid eventHeaders.getHeaderUUID(*event) : header.UUID); event->w[WEIGHTLIST_NORMALISATION] = eventHeaders.getNormalisation(headerUUID); event->w.at(WEIGHTLIST_RESCALED_EVENT_RATE) = weight; uuid_copy(event->header_uuid, headerUUID.uuid); outputFile.put(*event); } } Head newHead; copy(commonHeader, newHead); outputFile.put(newHead); outputFile.put(static_cast(eventHeaders)); JMultipleFileScanner::typelist> io(inputFiles); io >> outputFile; outputFile.close(); return 0; }