#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 "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" #include "JLang/JToken.hh" #include "JAAnet/JHead.hh" #include "JAAnet/JHeadToolkit.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JAAnet/JEvtWeightFactorMupage.hh" #include "JAAnet/JEvtWeightFactorFunction.hh" #include "JSupport/JMeta.hh" #include "JSupport/JSupport.hh" #include "JSupport/JFileRecorder.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JEvtWeightFileScanner.hh" #include "JGizmo/JGizmoToolkit.hh" /** * \file * Program for reweighting mupage data according to a specifiable TFormula. * * The TFormula may take any number of parameters, but is limited to the physical variables\n * defined in `JEvtWeightFactorMupage`. The indices of the physical variables are defined in\n * the enum `JEvtWeightFactorMupage::variables`. * * \author bjung */ int main(int argc, char **argv) { using namespace std; using namespace JPP; typedef typename JTYPELIST::typelist typelist; typedef JToken<';'> JToken_t; JMultipleFileScanner_t inputFiles; JFileRecorder outputFile; string formula; vector parameters; int debug; try { JParser<> zap; zap['f'] = make_field(inputFiles); zap['o'] = make_field(outputFile); zap['F'] = make_field(formula, "Reweighting formula, e.g.: \"[0] + [1]*x[0]\"") = ""; zap['@'] = make_field(parameters, "Parameter values, e.g.: \"p0 = 0.5; p1 = 0.2;\"") = JPARSER::initialised(); zap['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } // Create file scanner JEvtWeightFileScanner<> scanner(inputFiles); JHead& header = scanner.getHeader(); header.createUUID(); if (!is_mupage(header)) { FATAL("MC-header is incompatible with MUPAGE."); } // Define reweighting function JEvtWeightFactorMupage reweighter("reweighter", formula.c_str()); if (reweighter.GetNdim() > JEvtWeightFactorMupage::NUMBER_OF_VARIABLES) { FATAL("Invalid formula " << formula << " (number of dimensions must be < " << JEvtWeightFactorMupage::NUMBER_OF_VARIABLES << "; see `JAANET::JEvtWeightFactorMupage`)."); } if (size_t(reweighter.GetNpar()) != parameters.size()) { FATAL("Not all parameters have been specified."); } // Read parameter settings try { for (vector::const_iterator i = parameters.begin(); i != parameters.cend(); ++i) { const int index = getParameter(*i); reweighter.SetParameter(index, getValue(*i, 0)); } } catch (JLANG::JParseError& error) { FATAL(error << endl); } // Set reweighting factors for mupage-files if (!scanner.setEvtWeightFactor(TRACK_TYPE_MUON, make_weightFactor(reweighter)) || !scanner.setEvtWeightFactor(TRACK_TYPE_ANTIMUON, make_weightFactor(reweighter))) { FATAL("Setting of reweighting factor was unsuccessful."); } // Store events with new weights outputFile.open(); Head head; copy(header, head); outputFile.put(head); outputFile.put(JMeta(argc, argv)); STATUS(RIGHT(15) << "event" << RIGHT(15) << "weight" << endl); while (scanner.hasNext()) { const Evt* event = scanner.next(); if (event != NULL) { 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() << "!"); } Evt out = *event; uuid_copy(out.header_uuid, header.UUID.uuid); if (out.w.size() <= WEIGHTLIST_RESCALED_EVENT_RATE) { out.w.resize(WEIGHTLIST_RESCALED_EVENT_RATE + 1, 0.0); } out.w.at(WEIGHTLIST_NORMALISATION) = scanner.getNormalisation(*event); out.w.at(WEIGHTLIST_RESCALED_EVENT_RATE) = weight; outputFile.put(out); } else { WARNING("Event " << scanner.getCounter() << " is empty; skip."); } } JMultipleFileScanner::typelist> io(inputFiles); io >> outputFile; outputFile.close(); return 0; }