#include #include #include #include #include "km3net-dataformat/offline/Head.hh" #include "km3net-dataformat/offline/Evt.hh" #include "JLang/JException.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JAAnet/JEvtWeightFactorFunction.hh" #include "JAAnet/JEvtWeightFactorMultiParticle.hh" #include "JAAnet/JFlux.hh" #include "JSupport/JLimit.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JEvtWeightFileScannerSet.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" #include "Jeep/JProperties.hh" namespace { using JLANG::JValueOutOfRange; using JLANG::JNullPointerException; /** * Set flux of given event to zero. * * \param evt event * \return zero [GeV * m^-2 * sr^-1 * s^-1] */ inline double zeroFlux(const Evt& evt) { return 0.0; } /** * Example function object for constant flux. */ struct JFlatFlux { /** * Default constructor. */ JFlatFlux() : flux(0.0) {} /** * Constructor. * * \param flux [GeV * m^-2 * sr^-1 * s^-1] */ JFlatFlux(const double flux) : flux(flux) {} /** * Get flux of given event. * * \param evt event * \return flux [GeV * m^-2 * sr^-1 * s^-1] */ double operator()(const Evt& evt) const { return flux; } /** * Stream input. * * \param in input stream * \param object uniform flux * \return input stream */ friend inline std::istream& operator>>(std::istream& in, JFlatFlux& object) { in >> object.flux; return in; } /** * Write power-law neutrino fluxes to output stream. * * \param out output stream * \param object uniform * \return output stream */ friend inline std::ostream& operator<<(std::ostream& out, const JFlatFlux& object) { const JFormat format(out); out << FIXED(5,3) << object.flux; return out; } double flux; //!< flux [GeV * m^-2 * sr^-1 * s^-1] }; /** * Example function object for computing power-law flux. */ struct JPowerLawFlux { /** * Default constructor. */ JPowerLawFlux() : normalisation(1.0), spectralIndex(0.0) {} /** * Constructor. * * \param normalisation normalisation * \param spectralIndex spectral index */ JPowerLawFlux(const double normalisation, const double spectralIndex) : normalisation(normalisation), spectralIndex(spectralIndex) {} /** * Get flux of given event. * * \param evt event * \return flux [GeV * m^-2 * sr^-1 * s^-1] */ double operator()(const Evt& evt) const { using namespace std; using namespace JPP; if (has_neutrino(evt)) { const Trk& neutrino = get_neutrino(evt); return normalisation * pow(neutrino.E, -spectralIndex); } else { THROW(JNullPointerException, "JPowerLawFlux::operator(): No neutrino for event " << evt.id << '.' << endl); } } /** * Stream input. * * \param in input stream * \param object power-law flux * \return input stream */ inline friend std::istream& operator>>(std::istream& in, JPowerLawFlux& object) { return in >> object.normalisation >> object.spectralIndex; } /** * Write power-law parameters to output stream. * * \param out output stream * \param object power-law flux * \return output stream */ inline friend std::ostream& operator<<(std::ostream& out, const JPowerLawFlux& object) { const JFormat format(out); out << FIXED(5,3) << object.normalisation << ' ' << FIXED(5,3) << object.spectralIndex; return out; } double normalisation; //!< normalisation [GeV * m^-2 * sr^-1 * s^-1] double spectralIndex; //!< spectral index >= 0 }; } /** * \file * Example program for scanning event-weights of Monte Carlo files\n * containing a given set of primaries. * * The list of possible options for the flux includes: *
 *    -@ zero      \
 *    -@ flat      \  \
 *    -@ powerlaw  \  \  \
 * 
* * \author bjung */ int main(int argc, char **argv) { using namespace std; using namespace JPP; JMultipleFileScanner_t inputFiles; vector zeroFluxes; map flatFluxes; map powerlawFluxes; int debug; try { JProperties fluxMaps; fluxMaps["zero"] = zeroFluxes; fluxMaps["flat"] = flatFluxes; fluxMaps["powerlaw"] = powerlawFluxes; JParser<> zap; zap['f'] = make_field(inputFiles); zap['@'] = make_field(fluxMaps) = JPARSER::initialised(); zap['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } // Sanity checks for (map::const_iterator i = powerlawFluxes.cbegin(); i != powerlawFluxes.cend(); ++i) { if (abs(i->first) != TRACK_TYPE_NUE && abs(i->first) != TRACK_TYPE_NUMU && abs(i->first) != TRACK_TYPE_NUTAU) { FATAL("Particle type " << i->first << " does not correspond to a neutrino." << endl); } } // Create multi-particle flux function JFluxMultiParticle multiFlux; for (vector::const_iterator i = zeroFluxes.cbegin(); i != zeroFluxes.cend(); ++i) { multiFlux.insert(*i, make_fluxFunction(zeroFlux)); } for (map::const_iterator i = flatFluxes.cbegin(); i != flatFluxes.cend(); ++i) { multiFlux.insert(i->first, make_fluxFunction(i->second)); } for (map::const_iterator i = powerlawFluxes.cbegin(); i != powerlawFluxes.cend(); ++i) { multiFlux.insert(i->first, make_fluxFunction(i->second)); } // Set event weighter JEvtWeightFileScannerSet<> scanners(inputFiles); const size_t n = scanners.setFlux(multiFlux); if (n == 0) { WARNING("No file found containing all given primaries; Flux function not set." << endl); } // Scan events for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) { if (scanner->simul.size() > 0) { STATUS("Scanning " << scanner->simul[0].program << " files..." << endl); } STATUS(LEFT(15) << "event" << RIGHT(15) << "weight" << endl); while (scanner->hasNext()) { const Evt* event = scanner->next(); const double weight = scanner->getWeight(*event); STATUS(LEFT (15) << scanner->getCounter() << SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl); } } return 0; }