diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh --- a/include/Rivet/AnalysisHandler.hh +++ b/include/Rivet/AnalysisHandler.hh @@ -1,318 +1,316 @@ // -*- 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; // Needed to make smart pointers compare equivalent in the STL set struct CmpAnaHandle { bool operator() (const AnaHandle& a, const AnaHandle& b) const { return a.get() < b.get(); } }; /// 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 //@{ /// Get the name of this run. string runName() const { return _runname; } /// 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 { return _eventcounter.numEntries(); } /// @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(); } /// @brief Compatibility alias for sumOfWeights /// /// @deprecated Prefer sumW double sumOfWeights() const { return sumW(); } /// @brief Set the sum of weights /// /// This is useful if Rivet is steered externally and /// the analyses are run for a sub-contribution of the events /// (but of course have to be normalised to the total sum of weights) /// /// @todo What about the sumW2 term? That needs to be set coherently. Need a /// new version, with all three N,sumW,sumW2 numbers (or a counter) /// supplied. /// /// @deprecated Weight sums are no longer tracked this way... void setSumOfWeights(const double& sum) { //_sumOfWeights = sum; throw Error("Can't set sum of weights independently, since it's now tracked by a Counter. " "Please contact the Rivet authors if you need this."); } /// Is cross-section information required by at least one child analysis? /// @deprecated Should no-longer be an issue: does any generator not write the cross-section? bool needCrossSection() const; /// Whether the handler knows about a cross-section. /// @deprecated Should no-longer be an issue: does any generator not write the cross-section? bool hasCrossSection() const; /// Get the cross-section known to the handler. double crossSection() const { return _xs; } /// Set the cross-section for the process being generated. AnalysisHandler& setCrossSection(double xs, double xserr=0); /// 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); //@} /// @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::set& analyses() const { return _analyses; } /// Get a registered analysis by name. const AnaHandle analysis(const std::string& analysisname) const; /// 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 //@{ /// Add a vector of analysis objects to the current state. void addData(const std::vector& aos); /// Read analysis plots into the histo collection (via addData) from the named file. void readData(const std::string& filename); /// Get all analyses' plots as a vector of analysis objects. std::vector getData(bool includeorphans = false, - bool includetmps = false) const; - - /// Get all finalized analyses' plots as a vector of analysis objects. - std::vector getFinalizedData(bool includeorphans = false) const; - + bool includetmps = false, + bool usefinalized = true) const; + /// 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 anslysis. 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(AnalysisObjectPtr ao, const vector & delopts) const; //@} private: /// The collection of Analysis objects to be used. set _analyses; /// A vector of pre-loaded object which do not have a valid /// Analysis plugged in. vector _orphanedPreloads; /// A vector containing copies of analysis objects after /// finalize() has been run. vector _finalizedAOs; /// @name Run properties //@{ /// Run name std::string _runname; /// Event counter Counter _eventcounter; /// Cross-section known to AH double _xs, _xserr; /// 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; /// 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/src/Core/Analysis.cc b/src/Core/Analysis.cc --- a/src/Core/Analysis.cc +++ b/src/Core/Analysis.cc @@ -1,964 +1,964 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Analysis.hh" #include "Rivet/AnalysisHandler.hh" #include "Rivet/AnalysisInfo.hh" #include "Rivet/Tools/BeamConstraint.hh" #include "Rivet/Projections/ImpactParameterProjection.hh" #include "Rivet/Projections/GeneratedPercentileProjection.hh" #include "Rivet/Projections/UserCentEstimate.hh" #include "Rivet/Projections/CentralityProjection.hh" namespace Rivet { Analysis::Analysis(const string& name) : _crossSection(-1.0), _gotCrossSection(false), _analysishandler(NULL) { ProjectionApplier::_allowProjReg = false; _defaultname = name; unique_ptr ai = AnalysisInfo::make(name); assert(ai); _info = move(ai); assert(_info); } double Analysis::sqrtS() const { return handler().sqrtS(); } const ParticlePair& Analysis::beams() const { return handler().beams(); } const PdgIdPair Analysis::beamIds() const { return handler().beamIds(); } const string Analysis::histoDir() const { /// @todo Cache in a member variable string _histoDir; if (_histoDir.empty()) { _histoDir = "/" + name(); if (handler().runName().length() > 0) { _histoDir = "/" + handler().runName() + _histoDir; } replace_all(_histoDir, "//", "/"); //< iterates until none } return _histoDir; } const string Analysis::histoPath(const string& hname) const { const string path = histoDir() + "/" + hname; return path; } const string Analysis::histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return histoDir() + "/" + mkAxisCode(datasetId, xAxisId, yAxisId); } const string Analysis::mkAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { stringstream axisCode; axisCode << "d"; if (datasetId < 10) axisCode << 0; axisCode << datasetId; axisCode << "-x"; if (xAxisId < 10) axisCode << 0; axisCode << xAxisId; axisCode << "-y"; if (yAxisId < 10) axisCode << 0; axisCode << yAxisId; return axisCode.str(); } Log& Analysis::getLog() const { string logname = "Rivet.Analysis." + name(); return Log::getLog(logname); } /////////////////////////////////////////// size_t Analysis::numEvents() const { return handler().numEvents(); } double Analysis::sumW() const { return handler().sumW(); } double Analysis::sumW2() const { return handler().sumW2(); } /////////////////////////////////////////// bool Analysis::isCompatible(const ParticlePair& beams) const { return isCompatible(beams.first.pid(), beams.second.pid(), beams.first.energy(), beams.second.energy()); } bool Analysis::isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const { PdgIdPair beams(beam1, beam2); pair energies(e1, e2); return isCompatible(beams, energies); } // bool Analysis::beamIDsCompatible(const PdgIdPair& beams) const { // bool beamIdsOk = false; // for (const PdgIdPair& bp : requiredBeams()) { // if (compatible(beams, bp)) { // beamIdsOk = true; // break; // } // } // return beamIdsOk; // } // /// Check that the energies are compatible (within 1% or 1 GeV, whichever is larger, for a bit of UI forgiveness) // bool Analysis::beamEnergiesCompatible(const pair& energies) const { // /// @todo Use some sort of standard ordering to improve comparisons, esp. when the two beams are different particles // bool beamEnergiesOk = requiredEnergies().size() > 0 ? false : true; // typedef pair DoublePair; // for (const DoublePair& ep : requiredEnergies()) { // if ((fuzzyEquals(ep.first, energies.first, 0.01) && fuzzyEquals(ep.second, energies.second, 0.01)) || // (fuzzyEquals(ep.first, energies.second, 0.01) && fuzzyEquals(ep.second, energies.first, 0.01)) || // (abs(ep.first - energies.first) < 1*GeV && abs(ep.second - energies.second) < 1*GeV) || // (abs(ep.first - energies.second) < 1*GeV && abs(ep.second - energies.first) < 1*GeV)) { // beamEnergiesOk = true; // break; // } // } // return beamEnergiesOk; // } // bool Analysis::beamsCompatible(const PdgIdPair& beams, const pair& energies) const { bool Analysis::isCompatible(const PdgIdPair& beams, const pair& energies) const { // First check the beam IDs bool beamIdsOk = false; for (const PdgIdPair& bp : requiredBeams()) { if (compatible(beams, bp)) { beamIdsOk = true; break; } } if (!beamIdsOk) return false; // Next check that the energies are compatible (within 1% or 1 GeV, whichever is larger, for a bit of UI forgiveness) /// @todo Use some sort of standard ordering to improve comparisons, esp. when the two beams are different particles bool beamEnergiesOk = requiredEnergies().size() > 0 ? false : true; typedef pair DoublePair; for (const DoublePair& ep : requiredEnergies()) { if ((fuzzyEquals(ep.first, energies.first, 0.01) && fuzzyEquals(ep.second, energies.second, 0.01)) || (fuzzyEquals(ep.first, energies.second, 0.01) && fuzzyEquals(ep.second, energies.first, 0.01)) || (abs(ep.first - energies.first) < 1*GeV && abs(ep.second - energies.second) < 1*GeV) || (abs(ep.first - energies.second) < 1*GeV && abs(ep.second - energies.first) < 1*GeV)) { beamEnergiesOk = true; break; } } return beamEnergiesOk; } /////////////////////////////////////////// Analysis& Analysis::setCrossSection(double xs) { _crossSection = xs; _gotCrossSection = true; return *this; } double Analysis::crossSection() const { if (!_gotCrossSection || std::isnan(_crossSection)) { string errMsg = "You did not set the cross section for the analysis " + name(); throw Error(errMsg); } return _crossSection; } double Analysis::crossSectionPerEvent() const { const double sumW = sumOfWeights(); assert(sumW != 0.0); return _crossSection / sumW; } //////////////////////////////////////////////////////////// // Histogramming void Analysis::_cacheRefData() const { if (_refdata.empty()) { MSG_TRACE("Getting refdata cache for paper " << name()); _refdata = getRefData(getRefDataName()); } } vector Analysis::getAllData(bool includeorphans) const{ - return handler().getData(includeorphans); + return handler().getData(includeorphans, false, false); } CounterPtr Analysis::bookCounter(const string& cname, const string& title) { return addOrGetCompatAO(make_shared(histoPath(cname), title)); } CounterPtr Analysis::bookCounter(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title) { return bookCounter(mkAxisCode(datasetId, xAxisId, yAxisId), title); } Histo1DPtr Analysis::bookHisto1D(const string& hname, size_t nbins, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { Histo1DPtr hist = make_shared(nbins, lower, upper, histoPath(hname), title); hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); return addOrGetCompatAO(hist); } Histo1DPtr Analysis::bookHisto1D(const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { Histo1DPtr hist = make_shared(binedges, histoPath(hname), title); hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); return addOrGetCompatAO(hist); } Histo1DPtr Analysis::bookHisto1D(const string& hname, const initializer_list& binedges, const string& title, const string& xtitle, const string& ytitle) { return bookHisto1D(hname, vector{binedges}, title, xtitle, ytitle); } Histo1DPtr Analysis::bookHisto1D(const string& hname, const Scatter2D& refscatter, const string& title, const string& xtitle, const string& ytitle) { Histo1DPtr hist = make_shared(refscatter, histoPath(hname)); hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); return addOrGetCompatAO(hist); } Histo1DPtr Analysis::bookHisto1D(const string& hname, const string& title, const string& xtitle, const string& ytitle) { const Scatter2D& refdata = refData(hname); return bookHisto1D(hname, refdata, title, xtitle, ytitle); } Histo1DPtr Analysis::bookHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title, const string& xtitle, const string& ytitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookHisto1D(axisCode, title, xtitle, ytitle); } /// @todo Add booking methods which take a path, titles and *a reference Scatter from which to book* ///////////////// Histo2DPtr Analysis::bookHisto2D(const string& hname, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Histo2DPtr hist = make_shared(nxbins, xlower, xupper, nybins, ylower, yupper, path, title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); hist->setAnnotation("ZLabel", ztitle); return addOrGetCompatAO(hist); } Histo2DPtr Analysis::bookHisto2D(const string& hname, const vector& xbinedges, const vector& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Histo2DPtr hist = make_shared(xbinedges, ybinedges, path, title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); hist->setAnnotation("ZLabel", ztitle); return addOrGetCompatAO(hist); } Histo2DPtr Analysis::bookHisto2D(const string& hname, const initializer_list& xbinedges, const initializer_list& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { return bookHisto2D(hname, vector{xbinedges}, vector{ybinedges}, title, xtitle, ytitle, ztitle); } Histo2DPtr Analysis::bookHisto2D(const string& hname, const Scatter3D& refscatter, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Histo2DPtr hist( new Histo2D(refscatter, path) ); if (hist->hasAnnotation("IsRef")) hist->rmAnnotation("IsRef"); hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); hist->setAnnotation("ZLabel", ztitle); return addOrGetCompatAO(hist); } Histo2DPtr Analysis::bookHisto2D(const string& hname, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const Scatter3D& refdata = refData(hname); return bookHisto2D(hname, refdata, title, xtitle, ytitle, ztitle); } Histo2DPtr Analysis::bookHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookHisto2D(axisCode, title, xtitle, ytitle, ztitle); } ///////////////// Profile1DPtr Analysis::bookProfile1D(const string& hname, size_t nbins, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Profile1DPtr prof = make_shared(nbins, lower, upper, path, title); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); return addOrGetCompatAO(prof); } Profile1DPtr Analysis::bookProfile1D(const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Profile1DPtr prof = make_shared(binedges, path, title); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); return addOrGetCompatAO(prof); } Profile1DPtr Analysis::bookProfile1D(const string& hname, const initializer_list& binedges, const string& title, const string& xtitle, const string& ytitle) { return bookProfile1D(hname, vector{binedges}, title, xtitle, ytitle); } Profile1DPtr Analysis::bookProfile1D(const string& hname, const Scatter2D& refscatter, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Profile1DPtr prof = make_shared(refscatter, path); if (prof->hasAnnotation("IsRef")) prof->rmAnnotation("IsRef"); prof->setTitle(title); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); return addOrGetCompatAO(prof); } Profile1DPtr Analysis::bookProfile1D(const string& hname, const string& title, const string& xtitle, const string& ytitle) { const Scatter2D& refdata = refData(hname); return bookProfile1D(hname, refdata, title, xtitle, ytitle); } Profile1DPtr Analysis::bookProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title, const string& xtitle, const string& ytitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookProfile1D(axisCode, title, xtitle, ytitle); } /////////////////// Profile2DPtr Analysis::bookProfile2D(const string& hname, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Profile2DPtr prof = make_shared(nxbins, xlower, xupper, nybins, ylower, yupper, path, title); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); prof->setAnnotation("ZLabel", ztitle); return addOrGetCompatAO(prof); } Profile2DPtr Analysis::bookProfile2D(const string& hname, const vector& xbinedges, const vector& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Profile2DPtr prof = make_shared(xbinedges, ybinedges, path, title); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); prof->setAnnotation("ZLabel", ztitle); return addOrGetCompatAO(prof); } Profile2DPtr Analysis::bookProfile2D(const string& hname, const initializer_list& xbinedges, const initializer_list& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { return bookProfile2D(hname, vector{xbinedges}, vector{ybinedges}, title, xtitle, ytitle, ztitle); } Profile2DPtr Analysis::bookProfile2D(const string& hname, const Scatter3D& refscatter, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Profile2DPtr prof( new Profile2D(refscatter, path) ); if (prof->hasAnnotation("IsRef")) prof->rmAnnotation("IsRef"); prof->setTitle(title); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); prof->setAnnotation("ZLabel", ztitle); return addOrGetCompatAO(prof); } Profile2DPtr Analysis::bookProfile2D(const string& hname, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const Scatter3D& refdata = refData(hname); return bookProfile2D(hname, refdata, title, xtitle, ytitle, ztitle); } Profile2DPtr Analysis::bookProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookProfile2D(axisCode, title, xtitle, ytitle, ztitle); } ///////////////// Scatter2DPtr Analysis::bookScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, bool copy_pts, const string& title, const string& xtitle, const string& ytitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookScatter2D(axisCode, copy_pts, title, xtitle, ytitle); } Scatter2DPtr Analysis::bookScatter2D(const string& hname, bool copy_pts, const string& title, const string& xtitle, const string& ytitle) { Scatter2DPtr s; const string path = histoPath(hname); if (copy_pts) { const Scatter2D& refdata = refData(hname); s = make_shared(refdata, path); for (Point2D& p : s->points()) p.setY(0, 0); } else { s = make_shared(path); } if (s->hasAnnotation("IsRef")) s->rmAnnotation("IsRef"); s->setTitle(title); s->setAnnotation("XLabel", xtitle); s->setAnnotation("YLabel", ytitle); return addOrGetCompatAO(s); } Scatter2DPtr Analysis::bookScatter2D(const string& hname, size_t npts, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Scatter2DPtr s = make_shared(path); const double binwidth = (upper-lower)/npts; for (size_t pt = 0; pt < npts; ++pt) { const double bincentre = lower + (pt + 0.5) * binwidth; s->addPoint(bincentre, 0, binwidth/2.0, 0); } s->setTitle(title); s->setAnnotation("XLabel", xtitle); s->setAnnotation("YLabel", ytitle); return addOrGetCompatAO(s); } Scatter2DPtr Analysis::bookScatter2D(const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Scatter2DPtr s = make_shared(path); for (size_t pt = 0; pt < binedges.size()-1; ++pt) { const double bincentre = (binedges[pt] + binedges[pt+1]) / 2.0; const double binwidth = binedges[pt+1] - binedges[pt]; s->addPoint(bincentre, 0, binwidth/2.0, 0); } s->setTitle(title); s->setAnnotation("XLabel", xtitle); s->setAnnotation("YLabel", ytitle); return addOrGetCompatAO(s); } Scatter2DPtr Analysis::bookScatter2D(const Scatter2DPtr scPtr, const std::string& path, const std::string& title, const std::string& xtitle, const std::string& ytitle ) { Scatter2DPtr s = make_shared(*scPtr); s->setPath(path); s->setTitle(title); s->setAnnotation("XLabel",xtitle); s->setAnnotation("YLabel",ytitle); return addOrGetCompatAO(s); } ///////////////////// void Analysis::divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const { const string path = s->path(); *s = *c1 / *c2; s->setPath(path); } void Analysis::divide(const Counter& c1, const Counter& c2, Scatter1DPtr s) const { const string path = s->path(); *s = c1 / c2; s->setPath(path); } void Analysis::divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const { const string path = s->path(); *s = *h1 / *h2; s->setPath(path); } void Analysis::divide(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const { const string path = s->path(); *s = h1 / h2; s->setPath(path); } void Analysis::divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const { const string path = s->path(); *s = *p1 / *p2; s->setPath(path); } void Analysis::divide(const Profile1D& p1, const Profile1D& p2, Scatter2DPtr s) const { const string path = s->path(); *s = p1 / p2; s->setPath(path); } void Analysis::divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const { const string path = s->path(); *s = *h1 / *h2; s->setPath(path); } void Analysis::divide(const Histo2D& h1, const Histo2D& h2, Scatter3DPtr s) const { const string path = s->path(); *s = h1 / h2; s->setPath(path); } void Analysis::divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const { const string path = s->path(); *s = *p1 / *p2; s->setPath(path); } void Analysis::divide(const Profile2D& p1, const Profile2D& p2, Scatter3DPtr s) const { const string path = s->path(); *s = p1 / p2; s->setPath(path); } /// @todo Counter and Histo2D efficiencies and asymms void Analysis::efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const { const string path = s->path(); *s = YODA::efficiency(*h1, *h2); s->setPath(path); } void Analysis::efficiency(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const { const string path = s->path(); *s = YODA::efficiency(h1, h2); s->setPath(path); } void Analysis::asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const { const string path = s->path(); *s = YODA::asymm(*h1, *h2); s->setPath(path); } void Analysis::asymm(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const { const string path = s->path(); *s = YODA::asymm(h1, h2); s->setPath(path); } void Analysis::scale(CounterPtr cnt, double factor) { if (!cnt) { MSG_WARNING("Failed to scale counter=NULL in analysis " << name() << " (scale=" << factor << ")"); return; } if (std::isnan(factor) || std::isinf(factor)) { MSG_WARNING("Failed to scale counter=" << cnt->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")"); factor = 0; } MSG_TRACE("Scaling counter " << cnt->path() << " by factor " << factor); try { cnt->scaleW(factor); } catch (YODA::Exception& we) { MSG_WARNING("Could not scale counter " << cnt->path()); return; } } void Analysis::normalize(Histo1DPtr histo, double norm, bool includeoverflows) { if (!histo) { MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")"); return; } MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm); try { const double hint = histo->integral(includeoverflows); if (hint == 0) MSG_WARNING("Skipping histo with null area " << histo->path()); else histo->normalize(norm, includeoverflows); } catch (YODA::Exception& we) { MSG_WARNING("Could not normalize histo " << histo->path()); return; } } void Analysis::scale(Histo1DPtr histo, double factor) { if (!histo) { MSG_WARNING("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")"); return; } if (std::isnan(factor) || std::isinf(factor)) { MSG_WARNING("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")"); factor = 0; } MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor); try { histo->scaleW(factor); } catch (YODA::Exception& we) { MSG_WARNING("Could not scale histo " << histo->path()); return; } } void Analysis::normalize(Histo2DPtr histo, double norm, bool includeoverflows) { if (!histo) { MSG_ERROR("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")"); return; } MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm); try { const double hint = histo->integral(includeoverflows); if (hint == 0) MSG_WARNING("Skipping histo with null area " << histo->path()); else histo->normalize(norm, includeoverflows); } catch (YODA::Exception& we) { MSG_WARNING("Could not normalize histo " << histo->path()); return; } } void Analysis::scale(Histo2DPtr histo, double factor) { if (!histo) { MSG_ERROR("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")"); return; } if (std::isnan(factor) || std::isinf(factor)) { MSG_ERROR("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")"); factor = 0; } MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor); try { histo->scaleW(factor); } catch (YODA::Exception& we) { MSG_WARNING("Could not scale histo " << histo->path()); return; } } void Analysis::integrate(Histo1DPtr h, Scatter2DPtr s) const { // preserve the path info const string path = s->path(); *s = toIntegralHisto(*h); s->setPath(path); } void Analysis::integrate(const Histo1D& h, Scatter2DPtr s) const { // preserve the path info const string path = s->path(); *s = toIntegralHisto(h); s->setPath(path); } /// @todo 2D versions of integrate... defined how, exactly?!? ////////////////////////////////// void Analysis::addAnalysisObject(AnalysisObjectPtr ao) { _analysisobjects.push_back(ao); } void Analysis::removeAnalysisObject(const string& path) { for (vector::iterator it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) { if ((*it)->path() == path) { _analysisobjects.erase(it); break; } } } void Analysis::removeAnalysisObject(AnalysisObjectPtr ao) { for (vector::iterator it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) { if (*it == ao) { _analysisobjects.erase(it); break; } } } const CentralityProjection & Analysis::declareCentrality(const SingleValueProjection &proj, string calAnaName, string calHistName, const string projName, bool increasing) { CentralityProjection cproj; // Select the centrality variable from option. Use REF as default. // Other selections are "GEN", "IMP" and "USR" (USR only in HEPMC 3). string sel = getOption("cent","REF"); set done; if ( sel == "REF" ) { Scatter2DPtr refscat; auto refmap = getRefData(calAnaName); if ( refmap.find(calHistName) != refmap.end() ) refscat = dynamic_pointer_cast(refmap.find(calHistName)->second); if ( !refscat ) { MSG_WARNING("No reference calibration histogram for " << "CentralityProjection " << projName << " found " << "(requested histogram " << calHistName << " in " << calAnaName << ")"); } else { MSG_INFO("Found calibration histogram " << sel << " " << refscat->path()); cproj.add(PercentileProjection(proj, refscat, increasing), sel); } } else if ( sel == "GEN" ) { Histo1DPtr genhist; string histpath = "/" + calAnaName + "/" + calHistName; - for ( AnalysisObjectPtr ao : handler().getData(true) ) { + for ( AnalysisObjectPtr ao : handler().getData(true, false, false) ) { if ( ao->path() == histpath ) genhist = dynamic_pointer_cast(ao); } if ( !genhist || genhist->numEntries() <= 1 ) { MSG_WARNING("No generated calibration histogram for " << "CentralityProjection " << projName << " found " << "(requested histogram " << calHistName << " in " << calAnaName << ")"); } else { MSG_INFO("Found calibration histogram " << sel << " " << genhist->path()); cproj.add(PercentileProjection(proj, genhist, increasing), sel); } } else if ( sel == "IMP" ) { Histo1DPtr imphist = getAnalysisObject(calAnaName, calHistName + "_IMP"); if ( !imphist || imphist->numEntries() <= 1 ) { MSG_WARNING("No impact parameter calibration histogram for " << "CentralityProjection " << projName << " found " << "(requested histogram " << calHistName << "_IMP in " << calAnaName << ")"); } else { MSG_INFO("Found calibration histogram " << sel << " " << imphist->path()); cproj.add(PercentileProjection(ImpactParameterProjection(), imphist, true), sel); } } else if ( sel == "USR" ) { #if HEPMC_VERSION_CODE >= 3000000 Histo1DPtr usrhist = getAnalysisObject(calAnaName, calHistName + "_USR"); if ( !usrhist || usrhist->numEntries() <= 1 ) { MSG_WARNING("No user-defined calibration histogram for " << "CentralityProjection " << projName << " found " << "(requested histogram " << calHistName << "_USR in " << calAnaName << ")"); continue; } else { MSG_INFO("Found calibration histogram " << sel << " " << usrhist->path()); cproj.add((UserCentEstimate(), usrhist, true), sel); } #else MSG_WARNING("UserCentEstimate is only available with HepMC3."); #endif } else if ( sel == "RAW" ) { #if HEPMC_VERSION_CODE >= 3000000 cproj.add(GeneratedCentrality(), sel); #else MSG_WARNING("GeneratedCentrality is only available with HepMC3."); #endif } else MSG_WARNING("'" << sel << "' is not a valid PercentileProjection tag."); if ( cproj.empty() ) MSG_WARNING("CentralityProjection " << projName << " did not contain any valid PercentileProjections."); return declare(cproj, projName); } } diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc --- a/src/Core/AnalysisHandler.cc +++ b/src/Core/AnalysisHandler.cc @@ -1,591 +1,582 @@ // -*- 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" namespace Rivet { AnalysisHandler::AnalysisHandler(const string& runname) : _runname(runname), _eventcounter("/_EVTCOUNT"), _xs(NAN), _xserr(NAN), _initialised(false), _ignoreBeams(false), _dumpPeriod(0), _dumping(false) { } AnalysisHandler::~AnalysisHandler() { } Log& AnalysisHandler::getLog() const { return Log::getLog("Rivet.Analysis.Handler"); } void AnalysisHandler::init(const GenEvent& ge) { if (_initialised) throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!"); setRunBeams(Rivet::beams(ge)); MSG_DEBUG("Initialising the analysis handler"); _eventcounter.reset(); // 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 (toUpper(a->status()) == "PRELIMINARY") { MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!"); } else if (toUpper(a->status()) == "OBSOLETE") { MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!"); } else if (toUpper(a->status()).find("UNVALIDATED") != string::npos) { MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!"); } } // Initialize the remaining analyses 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()); } _initialised = true; MSG_DEBUG("Analysis handler initialised"); } 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 Event event(ge); // Weights /// @todo Drop this / just report first weight when we support multiweight events _eventcounter.fill(event.weight()); MSG_DEBUG("Event #" << _eventcounter.numEntries() << " weight = " << event.weight()); // Cross-section #ifdef HEPMC_HAS_CROSS_SECTION if (ge.cross_section()) { _xs = ge.cross_section()->cross_section(); _xserr = ge.cross_section()->cross_section_error(); } #endif // 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()%_dumpPeriod == 0 ) { MSG_INFO("Dumping intermediate results to " << _dumpFile << "."); _dumping = true; finalize(); _dumping = false; writeData(_dumpFile); } } 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; // First we make copies of all analysis objects. map backupAOs; - for (auto ao : getData(false, true) ) + for (auto ao : getData(false, true, false) ) backupAOs[ao->path()] = AnalysisObjectPtr(ao->newclone()); // Now we run the (re-entrant) finalize() functions for all analyses. MSG_INFO("Finalising analyses"); for (AnaHandle a : _analyses) { a->setCrossSection(_xs); try { if ( !_dumping || a->info().reentrant() ) a->finalize(); else if ( _dumping ) MSG_INFO("Skipping periodic dump of " << a->name() << " as it is not declared reentrant."); } catch (const Error& err) { cerr << "Error in " << a->name() << "::finalize method: " << err.what() << endl; exit(1); } } // Now we copy all analysis objects to the list of finalized // ones, and restore the value to their original ones. _finalizedAOs.clear(); - for ( auto ao : getData() ) + for ( auto ao : getData(false, false, false) ) _finalizedAOs.push_back(AnalysisObjectPtr(ao->newclone())); - for ( auto ao : getData(false, true) ) { + for ( auto ao : getData(false, true, false) ) { // TODO: This should be possible to do in a nicer way, with a flag etc. if (ao->path().find("/FINAL") != std::string::npos) continue; auto aoit = backupAOs.find(ao->path()); if ( aoit == backupAOs.end() ) { AnaHandle ana = analysis(split(ao->path(), "/")[0]); if ( ana ) ana->removeAnalysisObject(ao->path()); } else copyao(aoit->second, ao); } // Print out number of events processed const int nevts = _eventcounter.numEntries(); MSG_INFO("Processed " << nevts << " event" << (nevts != 1 ? "s" : "")); // // Delete analyses // MSG_DEBUG("Deleting analyses"); // _analyses.clear(); // 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.insert(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) { std::shared_ptr toremove; for (const AnaHandle a : _analyses) { if (a->name() == analysisname) { toremove = a; break; } } if (toremove.get() != 0) { MSG_DEBUG("Removing analysis '" << analysisname << "'"); _analyses.erase(toremove); } return *this; } ///////////////////////////// void AnalysisHandler::addData(const std::vector& aos) { for (const 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); a->addAnalysisObject(ao); /// @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(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) { vector< vector > aosv; vector xsecs; vector xsecerrs; vector sows; set ananames; _eventcounter.reset(); // First scan all files and extract analysis objects and add the // corresponding anayses.. for ( auto file : aofiles ) { Scatter1DPtr xsec; CounterPtr sow; // 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 (AnalysisObject* aor : aos_raw) { AnalysisObjectPtr ao = AnalysisObjectPtr(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 += *sow; sows.push_back(sow); aosv.push_back(aos); } catch (...) { //< YODA::ReadError& throw UserError("Unexpected error in reading file: " + file); } } // 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]; } 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); } // 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; + for ( auto ao : getData(false, true, false) ) 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(); } void AnalysisHandler::readData(const string& filename) { vector aos; try { /// @todo Use new YODA SFINAE to fill the smart ptr vector directly vector aos_raw; YODA::read(filename, aos_raw); for (AnalysisObject* aor : aos_raw) aos.push_back(AnalysisObjectPtr(aor)); } catch (...) { //< YODA::ReadError& throw UserError("Unexpected error in reading file: " + filename); } if (!aos.empty()) addData(aos); } vector AnalysisHandler:: - getData(bool includeorphans, bool includetmps) const { + getData(bool includeorphans, bool includetmps, bool usefinalized) const { vector rtn; // Event counter rtn.push_back( make_shared(_eventcounter) ); // Cross-section + err as scatter YODA::Scatter1D::Points pts; pts.insert(YODA::Point1D(_xs, _xserr)); rtn.push_back( make_shared(pts, "/_XSEC") ); // Analysis histograms - for (const AnaHandle a : analyses()) { - vector aos = a->analysisObjects(); - // MSG_WARNING(a->name() << " " << aos.size()); - for (const AnalysisObjectPtr ao : aos) { - // Exclude paths from final write-out if they contain a "TMP" layer (i.e. matching "/TMP/") - /// @todo This needs to be much more nuanced for re-entrant histogramming - if ( !includetmps && ao->path().find("/TMP/" ) != string::npos) continue; - rtn.push_back(ao); + vector aos; + if (usefinalized) + aos = _finalizedAOs; + else { + for (const AnaHandle a : analyses()) { + // MSG_WARNING(a->name() << " " << aos.size()); + for (const AnalysisObjectPtr ao : a->analysisObjects()) { + aos.push_back(ao); + } } } + for (const AnalysisObjectPtr ao : aos) { + // Exclude paths from final write-out if they contain a "TMP" layer (i.e. matching "/TMP/") + /// @todo This needs to be much more nuanced for re-entrant histogramming + if ( !includetmps && ao->path().find("/TMP/" ) != string::npos) continue; + rtn.push_back(ao); + } // Sort histograms alphanumerically by path before write-out sort(rtn.begin(), rtn.end(), [](AnalysisObjectPtr a, AnalysisObjectPtr b) {return a->path() < b->path();}); if ( includeorphans ) rtn.insert(rtn.end(), _orphanedPreloads.begin(), _orphanedPreloads.end()); return rtn; } - - vector AnalysisHandler:: - getFinalizedData(bool includeorphans) const { - // Finalized analysis histograms - vector rtn = _finalizedAOs; - // Event counter - rtn.push_back( make_shared(_eventcounter) ); - // Cross-section + err as scatter - YODA::Scatter1D::Points pts; pts.insert(YODA::Point1D(_xs, _xserr)); - rtn.push_back( make_shared(pts, "/_XSEC") ); - // Sort histograms alphanumerically by path before write-out - sort(rtn.begin(), rtn.end(), [](AnalysisObjectPtr a, AnalysisObjectPtr b) {return a->path() < b->path();}); - if ( includeorphans ) - rtn.insert(rtn.end(), _orphanedPreloads.begin(), _orphanedPreloads.end()); - return rtn; - } void AnalysisHandler::writeData(const string& filename) const { vector out = _finalizedAOs; out.reserve(2*out.size()); - vector aos = getData(false, true); + vector aos = getData(false, true, false); for ( auto ao : aos ) { ao = AnalysisObjectPtr(ao->newclone()); ao->setPath("/RAW" + ao->path()); out.push_back(ao); } try { YODA::write(filename, out.begin(), out.end()); } catch (...) { //< YODA::WriteError& throw UserError("Unexpected error in writing file: " + filename); } } std::vector AnalysisHandler::analysisNames() const { std::vector rtn; for (AnaHandle a : _analyses) { rtn.push_back(a->name()); } return rtn; } const AnaHandle AnalysisHandler::analysis(const std::string& analysisname) const { for (const AnaHandle a : analyses()) if (a->name() == analysisname) return a; throw Error("No analysis named '" + analysisname + "' registered in AnalysisHandler"); } 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; } bool AnalysisHandler::needCrossSection() const { bool rtn = false; for (const AnaHandle a : _analyses) { if (!rtn) rtn = a->needsCrossSection(); if (rtn) break; } return rtn; } AnalysisHandler& AnalysisHandler::setCrossSection(double xs, double xserr) { _xs = xs; _xserr = xserr; return *this; } bool AnalysisHandler::hasCrossSection() const { return (!std::isnan(crossSection())); } AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) { analysis->_analysishandler = this; _analyses.insert(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; } }