diff --git a/bin/rivet-nopy.cc b/bin/rivet-nopy.cc --- a/bin/rivet-nopy.cc +++ b/bin/rivet-nopy.cc @@ -1,41 +1,42 @@ #include "Rivet/AnalysisHandler.hh" #include "Rivet/AnalysisLoader.hh" //#include "HepMC/IO_GenEvent.h" #include "Rivet/Tools/RivetHepMC.hh" +#include using namespace std; int main(int argc, char** argv) { if (argc < 2) { cerr << "Usage: " << argv[0] << " [ ...]" << '\n'; cout << "Available analyses:\n"; for (const string& a : Rivet::AnalysisLoader::analysisNames()) cout << " " << a << "\n"; cout << endl; return 1; } Rivet::AnalysisHandler ah; for (int i = 2; i < argc; ++i) { ah.addAnalysis(argv[i]); } std::ifstream istr(argv[1], std::ios::in); std::shared_ptr reader = Rivet::HepMCUtils::makeReader(istr); std::shared_ptr evt = make_shared(); while(reader && Rivet::HepMCUtils::readEvent(reader, evt)){ ah.analyze(evt.get()); evt.reset(new Rivet::RivetHepMC::GenEvent()); } istr.close(); - ah.setCrossSection(1.0, 0.0); + ah.setCrossSection(make_pair(1.0, 0.0)); ah.finalize(); ah.writeData("Rivet.yoda"); return 0; } diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh --- a/include/Rivet/AnalysisHandler.hh +++ b/include/Rivet/AnalysisHandler.hh @@ -1,366 +1,370 @@ // -*- C++ -*- #ifndef RIVET_RivetHandler_HH #define RIVET_RivetHandler_HH #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Particle.hh" #include "Rivet/AnalysisLoader.hh" #include "Rivet/Tools/RivetYODA.hh" namespace Rivet { // Forward declaration and smart pointer for Analysis class Analysis; typedef std::shared_ptr AnaHandle; /// A class which handles a number of analysis objects to be applied to /// generated events. An {@link Analysis}' AnalysisHandler is also responsible /// for handling the final writing-out of histograms. class AnalysisHandler { public: /// @name Constructors and destructors. */ //@{ /// Preferred constructor, with optional run name. AnalysisHandler(const string& runname=""); /// @brief Destructor /// The destructor is not virtual, as this class should not be inherited from. ~AnalysisHandler(); //@} private: /// Get a logger object. Log& getLog() const; public: /// @name Run properties and weights //@{ /// Get the name of this run. string runName() const; /// Get the number of events seen. Should only really be used by external /// steering code or analyses in the finalize phase. size_t numEvents() const; /// @brief Access the sum of the event weights seen /// /// This is the weighted equivalent of the number of events. It should only /// be used by external steering code or analyses in the finalize phase. double sumW() const { return _eventCounter->sumW(); } /// Access to the sum of squared-weights double sumW2() const { return _eventCounter->sumW2(); } /// Names of event weight categories const vector& weightNames() const { return _weightNames; } /// Are any of the weights non-numeric? size_t numWeights() const { return _weightNames.size(); } /// Are any of the weights non-numeric? bool haveNamedWeights() const; /// Set the weight names from a GenEvent void setWeightNames(const GenEvent& ge); /// Get the index of the nominal weight-stream size_t defaultWeightIndex() const { return _defaultWeightIdx; } //@} /// @name Cross-sections //@{ /// Get the cross-section known to the handler. Scatter1DPtr crossSection() const { return _xs; } /// Set the cross-section for the process being generated. - AnalysisHandler& setCrossSection(double xs, double xserr); + void setCrossSection(pair xsec); /// Get the nominal cross-section double nominalCrossSection() const { _xs.get()->setActiveWeightIdx(_defaultWeightIdx); const YODA::Scatter1D::Points& ps = _xs->points(); if (ps.size() != 1) { string errMsg = "cross section missing when requesting nominal cross section"; throw Error(errMsg); } double xs = ps[0].x(); _xs.get()->unsetActiveWeight(); return xs; } //@} /// @name Beams //@{ /// Set the beam particles for this run AnalysisHandler& setRunBeams(const ParticlePair& beams) { _beams = beams; MSG_DEBUG("Setting run beams = " << beams << " @ " << sqrtS()/GeV << " GeV"); return *this; } /// Get the beam particles for this run, usually determined from the first event. const ParticlePair& beams() const { return _beams; } /// Get beam IDs for this run, usually determined from the first event. /// @deprecated Use standalone beamIds(ah.beams()), to clean AH interface PdgIdPair beamIds() const; /// Get energy for this run, usually determined from the first event. /// @deprecated Use standalone sqrtS(ah.beams()), to clean AH interface double sqrtS() const; /// Setter for _ignoreBeams void setIgnoreBeams(bool ignore=true); /// Setter for _skipWeights void skipMultiWeights(bool ignore=false); //@} /// @name Handle analyses //@{ /// Get a list of the currently registered analyses' names. std::vector analysisNames() const; /// Get the collection of currently registered analyses. const std::map& analysesMap() const { return _analyses; } /// Get the collection of currently registered analyses. std::vector analyses() const { std::vector rtn; rtn.reserve(_analyses.size()); for (const auto& apair : _analyses) rtn.push_back(apair.second); return rtn; } /// Get a registered analysis by name. AnaHandle analysis(const std::string& analysisname) { if ( _analyses.find(analysisname) == _analyses.end() ) throw LookupError("No analysis named '" + analysisname + "' registered in AnalysisHandler"); try { return _analyses[analysisname]; } catch (...) { throw LookupError("No analysis named '" + analysisname + "' registered in AnalysisHandler"); } } /// Add an analysis to the run list by object AnalysisHandler& addAnalysis(Analysis* analysis); /// @brief Add an analysis to the run list using its name. /// /// The actual Analysis to be used will be obtained via /// AnalysisLoader::getAnalysis(string). If no matching analysis is found, /// no analysis is added (i.e. the null pointer is checked and discarded. AnalysisHandler& addAnalysis(const std::string& analysisname); /// @brief Add an analysis with a map of analysis options. AnalysisHandler& addAnalysis(const std::string& analysisname, std::map pars); /// @brief Add analyses to the run list using their names. /// /// The actual {@link Analysis}' to be used will be obtained via /// AnalysisHandler::addAnalysis(string), which in turn uses /// AnalysisLoader::getAnalysis(string). If no matching analysis is found /// for a given name, no analysis is added, but also no error is thrown. AnalysisHandler& addAnalyses(const std::vector& analysisnames); /// Remove an analysis from the run list using its name. AnalysisHandler& removeAnalysis(const std::string& analysisname); /// Remove analyses from the run list using their names. AnalysisHandler& removeAnalyses(const std::vector& analysisnames); /// //@} /// @name Main init/execute/finalise //@{ /// Initialize a run, with the run beams taken from the example event. void init(const GenEvent& event); /// @brief Analyze the given \a event by reference. /// /// This function will call the AnalysisBase::analyze() function of all /// included analysis objects. void analyze(const GenEvent& event); /// @brief Analyze the given \a event by pointer. /// /// This function will call the AnalysisBase::analyze() function of all /// included analysis objects, after checking the event pointer validity. void analyze(const GenEvent* event); - + /// Finalize a run. This function calls the AnalysisBase::finalize() /// functions of all included analysis objects. void finalize(); //@} /// @name Histogram / data object access //@{ + /// After all subevents in an event group has been processed push + /// all histo fills to the relevant histograms. + void pushToPersistent(); + /// Read analysis plots into the histo collection (via addData) from the named file. void readData(const std::string& filename); /// Get all multi-weight Rivet analysis object wrappers vector getRivetAOs() const; /// Get a pointer to a preloaded yoda object with the given path, /// or null if path is not found. const YODA::AnalysisObjectPtr getPreload(string path) const { auto it = _preloads.find(path); if ( it == _preloads.end() ) return nullptr; return it->second; } /// Write all analyses' plots (via getData) to the named file. void writeData(const std::string& filename) const; /// Tell Rivet to dump intermediate result to a file named @a /// dumpfile every @a period'th event. If @period is not positive, /// no dumping will be done. void dump(string dumpfile, int period) { _dumpPeriod = period; _dumpFile = dumpfile; } /// Take the vector of yoda files and merge them together using /// the cross section and weight information provided in each /// file. Each file in @a aofiles is assumed to have been produced /// by Rivet. By default the files are assumed to contain /// different processes (or the same processs but mutually /// exclusive cuts), but if @a equiv if ture, the files are /// assumed to contain output of completely equivalent (but /// statistically independent) Rivet runs. The corresponding /// analyses will be loaded and their analysis objects will be /// filled with the merged result. finalize() will be run on each /// relevant analysis. The resulting YODA file can then be rwitten /// out by writeData(). If delopts is non-empty, it is assumed to /// contain names different options to be merged into the same /// analysis objects. void mergeYodas(const vector & aofiles, const vector & delopts = vector(), bool equiv = false); /// Helper function to strip specific options from data object paths. void stripOptions(YODA::AnalysisObjectPtr ao, const vector & delopts) const; //@} /// Indicate which Rivet stage we're in. /// At the moment, only INIT is used to enable booking. enum class Stage { OTHER, INIT, FINALIZE }; /// Which stage are we in? Stage stage() const { return _stage; } private: /// Current handler stage Stage _stage = Stage::OTHER; /// The collection of Analysis objects to be used. std::map _analyses; /// A vector of pre-loaded object which do not have a valid /// Analysis plugged in. map _preloads; /// A vector containing copies of analysis objects after /// finalize() has been run. vector _finalizedAOs; /// @name Run properties //@{ /// Weight names std::vector _weightNames; std::vector > _subEventWeights; size_t _numWeightTypes; // always == WeightVector.size() /// Run name std::string _runname; /// Event counter mutable CounterPtr _eventCounter; /// Cross-section known to AH Scatter1DPtr _xs; /// Beams used by this run. ParticlePair _beams; /// Flag to check if init has been called bool _initialised; /// Flag whether input event beams should be ignored in compatibility check bool _ignoreBeams; /// Flag to check if multiweights should be included bool _skipWeights; /// Current event number int _eventNumber; /// The index in the weight vector for the nominal weight stream size_t _defaultWeightIdx; /// Determines how often Rivet runs finalize() and writes the /// result to a YODA file. int _dumpPeriod; /// The name of a YODA file to which Rivet periodically dumps /// results. string _dumpFile; /// Flag to indicate periodic dumping is in progress bool _dumping; //@} private: /// The assignment operator is private and must never be called. /// In fact, it should not even be implemented. AnalysisHandler& operator=(const AnalysisHandler&); /// The copy constructor is private and must never be called. In /// fact, it should not even be implemented. AnalysisHandler(const AnalysisHandler&); }; } #endif diff --git a/include/Rivet/Tools/RivetHepMC.hh b/include/Rivet/Tools/RivetHepMC.hh --- a/include/Rivet/Tools/RivetHepMC.hh +++ b/include/Rivet/Tools/RivetHepMC.hh @@ -1,105 +1,105 @@ // -*- C++ -*- #ifndef RIVET_RivetHepMC_HH #define RIVET_RivetHepMC_HH #include #ifdef ENABLE_HEPMC_3 #include "HepMC3/HepMC3.h" #include "HepMC3/Relatives.h" #include "HepMC3/Reader.h" namespace Rivet{ namespace RivetHepMC = HepMC3; using RivetHepMC::ConstGenParticlePtr; using RivetHepMC::ConstGenVertexPtr; using RivetHepMC::Relatives; using RivetHepMC::ConstGenHeavyIonPtr; using HepMC_IO_type = RivetHepMC::Reader; using PdfInfo = RivetHepMC::GenPdfInfo; } #else #include "HepMC/GenEvent.h" #include "HepMC/GenParticle.h" #include "HepMC/HeavyIon.h" #include "HepMC/GenVertex.h" #include "HepMC/Version.h" #include "HepMC/GenRanges.h" #include "HepMC/IO_GenEvent.h" namespace Rivet{ namespace RivetHepMC = HepMC; // HepMC 2.07 provides its own #defines typedef const HepMC::GenParticle* ConstGenParticlePtr; typedef const HepMC::GenVertex* ConstGenVertexPtr; typedef const HepMC::HeavyIon* ConstGenHeavyIonPtr; /// @brief Replicated the HepMC3 Relatives syntax using HepMC2 IteratorRanges /// This is necessary mainly because of capitalisation differences class Relatives{ public: constexpr Relatives(HepMC::IteratorRange relo): _internal(relo){} constexpr HepMC::IteratorRange operator()() const {return _internal;} operator HepMC::IteratorRange() const {return _internal;} const static Relatives PARENTS; const static Relatives CHILDREN; const static Relatives ANCESTORS; const static Relatives DESCENDANTS; private: const HepMC::IteratorRange _internal; }; using HepMC_IO_type = HepMC::IO_GenEvent; using PdfInfo = RivetHepMC::PdfInfo; } #endif #include "Rivet/Tools/RivetSTL.hh" #include "Rivet/Tools/Exceptions.hh" namespace Rivet { using RivetHepMC::GenEvent; using ConstGenEventPtr = std::shared_ptr; /// @todo Use mcutils? namespace HepMCUtils{ ConstGenParticlePtr getParticlePtr(const RivetHepMC::GenParticle & gp); std::vector particles(ConstGenEventPtr ge); std::vector particles(const GenEvent *ge); std::vector vertices(ConstGenEventPtr ge); std::vector vertices(const GenEvent *ge); std::vector particles(ConstGenVertexPtr gv, const Relatives &relo); std::vector particles(ConstGenParticlePtr gp, const Relatives &relo); int uniqueId(ConstGenParticlePtr gp); int particles_size(ConstGenEventPtr ge); int particles_size(const GenEvent *ge); std::pair beams(const GenEvent *ge); std::shared_ptr makeReader(std::istream &istr, std::string * errm = 0); bool readEvent(std::shared_ptr io, std::shared_ptr evt); void strip(GenEvent & ge, const set & stripid = {1, -1, 2, -2, 3,-3, 21}); vector weightNames(const GenEvent & ge); - double crossSection(const GenEvent & ge); + pair crossSection(const GenEvent & ge); std::valarray weights(const GenEvent & ge); } } #endif diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc --- a/src/Core/AnalysisHandler.cc +++ b/src/Core/AnalysisHandler.cc @@ -1,696 +1,684 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/AnalysisHandler.hh" #include "Rivet/Analysis.hh" #include "Rivet/Tools/ParticleName.hh" #include "Rivet/Tools/BeamConstraint.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/Beam.hh" #include "YODA/IO.h" #include using std::cout; using std::cerr; namespace Rivet { AnalysisHandler::AnalysisHandler(const string& runname) : _runname(runname), _initialised(false), _ignoreBeams(false), _skipWeights(false), _defaultWeightIdx(0), _dumpPeriod(0), _dumping(false) { } AnalysisHandler::~AnalysisHandler() { } Log& AnalysisHandler::getLog() const { return Log::getLog("Rivet.Analysis.Handler"); } /// http://stackoverflow.com/questions/4654636/how-to-determine-if-a-string-is-a-number-with-c namespace { bool is_number(const std::string& s) { std::string::const_iterator it = s.begin(); while (it != s.end() && std::isdigit(*it)) ++it; return !s.empty() && it == s.end(); } } /// Check if any of the weightnames is not a number bool AnalysisHandler::haveNamedWeights() const { bool dec=false; for (unsigned int i=0;i<_weightNames.size();++i) { string s = _weightNames[i]; if (!is_number(s)) { dec=true; break; } } return dec; } void AnalysisHandler::init(const GenEvent& ge) { if (_initialised) throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!"); /// @todo Should the Rivet analysis objects know about weight names? setRunBeams(Rivet::beams(ge)); MSG_DEBUG("Initialising the analysis handler"); _eventNumber = ge.event_number(); setWeightNames(ge); if (_skipWeights) MSG_INFO("Only using nominal weight. Variation weights will be ignored."); else if (haveNamedWeights()) MSG_INFO("Using named weights"); else MSG_INFO("NOT using named weights. Using first weight as nominal weight"); _eventCounter = CounterPtr(weightNames(), Counter("_EVTCOUNT")); // Set the cross section based on what is reported by this event. - if (ge.cross_section()) { - MSG_TRACE("Getting cross section."); - double xs = ge.cross_section()->cross_section(); - double xserr = ge.cross_section()->cross_section_error(); - setCrossSection(xs, xserr); - } + if ( ge.cross_section() ) setCrossSection(HepMCUtils::crossSection(ge)); // Check that analyses are beam-compatible, and remove those that aren't const size_t num_anas_requested = analysisNames().size(); vector anamestodelete; for (const AnaHandle a : analyses()) { if (!_ignoreBeams && !a->isCompatible(beams())) { //MSG_DEBUG(a->name() << " requires beams " << a->requiredBeams() << " @ " << a->requiredEnergies() << " GeV"); anamestodelete.push_back(a->name()); } } for (const string& aname : anamestodelete) { MSG_WARNING("Analysis '" << aname << "' is incompatible with the provided beams: removing"); removeAnalysis(aname); } if (num_anas_requested > 0 && analysisNames().empty()) { cerr << "All analyses were incompatible with the first event's beams\n" << "Exiting, since this probably wasn't intentional!" << endl; exit(1); } // Warn if any analysis' status is not unblemished for (const AnaHandle a : analyses()) { if ( a->info().preliminary() ) { MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!"); } else if ( a->info().obsolete() ) { MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!"); } else if (( a->info().unvalidated() ) ) { MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!"); } } // Initialize the remaining analyses _stage = Stage::INIT; for (AnaHandle a : analyses()) { MSG_DEBUG("Initialising analysis: " << a->name()); try { // Allow projection registration in the init phase onwards a->_allowProjReg = true; a->init(); //MSG_DEBUG("Checking consistency of analysis: " << a->name()); //a->checkConsistency(); } catch (const Error& err) { cerr << "Error in " << a->name() << "::init method: " << err.what() << endl; exit(1); } MSG_DEBUG("Done initialising analysis: " << a->name()); } _stage = Stage::OTHER; _initialised = true; MSG_DEBUG("Analysis handler initialised"); } void AnalysisHandler::setWeightNames(const GenEvent& ge) { if (!_skipWeights) _weightNames = HepMCUtils::weightNames(ge); if ( _weightNames.empty() ) _weightNames.push_back(""); for ( int i = 0, N = _weightNames.size(); i < N; ++i ) if ( _weightNames[i] == "" ) _defaultWeightIdx = i; } void AnalysisHandler::analyze(const GenEvent& ge) { // Call init with event as template if not already initialised if (!_initialised) init(ge); assert(_initialised); // Ensure that beam details match those from the first event (if we're checking beams) if ( !_ignoreBeams ) { const PdgIdPair beams = Rivet::beamIds(ge); const double sqrts = Rivet::sqrtS(ge); if (!compatible(beams, _beams) || !fuzzyEquals(sqrts, sqrtS())) { cerr << "Event beams mismatch: " << PID::toBeamsString(beams) << " @ " << sqrts/GeV << " GeV" << " vs. first beams " << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl; exit(1); } } // Create the Rivet event wrapper /// @todo Filter/normalize the event here bool strip = ( getEnvParam("RIVET_STRIP_HEPMC", string("NOOOO") ) != "NOOOO" ); Event event(ge, strip); // set the cross section based on what is reported by this event. // if no cross section - MSG_TRACE("getting cross section."); - if (ge.cross_section()) { - MSG_TRACE("getting cross section from GenEvent."); - double xs = ge.cross_section()->cross_section(); - double xserr = ge.cross_section()->cross_section_error(); - setCrossSection(xs, xserr); - } + if ( ge.cross_section() ) setCrossSection(HepMCUtils::crossSection(ge)); // Won't happen for first event because _eventNumber is set in init() if (_eventNumber != ge.event_number()) { - /// @todo Can we get away with not passing a matrix? - MSG_TRACE("AnalysisHandler::analyze(): Pushing _eventCounter to persistent."); - _eventCounter.get()->pushToPersistent(_subEventWeights); - // if this is indeed a new event, push the temporary - // histograms and reset - for (const AnaHandle& a : analyses()) { - for (auto ao : a->analysisObjects()) { - MSG_TRACE("AnalysisHandler::analyze(): Pushing " << a->name() << "'s " << ao->name() << " to persistent."); - ao.get()->pushToPersistent(_subEventWeights); - } - MSG_TRACE("AnalysisHandler::analyze(): finished pushing " << a->name() << "'s objects to persistent."); - } + pushToPersistent(); - _eventNumber = ge.event_number(); + _eventNumber = ge.event_number(); - MSG_DEBUG("nominal event # " << _eventCounter.get()->_persistent[0]->numEntries()); - MSG_DEBUG("nominal sum of weights: " << _eventCounter.get()->_persistent[0]->sumW()); - MSG_DEBUG("Event has " << _subEventWeights.size() << " sub events."); - _subEventWeights.clear(); } MSG_TRACE("starting new sub event"); _eventCounter.get()->newSubEvent(); for (const AnaHandle& a : analyses()) { for (auto ao : a->analysisObjects()) { ao.get()->newSubEvent(); } } _subEventWeights.push_back(event.weights()); MSG_DEBUG("Analyzing subevent #" << _subEventWeights.size() - 1 << "."); _eventCounter->fill(); // Run the analyses for (AnaHandle a : analyses()) { MSG_TRACE("About to run analysis " << a->name()); try { a->analyze(event); } catch (const Error& err) { cerr << "Error in " << a->name() << "::analyze method: " << err.what() << endl; exit(1); } MSG_TRACE("Finished running analysis " << a->name()); } if ( _dumpPeriod > 0 && numEvents() > 0 && numEvents()%_dumpPeriod == 0 ) { MSG_INFO("Dumping intermediate results to " << _dumpFile << "."); _dumping = numEvents()/_dumpPeriod; finalize(); writeData(_dumpFile); _dumping = 0; } } void AnalysisHandler::analyze(const GenEvent* ge) { if (ge == nullptr) { MSG_ERROR("AnalysisHandler received null pointer to GenEvent"); //throw Error("AnalysisHandler received null pointer to GenEvent"); } analyze(*ge); } + void AnalysisHandler::pushToPersistent() { + if ( _subEventWeights.empty() ) return; + MSG_TRACE("AnalysisHandler::analyze(): Pushing _eventCounter to persistent."); + _eventCounter.get()->pushToPersistent(_subEventWeights); + for (const AnaHandle& a : analyses()) { + for (auto ao : a->analysisObjects()) { + MSG_TRACE("AnalysisHandler::analyze(): Pushing " << a->name() + << "'s " << ao->name() << " to persistent."); + ao.get()->pushToPersistent(_subEventWeights); + } + MSG_TRACE("AnalysisHandler::analyze(): finished pushing " + << a->name() << "'s objects to persistent."); + } + _subEventWeights.clear(); + } void AnalysisHandler::finalize() { if (!_initialised) return; MSG_INFO("Finalising analyses"); _stage = Stage::FINALIZE; // First push all analyses' objects to persistent and final MSG_TRACE("AnalysisHandler::finalize(): Pushing analysis objects to persistent."); - _eventCounter.get()->pushToPersistent(_subEventWeights); + pushToPersistent(); + + // Copy all histos to finalize versions. _eventCounter.get()->pushToFinal(); _xs.get()->pushToFinal(); - for (const AnaHandle& a : analyses()) { - for (auto ao : a->analysisObjects()) { - ao.get()->pushToPersistent(_subEventWeights); + for (const AnaHandle& a : analyses()) + for (auto ao : a->analysisObjects()) ao.get()->pushToFinal(); - } - } for (AnaHandle a : analyses()) { if ( _dumping && !a->info().reentrant() ) { if ( _dumping == 1 ) MSG_INFO("Skipping finalize in periodic dump of " << a->name() << " as it is not declared reentrant."); continue; } for (size_t iW = 0; iW < numWeights(); iW++) { _eventCounter.get()->setActiveFinalWeightIdx(iW); _xs.get()->setActiveFinalWeightIdx(iW); for (auto ao : a->analysisObjects()) ao.get()->setActiveFinalWeightIdx(iW); try { MSG_TRACE("running " << a->name() << "::finalize() for weight " << iW << "."); a->finalize(); } catch (const Error& err) { cerr << "Error in " << a->name() << "::finalize method: " << err.what() << '\n'; exit(1); } } } // Print out number of events processed const int nevts = numEvents(); MSG_INFO("Processed " << nevts << " event" << (nevts != 1 ? "s" : "")); _stage = Stage::OTHER; if ( _dumping ) return; // Print out MCnet boilerplate if (getLog().getLevel()<=20){ cout << endl; cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << endl; cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:1003.0694 (http://arxiv.org/abs/1003.0694)" << endl; } } AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname, std::map pars) { // Make an option handle. std::string parHandle = ""; for (map::iterator par = pars.begin(); par != pars.end(); ++par) { parHandle +=":"; parHandle += par->first + "=" + par->second; } return addAnalysis(analysisname + parHandle); } AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname) { // Check for a duplicate analysis /// @todo Might we want to be able to run an analysis twice, with different params? /// Requires avoiding histo tree clashes, i.e. storing the histos on the analysis objects. string ananame = analysisname; vector anaopt = split(analysisname, ":"); if ( anaopt.size() > 1 ) ananame = anaopt[0]; AnaHandle analysis( AnalysisLoader::getAnalysis(ananame) ); if (analysis.get() != 0) { // < Check for null analysis. MSG_DEBUG("Adding analysis '" << analysisname << "'"); map opts; for ( int i = 1, N = anaopt.size(); i < N; ++i ) { vector opt = split(anaopt[i], "="); if ( opt.size() != 2 ) { MSG_WARNING("Error in option specification. Skipping analysis " << analysisname); return *this; } if ( !analysis->info().validOption(opt[0], opt[1]) ) { MSG_WARNING("Cannot set option '" << opt[0] << "' to '" << opt[1] << "'. Skipping analysis " << analysisname); return *this; } opts[opt[0]] = opt[1]; } for ( auto opt: opts) { analysis->_options[opt.first] = opt.second; analysis->_optstring += ":" + opt.first + "=" + opt.second; } for (const AnaHandle& a : analyses()) { if (a->name() == analysis->name() ) { MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate"); return *this; } } analysis->_analysishandler = this; _analyses[analysisname] = analysis; } else { MSG_WARNING("Analysis '" << analysisname << "' not found."); } // MSG_WARNING(_analyses.size()); // for (const AnaHandle& a : _analyses) MSG_WARNING(a->name()); return *this; } AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) { MSG_DEBUG("Removing analysis '" << analysisname << "'"); if (_analyses.find(analysisname) != _analyses.end()) _analyses.erase(analysisname); // } return *this; } // void AnalysisHandler::addData(const std::vector& aos) { // for (const YODA::AnalysisObjectPtr ao : aos) { // string path = ao->path(); // if ( path.substr(0, 5) != "/RAW/" ) { // _orphanedPreloads.push_back(ao); // continue; // } // path = path.substr(4); // ao->setPath(path); // if (path.size() > 1) { // path > "/" // try { // const string ananame = ::split(path, "/")[0]; // AnaHandle a = analysis(ananame); // /// @todo FIXXXXX // //MultiweightAOPtr mao = ????; /// @todo generate right Multiweight object from ao // //a->addAnalysisObject(mao); /// @todo Need to statistically merge... // } catch (const Error& e) { // MSG_TRACE("Adding analysis object " << path << // " to the list of orphans."); // _orphanedPreloads.push_back(ao); // } // } // } // } void AnalysisHandler::stripOptions(YODA::AnalysisObjectPtr ao, const vector & delopts) const { string path = ao->path(); string ananame = split(path, "/")[0]; vector anaopts = split(ananame, ":"); for ( int i = 1, N = anaopts.size(); i < N; ++i ) for ( auto opt : delopts ) if ( opt == "*" || anaopts[i].find(opt + "=") == 0 ) path.replace(path.find(":" + anaopts[i]), (":" + anaopts[i]).length(), ""); ao->setPath(path); } void AnalysisHandler::mergeYodas(const vector & aofiles, const vector & delopts, bool equiv) { // Convenient typedef; typedef multimap AOMap; // Store all found weights here. set foundWeightNames; // Stor all found analyses. set foundAnalyses; // Store all analysis objects here. vector allaos; // Go through all files and collect information. for ( auto file : aofiles ) { allaos.push_back(AOMap()); AOMap & aomap = allaos.back(); vector aos_raw; try { YODA::read(file, aos_raw); } catch (...) { //< YODA::ReadError& throw UserError("Unexpected error in reading file: " + file); } for (YODA::AnalysisObject* aor : aos_raw) { YODA::AnalysisObjectPtr ao(aor); AOPath path(ao->path()); if ( !path ) throw UserError("Invalid path name in file: " + file); if ( !path.isRaw() ) continue; foundWeightNames.insert(path.weight()); // Now check if any options should be removed. for ( string delopt : delopts ) if ( path.hasOption(delopt) ) path.removeOption(delopt); path.setPath(); if ( path.analysisWithOptions() != "" ) foundAnalyses.insert(path.analysisWithOptions()); aomap.insert(make_pair(path.path(), ao)); } } // Now make analysis handler aware of the weight names present. _weightNames.clear(); _weightNames.push_back(""); _defaultWeightIdx = 0; for ( string name : foundWeightNames ) _weightNames.push_back(name); // Then we create and initialize all analyses for ( string ananame : foundAnalyses ) addAnalysis(ananame); for (AnaHandle a : analyses() ) { MSG_TRACE("Initialising analysis: " << a->name()); if ( !a->info().reentrant() ) MSG_WARNING("Analysis " << a->name() << " has not been validated to have " << "a reentrant finalize method. The merged result is unpredictable."); try { // Allow projection registration in the init phase onwards a->_allowProjReg = true; a->init(); } catch (const Error& err) { cerr << "Error in " << a->name() << "::init method: " << err.what() << endl; exit(1); } MSG_TRACE("Done initialising analysis: " << a->name()); } // Now get all booked analysis objects. vector raos; for (AnaHandle a : analyses()) { for (const auto & ao : a->analysisObjects()) { raos.push_back(ao); } } // Collect global weights and xcoss sections and fix scaling for // all files. _eventCounter = CounterPtr(weightNames(), Counter("_EVTCOUNT")); _xs = Scatter1DPtr(weightNames(), Scatter1D("_XSEC")); for (size_t iW = 0; iW < numWeights(); iW++) { _eventCounter.get()->setActiveWeightIdx(iW); _xs.get()->setActiveWeightIdx(iW); YODA::Counter & sumw = *_eventCounter; YODA::Scatter1D & xsec = *_xs; vector xsecs; vector sows; for ( auto & aomap : allaos ) { auto xit = aomap.find(xsec.path()); if ( xit == aomap.end() ) xsecs.push_back(dynamic_pointer_cast(xit->second)); else xsecs.push_back(YODA::Scatter1DPtr()); xit = aomap.find(sumw.path()); if ( xit == aomap.end() ) sows.push_back(dynamic_pointer_cast(xit->second)); else sows.push_back(YODA::CounterPtr()); } double xs = 0.0, xserr = 0.0; for ( int i = 0, N = sows.size(); i < N; ++i ) { if ( !sows[i] || !xsecs[i] ) continue; double xseci = xsecs[i]->point(0).x(); double xsecerri = sqr(xsecs[i]->point(0).xErrAvg()); sumw += *sows[i]; double effnent = sows[i]->effNumEntries(); xs += (equiv? effnent: 1.0)*xseci; xserr += (equiv? sqr(effnent): 1.0)*xsecerri; } vector scales(sows.size(), 1.0); if ( equiv ) { xs /= sumw.effNumEntries(); xserr = sqrt(xserr)/sumw.effNumEntries(); } else { xserr = sqrt(xserr); for ( int i = 0, N = sows.size(); i < N; ++i ) scales[i] = (sumw.sumW()/sows[i]->sumW())* (xsecs[i]->point(0).x()/xs); } xsec.point(0) = Point1D(xs, xserr); // Go through alla analyses and add stuff to their analysis objects; for (AnaHandle a : analyses()) { for (const auto & ao : a->analysisObjects()) { ao.get()->setActiveWeightIdx(iW); YODA::AnalysisObjectPtr yao = ao.get()->activeYODAPtr(); for ( int i = 0, N = sows.size(); i < N; ++i ) { if ( !sows[i] || !xsecs[i] ) continue; auto range = allaos[i].equal_range(yao->path()); for ( auto aoit = range.first; aoit != range.second; ++aoit ) if ( !addaos(yao, aoit->second, scales[i]) ) MSG_WARNING("Cannot merge objects with path " << yao->path() <<" of type " << yao->annotation("Type") ); } ao.get()->unsetActiveWeight(); } } _eventCounter.get()->unsetActiveWeight(); _xs.get()->unsetActiveWeight(); } // Finally we just have to finalize all analyses, leaving to the // controlling program to write it out some yoda-file. finalize(); } void AnalysisHandler::readData(const string& filename) { try { /// @todo Use new YODA SFINAE to fill the smart ptr vector directly vector aos_raw; YODA::read(filename, aos_raw); for (YODA::AnalysisObject* aor : aos_raw) _preloads[aor->path()] = YODA::AnalysisObjectPtr(aor); } catch (...) { //< YODA::ReadError& throw UserError("Unexpected error in reading file: " + filename); } } vector AnalysisHandler::getRivetAOs() const { vector rtn; for (AnaHandle a : analyses()) { for (const auto & ao : a->analysisObjects()) { rtn.push_back(ao); } } rtn.push_back(_eventCounter); rtn.push_back(_xs); return rtn; } void AnalysisHandler::writeData(const string& filename) const { // This is where we store the OAs to be written. vector output; // First get all multiwight AOs vector raos = getRivetAOs(); output.reserve(raos.size()*2*numWeights()); // Fix the oredering so that default weight is written out first. vector order; if ( _defaultWeightIdx >= 0 && _defaultWeightIdx < numWeights() ) order.push_back(_defaultWeightIdx); for ( size_t i = 0; i < numWeights(); ++i ) if ( i != _defaultWeightIdx ) order.push_back(i); // Then we go through all finalized AOs one weight at a time for (size_t iW : order ) { for ( auto rao : raos ) { rao.get()->setActiveFinalWeightIdx(iW); if ( rao->path().find("/TMP/") != string::npos ) continue; output.push_back(rao.get()->activeYODAPtr()); } } // Finally the RAW objects. for (size_t iW : order ) { for ( auto rao : raos ) { rao.get()->setActiveWeightIdx(iW); output.push_back(rao.get()->activeYODAPtr()); } } try { YODA::write(filename, output.begin(), output.end()); } catch (...) { //< YODA::WriteError& throw UserError("Unexpected error in writing file: " + filename); } } string AnalysisHandler::runName() const { return _runname; } size_t AnalysisHandler::numEvents() const { return _eventCounter->numEntries(); } std::vector AnalysisHandler::analysisNames() const { std::vector rtn; for (AnaHandle a : analyses()) { rtn.push_back(a->name()); } return rtn; } AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector& analysisnames) { for (const string& aname : analysisnames) { //MSG_DEBUG("Adding analysis '" << aname << "'"); addAnalysis(aname); } return *this; } AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector& analysisnames) { for (const string& aname : analysisnames) { removeAnalysis(aname); } return *this; } - AnalysisHandler& AnalysisHandler::setCrossSection(double xs, double xserr) { + void AnalysisHandler::setCrossSection(pair xsec) { _xs = Scatter1DPtr(weightNames(), Scatter1D("_XSEC")); _eventCounter.get()->setActiveWeightIdx(_defaultWeightIdx); double nomwgt = sumW(); // The cross section of each weight variation is the nominal cross section // times the sumW(variation) / sumW(nominal). // This way the cross section will work correctly for (size_t iW = 0; iW < numWeights(); iW++) { _eventCounter.get()->setActiveWeightIdx(iW); double s = sumW() / nomwgt; _xs.get()->setActiveWeightIdx(iW); - _xs->addPoint(xs*s, xserr*s); + _xs->addPoint(xsec.first*s, xsec.second*s); } _eventCounter.get()->unsetActiveWeight(); _xs.get()->unsetActiveWeight(); - return *this; + return; } AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) { analysis->_analysishandler = this; // _analyses.insert(AnaHandle(analysis)); _analyses[analysis->name()] = AnaHandle(analysis); return *this; } PdgIdPair AnalysisHandler::beamIds() const { return Rivet::beamIds(beams()); } double AnalysisHandler::sqrtS() const { return Rivet::sqrtS(beams()); } void AnalysisHandler::setIgnoreBeams(bool ignore) { _ignoreBeams=ignore; } void AnalysisHandler::skipMultiWeights(bool ignore) { _skipWeights=ignore; } } diff --git a/src/Core/Run.cc b/src/Core/Run.cc --- a/src/Core/Run.cc +++ b/src/Core/Run.cc @@ -1,142 +1,143 @@ // -*- C++ -*- #include "Rivet/Run.hh" #include "Rivet/AnalysisHandler.hh" #include "Rivet/Math/MathUtils.hh" #include "Rivet/Tools/RivetPaths.hh" #include "zstr/zstr.hpp" #include #include using std::cout; using std::endl; namespace Rivet { Run::Run(AnalysisHandler& ah) : _ah(ah), _fileweight(1.0), _xs(NAN) { } Run::~Run() { } Run& Run::setCrossSection(const double xs) { _xs = xs; return *this; } Run& Run::setListAnalyses(const bool dolist) { _listAnalyses = dolist; return *this; } // Fill event and check for a bad read state bool Run::readEvent() { /// @todo Clear rather than new the GenEvent object per-event? _evt.reset(new GenEvent()); if(!HepMCUtils::readEvent(_hepmcReader, _evt)){ Log::getLog("Rivet.Run") << Log::DEBUG << "Read failed. End of file?" << endl; return false; } // Rescale event weights by file-level weight, if scaling is non-trivial if (!fuzzyEquals(_fileweight, 1.0)) { for (size_t i = 0; i < (size_t) _evt->weights().size(); ++i) { _evt->weights()[i] *= _fileweight; } } return true; } bool Run::openFile(const std::string& evtfile, double weight) { // Set current weight-scaling member _fileweight = weight; // In case makeReader fails. std::string errormessage; // Set up HepMC input reader objects if (evtfile == "-") { #ifdef HAVE_LIBZ _istr = make_shared(std::cin); _hepmcReader = HepMCUtils::makeReader(*_istr, &errormessage); #else _hepmcReader = HepMCUtils::makeReader(std::cin, &errormessage); #endif } else { if ( !fileexists(evtfile) ) throw Error("Event file '" + evtfile + "' not found"); #ifdef HAVE_LIBZ // NB. zstr auto-detects if file is deflated or plain-text _istr = make_shared(evtfile.c_str()); #else _istr = make_shared(evtfile.c_str()); #endif _hepmcReader = HepMCUtils::makeReader(*_istr, &errormessage); } if (_hepmcReader == nullptr) { Log::getLog("Rivet.Run") << Log::ERROR << "Read error on file '" << evtfile << "' " << errormessage << endl; return false; } return true; } bool Run::init(const std::string& evtfile, double weight) { if (!openFile(evtfile, weight)) return false; // Read first event to define run conditions bool ok = readEvent(); if (!ok) return false; if(HepMCUtils::particles(_evt).size() == 0){ Log::getLog("Rivet.Run") << Log::ERROR << "Empty first event." << endl; return false; } // Initialise AnalysisHandler with beam information from first event _ah.init(*_evt); // Set cross-section from command line if (!std::isnan(_xs)) { Log::getLog("Rivet.Run") << Log::DEBUG << "Setting user cross-section = " << _xs << " pb" << endl; - _ah.setCrossSection(_xs, 0.0); + + _ah.setCrossSection(make_pair(_xs, 0.0)); } // List the chosen & compatible analyses if requested if (_listAnalyses) { for (const std::string& ana : _ah.analysisNames()) { cout << ana << endl; } } return true; } bool Run::processEvent() { // Analyze event _ah.analyze(*_evt); return true; } bool Run::finalize() { _evt.reset(); _ah.finalize(); return true; } } diff --git a/src/Tools/RivetHepMC_2.cc b/src/Tools/RivetHepMC_2.cc --- a/src/Tools/RivetHepMC_2.cc +++ b/src/Tools/RivetHepMC_2.cc @@ -1,200 +1,201 @@ // -*- C++ -*- //#include #include "Rivet/Tools/Utils.hh" #include "Rivet/Tools/RivetHepMC.hh" #include "Rivet/Tools/Logging.hh" /*namespace { inline std::vector split(const std::string& input, const std::string& regex) { // passing -1 as the submatch index parameter performs splitting std::regex re(regex); std::sregex_token_iterator first{input.begin(), input.end(), re, -1}, last; return {first, last}; } }*/ namespace Rivet{ const Relatives Relatives::PARENTS = HepMC::parents; const Relatives Relatives::CHILDREN = HepMC::children; const Relatives Relatives::ANCESTORS = HepMC::ancestors; const Relatives Relatives::DESCENDANTS = HepMC::descendants; namespace HepMCUtils{ ConstGenParticlePtr getParticlePtr(const RivetHepMC::GenParticle & gp) { return &gp; } std::vector particles(ConstGenEventPtr ge){ std::vector result; for(GenEvent::particle_const_iterator pi = ge->particles_begin(); pi != ge->particles_end(); ++pi){ result.push_back(*pi); } return result; } std::vector particles(const GenEvent *ge){ std::vector result; for(GenEvent::particle_const_iterator pi = ge->particles_begin(); pi != ge->particles_end(); ++pi){ result.push_back(*pi); } return result; } std::vector vertices(ConstGenEventPtr ge){ std::vector result; for(GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi){ result.push_back(*vi); } return result; } std::vector vertices(const GenEvent *ge){ std::vector result; for(GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi){ result.push_back(*vi); } return result; } std::vector particles(ConstGenVertexPtr gv, const Relatives &relo){ std::vector result; /// @todo A particle_const_iterator on GenVertex would be nice... // Before HepMC 2.7.0 there were no GV::particles_const_iterators and constness consistency was all screwed up :-/ #if HEPMC_VERSION_CODE >= 2007000 for (HepMC::GenVertex::particle_iterator pi = gv->particles_begin(relo); pi != gv->particles_end(relo); ++pi) result.push_back(*pi); #else HepMC::GenVertex* gv2 = const_cast(gv); for (HepMC::GenVertex::particle_iterator pi = gv2->particles_begin(relo); pi != gv2->particles_end(relo); ++pi) result.push_back(const_cast(*pi)); #endif return result; } std::vector particles(ConstGenParticlePtr gp, const Relatives &relo){ ConstGenVertexPtr vtx; switch(relo){ case HepMC::parents: case HepMC::ancestors: vtx = gp->production_vertex(); break; case HepMC::children: case HepMC::descendants: vtx = gp->end_vertex(); break; default: throw std::runtime_error("Not implemented!"); break; } return particles(vtx, relo); } int uniqueId(ConstGenParticlePtr gp){ return gp->barcode(); } int particles_size(ConstGenEventPtr ge){ return ge->particles_size(); } int particles_size(const GenEvent *ge){ return ge->particles_size(); } std::pair beams(const GenEvent *ge){ return ge->beam_particles(); } std::shared_ptr makeReader(std::istream &istr, std::string *) { return make_shared(istr); } bool readEvent(std::shared_ptr io, std::shared_ptr evt){ if(io->rdstate() != 0) return false; if(!io->fill_next_event(evt.get())) return false; return true; } // This functions could be filled with code doing the same stuff as // in the HepMC3 version of This file. void strip(GenEvent &, const set &) {} vector weightNames(const GenEvent & ge) { /// reroute the print output to a std::stringstream and process /// The iteration is done over a map in hepmc2 so this is safe vector ret; /// Obtaining weight names using regex probably neater, but regex /// is not defined in GCC4.8, which is currently used by Lxplus. /// Attempt an alternative solution based on stringstreams: std::stringstream stream; ge.weights().print(stream); std::string pair; // placeholder for subtsring matches while (std::getline(stream, pair, ' ')) { if ( pair.size() < 2 ) continue; pair.erase(pair.begin()); // removes the "(" on the LHS pair.pop_back(); // removes the ")" on the RHS if (pair.empty()) continue; std::stringstream spair(pair); vector temp; while (std::getline(spair, pair, ',')) { temp.push_back(std::move(pair)); } if (temp.size() == 2) { // store the default weight based on weight names if (temp[0] == "Weight" || temp[0] == "0" || temp[0] == "Default") { ret.push_back(""); } else ret.push_back(temp[0]); } } /// Possible future solution based on regex /*std::ostringstream stream; ge.weights().print(stream); // Super lame, I know string str = stream.str(); std::regex re("(([^()]+))"); // Regex for stuff enclosed by parentheses () for (std::sregex_iterator i = std::sregex_iterator(str.begin(), str.end(), re); i != std::sregex_iterator(); ++i ) { std::smatch m = *i; vector temp = ::split(m.str(), "[,]"); if (temp.size() ==2) { // store the default weight based on weight names if (temp[0] == "Weight" || temp[0] == "0" || temp[0] == "Default") { ret.push_back(""); } else ret.push_back(temp[0]); } }*/ return ret; } - double crossSection(const GenEvent & ge) { - return ge.cross_section()->cross_section(); + pair crossSection(const GenEvent & ge) { + return make_pair(ge.cross_section()->cross_section(), + ge.cross_section()->cross_section_error()); } std::valarray weights(const GenEvent & ge) { const size_t W = ge.weights().size(); std::valarray wts(W); for (unsigned int iw = 0; iw < W; ++iw) wts[iw] = ge.weights()[iw]; return wts; } } } diff --git a/src/Tools/RivetHepMC_3.cc b/src/Tools/RivetHepMC_3.cc --- a/src/Tools/RivetHepMC_3.cc +++ b/src/Tools/RivetHepMC_3.cc @@ -1,159 +1,179 @@ // -*- C++ -*- #include "Rivet/Tools/RivetHepMC.hh" #include "Rivet/Tools/ReaderCompressedAscii.hh" #include "HepMC3/ReaderAscii.h" #include "HepMC3/ReaderAsciiHepMC2.h" +#include "HepMC3/GenCrossSection.h" +#include namespace Rivet{ namespace HepMCUtils{ ConstGenParticlePtr getParticlePtr(const RivetHepMC::GenParticle & gp) { return gp.shared_from_this(); } std::vector particles(ConstGenEventPtr ge){ return ge->particles(); } std::vector particles(const GenEvent *ge){ assert(ge); return ge->particles(); } std::vector vertices(ConstGenEventPtr ge){ return ge->vertices(); } std::vector vertices(const GenEvent *ge){ assert(ge); return ge->vertices(); } std::vector particles(ConstGenVertexPtr gv, const Relatives &relo){ return relo(gv); } std::vector particles(ConstGenParticlePtr gp, const Relatives &relo){ return relo(gp); } int particles_size(ConstGenEventPtr ge){ return particles(ge).size(); } int particles_size(const GenEvent *ge){ return particles(ge).size(); } int uniqueId(ConstGenParticlePtr gp){ return gp->id(); } std::pair beams(const GenEvent *ge) { std::vector beamlist = ge->beams(); if ( beamlist.size() < 2 ) { std::cerr << "CANNOT FIND ANY BEAMS!" << std::endl; return std::pair(); } return std::make_pair(beamlist[0], beamlist[1]); } bool readEvent(std::shared_ptr io, std::shared_ptr evt){ return io->read_event(*evt) && !io->failed(); } shared_ptr makeReader(std::istream & istr, std::string * errm) { shared_ptr ret; // First scan forward and check if there is some hint as to what // kind of file we are looking att. int ntry = 10; std::string header; int filetype = -1; while ( ntry ) { std::getline(istr, header); if ( header.empty() ) continue; if ( header.substr(0, 34) == "HepMC::Asciiv3-START_EVENT_LISTING" ) { filetype = 3; break; } if ( header.substr(0, 44) == "HepMC::CompressedAsciiv3-START_EVENT_LISTING" ) { filetype = 4; break; } if ( header.substr(0, 38) == "HepMC::IO_GenEvent-START_EVENT_LISTING" ) { filetype = 2; break; } --ntry; } if ( filetype == 3 ) ret = make_shared(istr); else if ( filetype == 4 ) ret = make_shared(istr); else ret = make_shared(istr); if ( filetype == 0 && errm ) *errm += "Could not determine file type. Assuming HepMC2 file. "; // Check that everything was ok. if ( ret->failed() ) { if ( errm ) *errm = "Problems reading from HepMC file."; ret = shared_ptr(); } return ret; } void strip(GenEvent & ge, const set & stripid) { // std::cerr << "Stripping event " << ge.event_number() << std::endl; vector allparticles = ge.particles(); for ( auto & p : allparticles ) { if ( !p->production_vertex() || !p->end_vertex() || stripid.count(p->pid()) == 0 || p->production_vertex()->id() == 0 ) continue; // std::cout << "Removing particle " << p->id() // << " (" << p->pid() << ")" << std::endl; HepMC3::GenVertexPtr vp = p->production_vertex(); HepMC3::GenVertexPtr ve = p->end_vertex(); if ( !vp || !ve ) continue; if ( vp == ve ) continue; // Check if the vertices would leave particles with the sam // production as decay vertex - we don't want that. if ( ( vp->particles_out().size() == 1 && vp->particles_out()[0] == p ) || ( ve->particles_in().size() == 1 && ve->particles_in()[0] == p ) ) { bool loop = false; for ( auto pi : vp->particles_in() ) for ( auto po : ve->particles_out() ) if ( pi == po ) loop = true; if ( loop ) continue; } if ( vp->particles_in().size() == 1 && ( vp->particles_in()[0]->pid() > 21 && vp->particles_in()[0]->pid() < 30 ) ) continue; vp->remove_particle_out(p); ve->remove_particle_in(p); if ( ve->particles_in().empty() ) { auto prem = ve->particles_out(); for ( auto po : prem ) vp->add_particle_out(po); ge.remove_vertex(ve); } else if ( vp->particles_out().empty() ) { auto prem = vp->particles_in(); for ( auto pi : prem ) ve->add_particle_in(pi); ge.remove_vertex(vp); } ge.remove_particle(p); } } + pair crossSection(const GenEvent & ge) { + // Work-around since access functions are not const. + HepMC3::GenCrossSection xs = *ge.cross_section(); + return make_pair(xs.xsec(), xs.xsec_err()); + } + + vector weightNames(const GenEvent & ge) { + try { + return ge.weight_names(""); + } catch (HepMC3::WeightError & w) { + return vector(); + } + } + + std::valarray weights(const GenEvent & ge) { + return std::valarray(&ge.weights()[0], ge.weights().size()); + + } } } diff --git a/test/testApi.cc b/test/testApi.cc --- a/test/testApi.cc +++ b/test/testApi.cc @@ -1,34 +1,35 @@ #include "Rivet/AnalysisHandler.hh" #include "HepMC/GenEvent.h" #include "Rivet/Tools/RivetHepMC.hh" +#include using namespace std; int main(int argc, char* argv[]) { assert(argc > 1); Rivet::AnalysisHandler ah; Rivet::Log::setLevel("Rivet", Rivet::Log::DEBUG); // Specify the analyses to be used ah.addAnalysis("EXAMPLE"); ah.addAnalyses({{ "MC_JETS", "EXAMPLE_CUTS", "EXAMPLE_SMEAR" }}); std::ifstream file("testApi.hepmc"); shared_ptr reader = Rivet::HepMCUtils::makeReader(file); std::shared_ptr evt = make_shared(); double sum_of_weights = 0.0; while ( Rivet::HepMCUtils::readEvent(reader, evt) ) { // Analyse current event ah.analyze(*evt); sum_of_weights += evt->weights()[0]; } file.close(); - ah.setCrossSection(1.0, 0.1); + ah.setCrossSection(make_pair(1.0, 0.1)); ah.finalize(); ah.writeData("out.yoda"); return 0; } diff --git a/test/testNaN.cc b/test/testNaN.cc --- a/test/testNaN.cc +++ b/test/testNaN.cc @@ -1,74 +1,75 @@ #include "Rivet/AnalysisHandler.hh" #include "Rivet/Analysis.hh" #include "Rivet/Tools/RivetYODA.hh" #include #include #include +#include using namespace std; class Test : public Rivet::Analysis { public: Test() : Analysis("Test") {} void init() { book(_h_test, "test", 50, 66.0, 116.0); } void analyze(const Rivet::Event & e) { cout << "Normal fill" << endl; _h_test->fill(90., 1.); cout << "Underflow fill" << endl; _h_test->fill(30.,1.); cout << "Overflow fill" << endl; _h_test->fill(130.,1.); cout << "Inf fill" << endl; try { _h_test->fill(numeric_limits::infinity(), 1.); } catch (YODA::RangeError & e) { cerr << e.what() << '\n'; if ( string(e.what()) != string("X is Inf") ) throw; } cout << "NaN fill" << endl; try { _h_test->fill(numeric_limits::quiet_NaN(), 1.); } catch (YODA::RangeError & e) { cerr << e.what() << '\n'; if ( string(e.what()) != string("X is NaN") ) throw; } } private: Rivet::Histo1DPtr _h_test; }; DECLARE_RIVET_PLUGIN(Test); int main(int argc, char* argv[]) { assert(argc > 1); Rivet::AnalysisHandler rivet; rivet.addAnalysis("Test"); std::ifstream file("testApi.hepmc"); shared_ptr reader = Rivet::HepMCUtils::makeReader(file); std::shared_ptr evt = make_shared(); double sum_of_weights = 0.0; while ( Rivet::HepMCUtils::readEvent(reader, evt) ) { // Analyse current event rivet.analyze(*evt); sum_of_weights += evt->weights()[0]; } file.close(); - rivet.setCrossSection(1.0, 0.1); + rivet.setCrossSection(make_pair(1.0, 0.1)); rivet.finalize(); rivet.writeData("NaN.aida"); return 0; }