#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 "Jeep/JProperties.hh" #include "JLang/JToken.hh" #include "JAAnet/JHead.hh" #include "JAAnet/JMultiHead.hh" #include "JAAnet/JHeadToolkit.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JAAnet/JParticleTypes.hh" #include "JAAnet/JEvtWeightFactorGSeaGen.hh" #include "JAAnet/JEvtWeightFactorFunction.hh" #include "JAAnet/JEvtWeightFactorMultiParticle.hh" #include "JAAnet/JAtmosphericNeutrinoFlux.hh" #include "JAAnet/JFlux.hh" #include "JOscProb/JOscParameters.hh" #include "JOscProb/JOscProbInterpolator.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" namespace { /** * Neutrino labels for command line interface. */ const std::map nuLabels { { JAANET::TRACK_TYPE_ANTINUE, "anue" }, { JAANET::TRACK_TYPE_NUE, "nue" }, { JAANET::TRACK_TYPE_ANTINUMU, "anumu" }, { JAANET::TRACK_TYPE_NUMU, "numu" }, { JAANET::TRACK_TYPE_ANTINUTAU, "anutau" }, { JAANET::TRACK_TYPE_NUTAU, "nutau" } }; } /** * \file * Program for reweighting gSeaGen data according to a specifiable TFormula. * * One TFormula can be specified for each neutrino type. * The TFormula may take any number of parameters, but is limited to the physical variables\n * defined in `JEvtWeightFactorGSeaGen`. The indices of the physical variables are defined in\n * the enum `JEvtWeightFactorGSeaGen::variables`.\n * Parameter settings can be specified behind the TFormula separated by semicolons, e.g.:\n * -\#anumu="[0]*(x[1]-[1]); p0=1; p1=2.5;" * * \author bjjung */ int main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace KM3NETDAQ; typedef JTYPELIST::typelist, JFIT::JEvt>::typelist typelist; typedef JEvtWeightFactorFunction JFluxFunction_t; JMultipleFileScanner_t inputFiles; JFileRecorder outputFile; string oscProbTable; JOscParameters oscParameters; map formulae; int debug; try { JProperties formulaMap(JEquationParameters("=","\n\r","./","#")); for (map::const_iterator i = nuLabels.cbegin(); i != nuLabels.cend(); ++i) { formulaMap[i->second] = formulae[i->first]; } JParser<> zap; zap['f'] = make_field(inputFiles); zap['o'] = make_field(outputFile); zap['P'] = make_field(oscProbTable); zap['@'] = make_field(oscParameters) = JPARSER::initialised(); zap['#'] = make_field(formulaMap) = JPARSER::initialised(); zap['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } // Create flux function JFluxMultiParticle fluxFunction; const JFluxFunction_t atmFlux = make_atmosphericNeutrinoFluxFunction(oscProbTable.c_str(), oscParameters); for (map::iterator i = formulae.begin(); i != formulae.end(); ++i) { const TString expression = i->second.GetExpFormula(); if (!expression.IsNull() && i->second.IsValid() && i->second.GetNdim() <= JEvtWeightFactorGSeaGen::NUMBER_OF_VARIABLES) { i->second.configure(atmFlux); fluxFunction.insert(i->first, make_fluxFunction(i->second)); } else if (expression.IsNull()) { fluxFunction.insert(i->first, atmFlux); } else { FATAL("Invalid formula " << expression << " (number of dimensions must be < " << JEvtWeightFactorGSeaGen::NUMBER_OF_VARIABLES << "; see `JAANET::JEvtWeightFactorGSeaGen`)."); } } // Reweight events outputFile.open(); outputFile.put(JMeta(argc, argv)); JEvtWeightFileScannerSet<> scanners(inputFiles); scanners.setFlux(fluxFunction); JMultiHead eventHeaders = getMultiHeader(inputFiles); JHead commonHeader = scanners.begin()->getHeader(); 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; }