diff --git a/include/Rivet/Tools/RivetYODA.hh b/include/Rivet/Tools/RivetYODA.hh --- a/include/Rivet/Tools/RivetYODA.hh +++ b/include/Rivet/Tools/RivetYODA.hh @@ -1,609 +1,697 @@ #ifndef RIVET_RIVETYODA_HH #define RIVET_RIVETYODA_HH #include "Rivet/Config/RivetCommon.hh" #include "YODA/AnalysisObject.h" #include "YODA/Counter.h" #include "YODA/Histo1D.h" #include "YODA/Histo2D.h" #include "YODA/Profile1D.h" #include "YODA/Profile2D.h" #include "YODA/Scatter1D.h" #include "YODA/Scatter2D.h" #include "YODA/Scatter3D.h" #include #include namespace YODA { typedef std::shared_ptr AnalysisObjectPtr; typedef std::shared_ptr CounterPtr; typedef std::shared_ptr Histo1DPtr; typedef std::shared_ptr Histo2DPtr; typedef std::shared_ptr Profile1DPtr; typedef std::shared_ptr Profile2DPtr; typedef std::shared_ptr Scatter1DPtr; typedef std::shared_ptr Scatter2DPtr; typedef std::shared_ptr Scatter3DPtr; } namespace Rivet { class AnalysisObjectWrapper { public: virtual ~AnalysisObjectWrapper() {} virtual YODA::AnalysisObject* operator->() = 0; virtual YODA::AnalysisObject* operator->() const = 0; virtual const YODA::AnalysisObject & operator*() const = 0; /// @todo Rename to setActive(idx) virtual void setActiveWeightIdx(unsigned int iWeight) = 0; /// @todo Set active object for finalize virtual void setActiveFinalWeightIdx(unsigned int iWeight) = 0; virtual void unsetActiveWeight() = 0; bool operator ==(const AnalysisObjectWrapper& p) { return (this == &p); } bool operator !=(const AnalysisObjectWrapper& p) { return (this != &p); } protected: /// @todo do we need this? // virtual void reset() = 0; }; /// @todo /// implement scatter1dptr and scatter2dptr here /// these need to be multi-weighted eventually. /* class Scatter1DPtr : public AnalysisObjectPtr { public: Scatter1DPtr() : _persistent() { } Scatter1DPtr(size_t len_of_weightvec, const YODA::Scatter1D& p) { for (size_t m = 0; m < len_of_weightvec; ++m) _persistent.push_back(make_shared(p)); } bool operator!() const { return !_persistent; } explicit operator bool() const { return bool(_persistent); } YODA::Scatter1D* operator->() { return _persistent.get(); } YODA::Scatter1D* operator->() const { return _persistent.get(); } YODA::Scatter1D & operator*() { return *_persistent; } const YODA::Scatter1D & operator*() const { return *_persistent; } protected: vector _persistent; }; class Scatter2DPtr : public AnalysisObjectPtr { public: Scatter2DPtr(size_t len_of_weightvec, const YODA::Scatter2D& p) { for (size_t m = 0; m < len_of_weightvec; ++m) _persistent.push_back(make_shared(p)); } Scatter2DPtr() : _persistent() { } bool operator!() { return !_persistent; } explicit operator bool() { return bool(_persistent); } YODA::Scatter2D* operator->() { return _persistent.get(); } YODA::Scatter2D* operator->() const { return _persistent.get(); } YODA::Scatter2D & operator*() { return *_persistent; } const YODA::Scatter2D & operator*() const { return *_persistent; } protected: vector _persistent; }; class Scatter3DPtr : public AnalysisObjectPtr { public: Scatter3DPtr(size_t len_of_weightvec, const YODA::Scatter3D& p) { for (size_t m = 0; m < len_of_weightvec; ++m) _persistent.push_back(make_shared(p)); } Scatter3DPtr() : _persistent() { } bool operator!() { return !_persistent; } explicit operator bool() { return bool(_persistent); } YODA::Scatter3D* operator->() { return _persistent.get(); } YODA::Scatter3D* operator->() const { return _persistent.get(); } YODA::Scatter3D & operator*() { return *_persistent; } const YODA::Scatter3D & operator*() const { return *_persistent; } protected: vector _persistent; }; */ class MultiweightAOWrapper : public AnalysisObjectWrapper { public: using Inner = YODA::AnalysisObject; virtual void newSubEvent() = 0; virtual void pushToPersistent(const vector >& weight) = 0; virtual void pushToFinal() = 0; virtual YODA::AnalysisObjectPtr activeYODAPtr() const = 0; }; using Weight = double; template using Fill = pair; template using Fills = multiset>; // TODO TODO TODO // need to override the old fill method too! // otherwise we can't intercept existing fill calls in analysis code // TODO TODO TODO /// Wrappers for analysis objects to store all fills unaggregated, until collapsed template class TupleWrapper; template<> class TupleWrapper : public YODA::Counter { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Counter & h) : YODA::Counter(h) {} // todo: do we need to deal with users using fractions directly? void fill( double weight=1.0, double fraction=1.0 ) { fills_.insert( {YODA::Counter::FillType(),weight} ); } void reset() { fills_.clear(); } const Fills & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Histo1D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Histo1D & h) : YODA::Histo1D(h) {} // todo: do we need to deal with users using fractions directly? void fill( double x, double weight=1.0, double fraction=1.0 ) { if ( std::isnan(x) ) throw YODA::RangeError("X is NaN"); fills_.insert( { x , weight } ); } void reset() { fills_.clear(); } const Fills & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Profile1D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Profile1D & h) : YODA::Profile1D(h) {} // todo: do we need to deal with users using fractions directly? void fill( double x, double y, double weight=1.0, double fraction=1.0 ) { if ( std::isnan(x) ) throw YODA::RangeError("X is NaN"); if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN"); fills_.insert( { YODA::Profile1D::FillType{x,y}, weight } ); } void reset() { fills_.clear(); } const Fills & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Histo2D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Histo2D & h) : YODA::Histo2D(h) {} // todo: do we need to deal with users using fractions directly? void fill( double x, double y, double weight=1.0, double fraction=1.0 ) { if ( std::isnan(x) ) throw YODA::RangeError("X is NaN"); if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN"); fills_.insert( { YODA::Histo2D::FillType{x,y}, weight } ); } void reset() { fills_.clear(); } const Fills & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Profile2D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Profile2D & h) : YODA::Profile2D(h) {} // todo: do we need to deal with users using fractions directly? void fill( double x, double y, double z, double weight=1.0, double fraction=1.0 ) { if ( std::isnan(x) ) throw YODA::RangeError("X is NaN"); if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN"); if ( std::isnan(z) ) throw YODA::RangeError("Z is NaN"); fills_.insert( { YODA::Profile2D::FillType{x,y,z}, weight } ); } void reset() { fills_.clear(); } const Fills & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Scatter1D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Scatter1D & h) : YODA::Scatter1D(h) {} }; template<> class TupleWrapper : public YODA::Scatter2D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Scatter2D & h) : YODA::Scatter2D(h) {} }; template<> class TupleWrapper : public YODA::Scatter3D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Scatter3D & h) : YODA::Scatter3D(h) {} }; template class Wrapper : public MultiweightAOWrapper { friend class Analysis; public: using Inner = T; /* @todo * some things are not really well-defined here * for instance: fill() in the finalize() method and integral() in * the analyze() method. */ Wrapper() = default; Wrapper(const vector& weightnames, const T & p); ~Wrapper(); typename T::Ptr active() const; /* @todo this probably need to loop over all? */ bool operator!() const { return !_active; } // Don't use active() here, assert will catch explicit operator bool() const { return static_cast(_active); } // Don't use active() here, assert will catch T * operator->() { return active().get(); } T * operator->() const { return active().get(); } T & operator*() { return *active(); } const T & operator*() const { return *active(); } /* @todo * these need to be re-thought out. void reset() { active()->reset(); } */ /* @todo * these probably need to loop over all? * do we even want to provide equality? */ /* @todo * how about no. friend bool operator==(Wrapper a, Wrapper b){ if (a._persistent.size() != b._persistent.size()) return false; for (size_t i = 0; i < a._persistent.size(); i++) { if (a._persistent.at(i) != b._persistent.at(i)) { return false; } } return true; } friend bool operator!=(Wrapper a, Wrapper b){ return !(a == b); } friend bool operator<(Wrapper a, Wrapper b){ if (a._persistent.size() >= b._persistent.size()) return false; for (size_t i = 0; i < a._persistent.size(); i++) { if (*(a._persistent.at(i)) >= *(b._persistent.at(i))) { return false; } } return true; } */ private: void setActiveWeightIdx(unsigned int iWeight) { _active = _persistent.at(iWeight); } void setActiveFinalWeightIdx(unsigned int iWeight) { _active = _final.at(iWeight); } /* this is for dev only---we shouldn't need this in real runs. */ void unsetActiveWeight() { _active.reset(); } void newSubEvent(); virtual YODA::AnalysisObjectPtr activeYODAPtr() const { return _active; } const vector & persistent() const { return _persistent; } const vector & final() const { return _final; } /* to be implemented for each type */ void pushToPersistent(const vector >& weight); void pushToFinal(); /* M of these, one for each weight */ vector _persistent; /* This is the copy of _persistent that will be passed to finalize(). */ vector _final; /* N of these, one for each event in evgroup */ vector::Ptr> _evgroup; typename T::Ptr _active; // do we need implicit cast? // operator typename T::Ptr () { // return _active; // } friend class AnalysisHandler; }; /// We need our own shared_ptr class, so we can dispatch -> and * /// all the way down to the inner YODA analysis objects /// /// TODO: provide remaining functionality that shared_ptr has (not needed right now) /// template class rivet_shared_ptr { public: typedef T value_type; rivet_shared_ptr() = default; rivet_shared_ptr(decltype(nullptr)) : _p(nullptr) {} /// Convenience constructor, pass through to the Wrapper constructor rivet_shared_ptr(const vector& weightNames, const typename T::Inner & p) : _p( make_shared(weightNames, p) ) {} template rivet_shared_ptr(const shared_ptr & p) : _p(p) {} template rivet_shared_ptr(const rivet_shared_ptr & p) : _p(p.get()) {} // Goes right through to the active YODA object's members T & operator->() { return *_p; } const T & operator->() const { return *_p; } // The active YODA object typename T::Inner & operator*() { return **_p; } const typename T::Inner & operator*() const { return **_p; } bool operator!() const { return !_p || !(*_p); } explicit operator bool() const { return _p && bool(*_p); } template bool operator==(const rivet_shared_ptr & other) const { return _p == other._p; } template bool operator!=(const rivet_shared_ptr & other) const { return _p != other._p; } template bool operator<(const rivet_shared_ptr & other) const { return _p < other._p; } template bool operator>(const rivet_shared_ptr & other) const { return _p > other._p; } template bool operator<=(const rivet_shared_ptr & other) const { return _p <= other._p; } template bool operator>=(const rivet_shared_ptr & other) const { return _p >= other._p; } shared_ptr get() const { return _p; } private: shared_ptr _p; }; // every object listed here needs a virtual fill method in YODA, // otherwise the Tuple fakery won't work. using MultiweightAOPtr = rivet_shared_ptr; using Histo1DPtr = rivet_shared_ptr>; using Histo2DPtr = rivet_shared_ptr>; using Profile1DPtr = rivet_shared_ptr>; using Profile2DPtr = rivet_shared_ptr>; using CounterPtr = rivet_shared_ptr>; using Scatter1DPtr = rivet_shared_ptr>; using Scatter2DPtr = rivet_shared_ptr>; using Scatter3DPtr = rivet_shared_ptr>; using YODA::Counter; using YODA::Histo1D; using YODA::HistoBin1D; using YODA::Histo2D; using YODA::HistoBin2D; using YODA::Profile1D; using YODA::ProfileBin1D; using YODA::Profile2D; using YODA::ProfileBin2D; using YODA::Scatter1D; using YODA::Point1D; using YODA::Scatter2D; using YODA::Point2D; using YODA::Scatter3D; using YODA::Point3D; /// Function to get a map of all the refdata in a paper with the /// given @a papername. map getRefData(const string& papername); /// @todo Also provide a Scatter3D getRefData() version? /// Get the file system path to the reference file for this paper. string getDatafilePath(const string& papername); /// Traits class to access the type of the AnalysisObject in the /// reference files. template struct ReferenceTraits {}; template<> struct ReferenceTraits { typedef Counter RefT; }; template<> struct ReferenceTraits { typedef Scatter1D RefT; }; template<> struct ReferenceTraits { typedef Scatter2D RefT; }; template<> struct ReferenceTraits { typedef Scatter2D RefT; }; template<> struct ReferenceTraits { typedef Scatter2D RefT; }; template<> struct ReferenceTraits { typedef Scatter3D RefT; }; template<> struct ReferenceTraits { typedef Scatter3D RefT; }; template<> struct ReferenceTraits { typedef Scatter3D RefT; }; /// If @a dst and @a src both are of same subclass T, copy the /// contents of @a src into @a dst and return true. Otherwise return /// false. template inline bool aocopy(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst) { shared_ptr tsrc = dynamic_pointer_cast(src); if ( !tsrc ) return false; shared_ptr tdst = dynamic_pointer_cast(dst); if ( !tdst ) return false; *tdst = *tsrc; return true; } /// If @a dst and @a src both are of same subclass T, add the /// contents of @a src into @a dst and return true. Otherwise return /// false. template inline bool aoadd(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale) { shared_ptr tsrc = dynamic_pointer_cast(src); if ( !tsrc ) return false; shared_ptr tdst = dynamic_pointer_cast(dst); if ( !tdst ) return false; tsrc->scaleW(scale); *tdst += *tsrc; return true; } /// If @a dst is the same subclass as @a src, copy the contents of @a /// src into @a dst and return true. Otherwise return false. bool copyao(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst); /// If @a dst is the same subclass as @a src, scale the contents of /// @a src with @a scale and add it to @a dst and return true. Otherwise /// return false. bool addaos(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale); /// Check if two analysis objects have the same binning or, if not /// binned, are in other ways compatible. template inline bool bookingCompatible(TPtr a, TPtr b) { return a->sameBinning(*b); } inline bool bookingCompatible(CounterPtr, CounterPtr) { return true; } inline bool bookingCompatible(Scatter1DPtr a, Scatter1DPtr b) { return a->numPoints() == b->numPoints(); } inline bool bookingCompatible(Scatter2DPtr a, Scatter2DPtr b) { return a->numPoints() == b->numPoints(); } inline bool bookingCompatible(Scatter3DPtr a, Scatter3DPtr b) { return a->numPoints() == b->numPoints(); } -inline bool bookingCompatible(YODA::CounterPtr, YODA::CounterPtr) { + inline bool bookingCompatible(YODA::CounterPtr, YODA::CounterPtr) { return true; } inline bool bookingCompatible(YODA::Scatter1DPtr a, YODA::Scatter1DPtr b) { return a->numPoints() == b->numPoints(); } inline bool bookingCompatible(YODA::Scatter2DPtr a, YODA::Scatter2DPtr b) { return a->numPoints() == b->numPoints(); } inline bool bookingCompatible(YODA::Scatter3DPtr a, YODA::Scatter3DPtr b) { return a->numPoints() == b->numPoints(); } + /// class representing a YODA path with all its components. + class AOPath { + + public: + + /// Constructor + AOPath(string fullpath); + + /// The full path. + string path() const { return _path; } + + /// The analysis name. + string analysis() const { return _analysis; } + + /// The analysis name with options. + string analysisWithOptions() const { return _analysis + _optionstring; } + + /// The base name of the analysis object. + string name() const { return _name; } + + /// The weight name. + string weight() const { return _weight; } + + /// Is This a RAW (filling) object? + bool isRaw() const { return _raw; } + + // Is This a temporary (filling) object? + bool isTmp() const { return _tmp; } + + /// Is This a reference object? + bool isRef() const { return _ref; } + + /// The string describing the options passed to the analysis. + string optionString() const { return _optionstring; } + + /// Are there options passed to the analysis? + bool hasOptions() const { return !_options.empty(); } + + /// Don't pass This optionto the analysis + void removeOption(string opt) { _options.erase(opt); fixOptionString(); } + + /// Pass this option to the analysis. + void setOption(string opt, string val) { _options[opt] = val; fixOptionString();} + + /// Was This option passed to the analyisi. + bool hasOption(string opt) const { return _options.find(opt) != _options.end(); } + + /// Get the value of this option. + string getOption(string opt) const { + auto it = _options.find(opt); + if ( it != _options.end() ) return it->second; + return ""; + } + + /// Reset the option string after changes; + void fixOptionString(); + + /// Creat a full path (and set) for this. + string mkPath() const; + string setPath() { return _path = mkPath(); } + + /// Print out information + void debug() const; + + /// Make this class ordered. + bool operator<(const AOPath & other) const { + return _path < other._path; + } + + /// Check if path is valid. + bool valid() const { return _valid; }; + bool operator!() const { return !valid(); } + + private: + + bool _valid; + string _path; + string _analysis; + string _optionstring; + string _name; + string _weight; + bool _raw; + bool _tmp; + bool _ref; + map _options; + + }; + } #endif diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc --- a/src/Core/AnalysisHandler.cc +++ b/src/Core/AnalysisHandler.cc @@ -1,751 +1,780 @@ // -*- 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 #include using std::cout; using std::cerr; 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 { AnalysisHandler::AnalysisHandler(const string& runname) : _runname(runname), // *** LEIF *** temporarily removed this // _eventcounter("/_EVTCOUNT"), // _xs(NAN), _xserr(NAN), _initialised(false), _ignoreBeams(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 (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); } // 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"); } // *** LEIF *** Reinstated this. void AnalysisHandler::setWeightNames(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 std::ostringstream stream; ge.weights().print(stream); // Super lame, I know string str = stream.str(); std::regex re("(([^()]+))"); // Regex for stuff enclosed by parentheses () size_t idx = 0; 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) { MSG_DEBUG("Name of weight #" << _weightNames.size() << ": " << temp[0]); // store the default weight based on weight names if (temp[0] == "Weight" || temp[0] == "0" || temp[0] == "Default") { MSG_DEBUG(_weightNames.size() << " is being used as the nominal."); _weightNames.push_back(""); _defaultWeightIdx = idx; } else _weightNames.push_back(temp[0]); idx++; } } } 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); // Weights /// @todo Drop this / just report first weight when we support multiweight events // *** LEIF *** temporarily removed this // _eventcounter.fill(event.weight()); // MSG_DEBUG("Event #" << _eventcounter.numEntries() << " weight = " << event.weight()); // *** LEIF *** temporarily removed this // // Cross-section // #if defined ENABLE_HEPMC_3 // if (ge.cross_section()) { // //@todo HepMC3::GenCrossSection methods aren't const accessible :( // RivetHepMC::GenCrossSection gcs = *(event.genEvent()->cross_section()); // _xs = gcs.xsec(); // _xserr = gcs.xsec_err(); // } // #elif defined HEPMC_HAS_CROSS_SECTION // if (ge.cross_section()) { // _xs = ge.cross_section()->cross_section(); // _xserr = ge.cross_section()->cross_section_error(); // } // #endif // 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."); } _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::finalize() { if (!_initialised) return; MSG_INFO("Finalising analyses"); // First push all analyses' objects to persistent and final MSG_TRACE("AnalysisHandler::finalize(): Pushing analysis objects to persistent."); _eventCounter.get()->pushToPersistent(_subEventWeights); _eventCounter.get()->pushToFinal(); _xs.get()->pushToFinal(); for (const AnaHandle& a : analyses()) { for (auto ao : a->analysisObjects()) { ao.get()->pushToPersistent(_subEventWeights); 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" : "")); if ( _dumping ) return; // Print out MCnet boilerplate 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 & , - const vector & , bool ) {} - // void AnalysisHandler::mergeYodas(const vector & aofiles, - // const vector & delopts, bool equiv) { - // vector< vector > aosv; - // vector xsecs; - // vector xsecerrs; - // vector sows; - // set ananames; - // _eventCounter->reset(); + void AnalysisHandler::mergeYodas(const vector & aofiles, + const vector & delopts, bool equiv) { - // // First scan all files and extract analysis objects and add the - // // corresponding analyses.. - // for (const string& file : aofiles ) { - // Scatter1DPtr xsec; - // CounterPtr sow; + // Convenient typedef; + typedef multimap AOMap; - // // For each file make sure that cross section and sum-of-weights - // // objects are present and stor all RAW ones in a vector; - // vector aos; - // try { - // /// @todo Use new YODA SFINAE to fill the smart ptr vector directly - // vector aos_raw; - // YODA::read(file, aos_raw); - // for (YODA::AnalysisObject* aor : aos_raw) { - // YODA::AnalysisObjectPtr ao(aor); - // if ( ao->path().substr(0, 5) != "/RAW/" ) continue; - // ao->setPath(ao->path().substr(4)); - // if ( ao->path() == "/_XSEC" ) - // xsec = dynamic_pointer_cast(ao); - // else if ( ao->path() == "/_EVTCOUNT" ) - // sow = dynamic_pointer_cast(ao); - // else { - // stripOptions(ao, delopts); - // string ananame = split(ao->path(), "/")[0]; - // if ( ananames.insert(ananame).second ) addAnalysis(ananame); - // aos.push_back(ao); - // } - // } - // if ( !xsec || !sow ) { - // MSG_ERROR( "Error in AnalysisHandler::mergeYodas: The file " << file - // << " did not contain weights and cross section info."); - // exit(1); - // } - // xsecs.push_back(xsec->point(0).x()); - // sows.push_back(sow); - // xsecerrs.push_back(sqr(xsec->point(0).xErrAvg())); - // _eventCounter->operator+=(*sow); //< HAHAHAHAHA! - // sows.push_back(sow); - // aosv.push_back(aos); - // } catch (...) { //< YODA::ReadError& - // throw UserError("Unexpected error in reading file: " + file); - // } - // } + // Store all found weights here. + set foundWeightNames; - // // Now calculate the scale to be applied for all bins in a file - // // and get the common cross section and sum of weights. - // _xs = _xserr = 0.0; - // for ( int i = 0, N = sows.size(); i < N; ++i ) { - // double effnent = sows[i]->effNumEntries(); - // _xs += (equiv? effnent: 1.0)*xsecs[i]; - // _xserr += (equiv? sqr(effnent): 1.0)*xsecerrs[i]; - // } + // Stor all found analyses. + set foundAnalyses; - // vector scales(sows.size(), 1.0); - // if ( equiv ) { - // _xs /= _eventCounter.effNumEntries(); - // _xserr = sqrt(_xserr)/_eventCounter.effNumEntries(); - // } else { - // _xserr = sqrt(_xserr); - // for ( int i = 0, N = sows.size(); i < N; ++i ) - // scales[i] = (_eventCounter.sumW()/sows[i]->sumW())*(xsecs[i]/_xs); - // } + // Store all analysis objects here. + vector allaos; - // // Initialize the analyses allowing them to book analysis objects. - // for (AnaHandle a : _analyses) { - // MSG_DEBUG("Initialising analysis: " << a->name()); - // if ( !a->info().reentrant() ) - // MSG_WARNING("Analysis " << a->name() << " has not been validated to have " - // << "a reentrant finalize method. The result is unpredictable."); - // try { - // // Allow projection registration in the init phase onwards - // a->_allowProjReg = true; - // cerr << "sqrtS " << sqrtS() << endl; - // 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()); - // } - // _initialised = true; - // // Get a list of all anaysis objects to handle. - // map current; - // for ( auto ao : getData(false, true) ) current[ao->path()] = ao; - // // Go through all objects to be merged and add them to current - // // after appropriate scaling. - // for ( int i = 0, N = aosv.size(); i < N; ++i) - // for ( auto ao : aosv[i] ) { - // if ( ao->path() == "/_XSEC" || ao->path() == "_EVTCOUNT" ) continue; - // auto aoit = current.find(ao->path()); - // if ( aoit == current.end() ) { - // MSG_WARNING("" << ao->path() << " was not properly booked."); - // continue; - // } - // if ( !addaos(aoit->second, ao, scales[i]) ) - // MSG_WARNING("Cannot merge objects with path " << ao->path() - // <<" of type " << ao->annotation("Type") ); - // } - // // Now we can simply finalize() the analysis, leaving the - // // controlling program to write it out some yoda-file. - // finalize(); - // } + // 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(""); + 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; } // *** LEIF *** thinks This is not needed anymore // vector AnalysisHandler::getYodaAOs(bool includeorphans, // bool includetmps, // bool usefinalized) const { // vector rtn; // if (usefinalized) // rtn = _finalizedAOs; // else { // for (auto rao : getRivetAOs()) { // // need to set the index // // before we can search the PATH // rao.get()->setActiveWeightIdx(_defaultWeightIdx); // // Exclude paths from final write-out if they contain a "TMP" layer (i.e. matching "/TMP/") // if (!includetmps && rao->path().find("/TMP/") != string::npos) // continue; // for (size_t iW = 0; iW < numWeights(); iW++) { // rao.get()->setActiveWeightIdx(iW); // rtn.push_back(rao.get()->activeYODAPtr()); // } // } // } // // Sort histograms alphanumerically by path before write-out // sort(rtn.begin(), rtn.end(), // [](YODA::AnalysisObjectPtr a, YODA::AnalysisObjectPtr b) { // return a->path() < b->path(); // } // ); // return rtn; // } // vector AnalysisHandler::getData(bool includeorphans, // bool includetmps, // bool usefinalized) const { // return getYodaAOs(includeorphans, includetmps, usefinalized); // } 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); + output.reserve(raos.size()*2*numWeights()); // Then we go through all finalized AOs one weight at a time for (size_t iW = 0; iW < numWeights(); iW++) { 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 = 0; iW < numWeights(); iW++) { for ( auto rao : raos ) { rao.get()->setActiveFinalWeightIdx(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) { _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); } _eventCounter.get()->unsetActiveWeight(); _xs.get()->unsetActiveWeight(); return *this; } // *** LEIF *** temporarily removed this // bool AnalysisHandler::hasCrossSection() const { // return (!std::isnan(crossSection())); // } 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; } } diff --git a/src/Tools/RivetYODA.cc b/src/Tools/RivetYODA.cc --- a/src/Tools/RivetYODA.cc +++ b/src/Tools/RivetYODA.cc @@ -1,428 +1,504 @@ #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Tools/RivetYODA.hh" #include "Rivet/Tools/RivetPaths.hh" #include "YODA/ReaderYODA.h" #include "YODA/ReaderAIDA.h" // use execinfo for backtrace if available #include "DummyConfig.hh" #ifdef HAVE_EXECINFO_H #include #endif +#include +#include + using namespace std; namespace Rivet { template Wrapper::~Wrapper() {} template Wrapper::Wrapper(const vector& weightNames, const T & p) { for (const string& weightname : weightNames) { _persistent.push_back(make_shared(p)); _final.push_back(make_shared(p)); auto obj = _persistent.back(); auto final = _final.back(); if (weightname != "") { obj->setPath("/RAW" + obj->path() + "[" + weightname + "]"); final->setPath(final->path() + "[" + weightname + "]"); } } } template typename T::Ptr Wrapper::active() const { if ( !_active ) { #ifdef HAVE_BACKTRACE void * buffer[4]; backtrace(buffer, 4); backtrace_symbols_fd(buffer, 4 , 1); #endif assert(false && "No active pointer set. Was this object booked in init()?"); } return _active; } template void Wrapper::newSubEvent() { typename TupleWrapper::Ptr tmp = make_shared>(_persistent[0]->clone()); tmp->reset(); _evgroup.push_back( tmp ); _active = _evgroup.back(); assert(_active); } string getDatafilePath(const string& papername) { /// Try to find YODA otherwise fall back to try AIDA const string path1 = findAnalysisRefFile(papername + ".yoda"); if (!path1.empty()) return path1; const string path2 = findAnalysisRefFile(papername + ".aida"); if (!path2.empty()) return path2; throw Rivet::Error("Couldn't find ref data file '" + papername + ".yoda" + " in data path, '" + getRivetDataPath() + "', or '.'"); } map getRefData(const string& papername) { const string datafile = getDatafilePath(papername); // Make an appropriate data file reader and read the data objects /// @todo Remove AIDA support some day... YODA::Reader& reader = (datafile.find(".yoda") != string::npos) ? \ YODA::ReaderYODA::create() : YODA::ReaderAIDA::create(); vector aovec; reader.read(datafile, aovec); // Return value, to be populated map rtn; for ( YODA::AnalysisObject* ao : aovec ) { YODA::AnalysisObjectPtr refdata(ao); if (!refdata) continue; const string plotpath = refdata->path(); // Split path at "/" and only return the last field, i.e. the histogram ID const size_t slashpos = plotpath.rfind("/"); const string plotname = (slashpos+1 < plotpath.size()) ? plotpath.substr(slashpos+1) : ""; rtn[plotname] = refdata; } return rtn; } } namespace { using Rivet::Fill; using Rivet::Fills; using Rivet::TupleWrapper; template double get_window_size(const typename T::Ptr & histo, typename T::BinType x) { // the bin index we fall in const auto binidx = histo->binIndexAt(x); // gaps, overflow, underflow don't contribute if ( binidx == -1 ) return 0; const auto & b = histo->bin(binidx); // if we don't have a valid neighbouring bin, // we use infinite width typename T::Bin b1(-1.0/0.0, 1.0/0.0); // points in the top half compare to the upper neighbour if ( x > b.xMid() ) { size_t nextidx = binidx + 1; if ( nextidx < histo->bins().size() ) b1 = histo->bin(nextidx); } else { // compare to the lower neighbour int nextidx = binidx - 1; if ( nextidx >= 0 ) b1 = histo->bin(nextidx); } // the factor 2 is arbitrary, could poss. be smaller return min( b.width(), b1.width() ) / 2.0; } template typename T::BinType fillT2binT(typename T::FillType a) { return a; } template <> YODA::Profile1D::BinType fillT2binT(YODA::Profile1D::FillType a) { return get<0>(a); } template <> YODA::Profile2D::BinType fillT2binT(YODA::Profile2D::FillType a) { return YODA::Profile2D::BinType{ get<0>(a), get<1>(a) }; } template void commit(vector & persistent, const vector< vector> > & tuple, const vector> & weights ) { // TODO check if all the xs are in the same bin anyway! // Then no windowing needed assert(persistent.size() == weights[0].size()); for ( const auto & x : tuple ) { double maxwindow = 0.0; for ( const auto & xi : x ) { // TODO check for NOFILL here // persistent[0] has the same binning as all the other persistent objects double window = get_window_size(persistent[0], fillT2binT(xi.first)); if ( window > maxwindow ) maxwindow = window; } const double wsize = maxwindow; // all windows have same size set edgeset; // bin edges need to be in here! for ( const auto & xi : x ) { edgeset.insert(fillT2binT(xi.first) - wsize); edgeset.insert(fillT2binT(xi.first) + wsize); } vector< std::tuple,double> > hfill; double sumf = 0.0; auto edgit = edgeset.begin(); double ehi = *edgit; while ( ++edgit != edgeset.end() ) { double elo = ehi; ehi = *edgit; valarray sumw(0.0, persistent.size()); // need m copies of this bool gap = true; // Check for gaps between the sub-windows. for ( size_t i = 0; i < x.size(); ++i ) { // check equals comparisons here! if ( fillT2binT(x[i].first) + wsize >= ehi && fillT2binT(x[i].first) - wsize <= elo ) { sumw += x[i].second * weights[i]; gap = false; } } if ( gap ) continue; hfill.push_back( make_tuple( (ehi + elo)/2.0, sumw, (ehi - elo) ) ); sumf += ehi - elo; } for ( auto f : hfill ) for ( size_t m = 0; m < persistent.size(); ++m ) persistent[m]->fill( get<0>(f), get<1>(f)[m], get<2>(f)/sumf ); // Note the scaling to one single fill } } template<> void commit(vector & persistent, const vector< vector> > & tuple, const vector> & weights) {} template<> void commit(vector & persistent, const vector< vector> > & tuple, const vector> & weights) {} template double distance(T a, T b) { return abs(a - b); } template <> double distance >(tuple a, tuple b) { return Rivet::sqr(get<0>(a) - get<0>(b)) + Rivet::sqr(get<1>(a) - get<1>(b)); } } namespace Rivet { bool copyao(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst) { for (const std::string& a : src->annotations()) dst->setAnnotation(a, src->annotation(a)); if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; return false; } bool addaos(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale) { if ( aoadd(dst,src,scale) ) return true; if ( aoadd(dst,src,scale) ) return true; if ( aoadd(dst,src,scale) ) return true; if ( aoadd(dst,src,scale) ) return true; if ( aoadd(dst,src,scale) ) return true; return false; } } namespace { /// fills is a vector of sub-event with an ordered set of x-values of /// the fills in each sub-event. NOFILL should be an "impossible" /// value for this histogram. Returns a vector of sub-events with /// an ordered vector of fills (including NOFILLs) for each sub-event. template vector< vector > > match_fills(const vector::Ptr> & evgroup, const Fill & NOFILL) { vector< vector > > matched; // First just copy subevents into vectors and find the longest vector. unsigned int maxfill = 0; // length of biggest vector int imax = 0; // index position of biggest vector for ( const auto & it : evgroup ) { const auto & subev = it->fills(); if ( subev.size() > maxfill ) { maxfill = subev.size(); imax = matched.size(); } matched.push_back(vector >(subev.begin(), subev.end())); } // Now, go through all subevents with missing fills. const vector> & full = matched[imax]; // the longest one for ( auto & subev : matched ) { if ( subev.size() == maxfill ) continue; // Add NOFILLs to the end; while ( subev.size() < maxfill ) subev.push_back(NOFILL); // Iterate from the back and shift all fill values backwards by // swapping with NOFILLs so that they better match the full // subevent. for ( int i = maxfill - 1; i >= 0; --i ) { if ( subev[i] == NOFILL ) continue; size_t j = i; while ( j + 1 < maxfill && subev[j + 1] == NOFILL && distance(fillT2binT(subev[j].first), fillT2binT(full[j].first)) > distance(fillT2binT(subev[j].first), fillT2binT(full[j + 1].first)) ) { swap(subev[j], subev[j + 1]); ++j; } } } // transpose vector>> result(maxfill,vector>(matched.size())); for (size_t i = 0; i < matched.size(); ++i) for (size_t j = 0; j < maxfill; ++j) result.at(j).at(i) = matched.at(i).at(j); return result; } } namespace Rivet { template void Wrapper::pushToPersistent(const vector >& weight) { assert( _evgroup.size() == weight.size() ); // have we had subevents at all? const bool have_subevents = _evgroup.size() > 1; if ( ! have_subevents ) { // simple replay of all tuple entries // each recorded fill is inserted into all persistent weightname histos for ( size_t m = 0; m < _persistent.size(); ++m ) for ( const auto & f : _evgroup[0]->fills() ) _persistent[m]->fill( f.first, f.second * weight[0][m] ); } else { // outer index is subevent, inner index is jets in the event vector>> linedUpXs = match_fills(_evgroup, {typename T::FillType(), 0.0}); commit( _persistent, linedUpXs, weight ); } _evgroup.clear(); _active.reset(); } template void Wrapper::pushToFinal() { for ( size_t m = 0; m < _persistent.size(); ++m ) { copyao(_persistent.at(m), _final.at(m)); } } template <> void Wrapper::pushToPersistent(const vector >& weight) { for ( size_t m = 0; m < _persistent.size(); ++m ) { for ( size_t n = 0; n < _evgroup.size(); ++n ) { for ( const auto & f : _evgroup[n]->fills() ) { _persistent[m]->fill( f.second * weight[n][m] ); } } } _evgroup.clear(); _active.reset(); } template <> void Wrapper::pushToPersistent(const vector >& weight) { _evgroup.clear(); _active.reset(); } template <> void Wrapper::pushToPersistent(const vector >& weight) { _evgroup.clear(); _active.reset(); } template <> void Wrapper::pushToPersistent(const vector >& weight) { _evgroup.clear(); _active.reset(); } // explicitly instantiate all wrappers template class Wrapper; template class Wrapper; template class Wrapper; template class Wrapper; template class Wrapper; template class Wrapper; template class Wrapper; template class Wrapper; + AOPath::AOPath(string fullpath) { + // First checck if This is a global system path + _path = fullpath; + std::regex resys("^(/RAW)?/([^\\[/]+)(\\[(.+)\\])?$"); + smatch m; + _valid = regex_search(fullpath, m, resys); + if ( _valid ) { + _raw = (m[1] == "/RAW"); + _name = m[2]; + _weight = m[4]; + return; + } + // If not, assume it is a normal analysis path. + std::regex repath("^(/RAW)?(/REF)?/([^/:]+)(:[^/]+)?(/TMP)?/([^\\[]+)(\\[(.+)\\])?"); + _valid = regex_search(fullpath, m, repath); + if ( !_valid ) return; + _raw = (m[1] == "/RAW"); + _ref = (m[2] == "/REF"); + _analysis = m[3]; + _optionstring = m[4]; + _tmp = (m[5] == "/TMP"); + _name = m[6]; + _weight = m[8]; + std::regex reopt(":([^=]+)=([^:]+)"); + string s = _optionstring; + while ( regex_search(s, m, reopt) ) { + _options[m[1]] = m[2]; + s = m.suffix(); + } + } + + void AOPath::fixOptionString() { + ostringstream oss; + for ( auto optval : _options ) + oss << ":" << optval.first << "=" << optval.second; + _optionstring = oss.str(); + } + + string AOPath::mkPath() const { + ostringstream oss; + if ( isRaw() ) oss << "/RAW"; + else if ( isRef() ) oss << "/REF"; + if ( _analysis != "" ) oss << "/" << analysis(); + for ( auto optval : _options ) + oss << ":" << optval.first << "=" << optval.second; + if ( isTmp() ) oss << "/TMP"; + oss << "/" << name(); + if ( weight() != "" ) + oss << "[" << weight() << "]"; + return oss.str(); + } + + void AOPath::debug() const { + cout << "Full path: " << _path << endl; + if ( !_valid ) { + cout << "This is not a valid analysis object path" << endl << endl; + return; + } + cout << "Check path: " << mkPath() << endl; + cout << "Analysis: " << _analysis << endl; + cout << "Name: " << _name << endl; + cout << "Weight: " << _weight << endl; + cout << "Properties: "; + if ( _raw ) cout << "raw "; + if ( _tmp ) cout << "tmp "; + if ( _ref ) cout << "ref "; + cout << endl; + cout << "Options: "; + for ( auto opt : _options ) + cout << opt.first << "->" << opt.second << " "; + cout << endl << endl; + } + }