diff --git a/include/Rivet/Analysis.hh b/include/Rivet/Analysis.hh --- a/include/Rivet/Analysis.hh +++ b/include/Rivet/Analysis.hh @@ -1,1040 +1,1055 @@ // -*- C++ -*- #ifndef RIVET_Analysis_HH #define RIVET_Analysis_HH #include "Rivet/Config/RivetCommon.hh" #include "Rivet/AnalysisInfo.hh" #include "Rivet/Event.hh" #include "Rivet/Projection.hh" #include "Rivet/ProjectionApplier.hh" #include "Rivet/ProjectionHandler.hh" #include "Rivet/AnalysisLoader.hh" #include "Rivet/Tools/RivetYODA.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Tools/ParticleUtils.hh" #include "Rivet/Tools/Cuts.hh" /// @def vetoEvent /// Preprocessor define for vetoing events, including the log message and return. #define vetoEvent \ do { MSG_DEBUG("Vetoing event on line " << __LINE__ << " of " << __FILE__); return; } while(0) namespace Rivet { // Forward declaration class AnalysisHandler; /// @brief This is the base class of all analysis classes in Rivet. /// /// There are /// three virtual functions which should be implemented in base classes: /// /// void init() is called by Rivet before a run is started. Here the /// analysis class should book necessary histograms. The needed /// projections should probably rather be constructed in the /// constructor. /// /// void analyze(const Event&) is called once for each event. Here the /// analysis class should apply the necessary Projections and fill the /// histograms. /// /// void finalize() is called after a run is finished. Here the analysis /// class should do whatever manipulations are necessary on the /// histograms. Writing the histograms to a file is, however, done by /// the Rivet class. class Analysis : public ProjectionApplier { /// The AnalysisHandler is a friend. friend class AnalysisHandler; public: /// @name Standard constructors and destructors. //@{ // /// The default constructor. // Analysis(); /// Constructor Analysis(const std::string& name); /// The destructor. virtual ~Analysis() {} //@} public: /// @name Main analysis methods //@{ /// Initialize this analysis object. A concrete class should here /// book all necessary histograms. An overridden function must make /// sure it first calls the base class function. virtual void init() { } /// Analyze one event. A concrete class should here apply the /// necessary projections on the \a event and fill the relevant /// histograms. An overridden function must make sure it first calls /// the base class function. virtual void analyze(const Event& event) = 0; /// Finalize this analysis object. A concrete class should here make /// all necessary operations on the histograms. Writing the /// histograms to a file is, however, done by the Rivet class. An /// overridden function must make sure it first calls the base class /// function. virtual void finalize() { } - + + /// A new method for postprocessing analyses is introduced. + /// This function can be called to handle e. g. merging analyses from + /// different runs. Besides finalize() this method can be called as + /// often as wanted because it only depends on finished processed + /// histograms which stay itself untouched but are used to produce + /// further output. This method is virtual because it is considered + /// to be overwritten in inherited classes. + virtual void post() { } + + /// Replaces the content of booked histograms if there were those + /// with the same path in _analysisObjectsRead. So far, the method is + /// called in AnalysisHandler::init method such that the replacement + /// is done automatically. This can of course be changed. + void replaceByData( std::map< std::string, AnalysisObjectPtr > readObjects); + //@} public: /// @name Metadata /// Metadata is used for querying from the command line and also for /// building web pages and the analysis pages in the Rivet manual. //@{ /// Get the actual AnalysisInfo object in which all this metadata is stored. const AnalysisInfo& info() const { assert(_info && "No AnalysisInfo object :O"); return *_info; } /// @brief Get the name of the analysis. /// /// By default this is computed by combining the results of the experiment, /// year and Spires ID metadata methods and you should only override it if /// there's a good reason why those won't work. virtual std::string name() const { return (info().name().empty()) ? _defaultname : info().name(); } /// Get the Inspire ID code for this analysis. virtual std::string inspireId() const { return info().inspireId(); } /// Get the SPIRES ID code for this analysis (~deprecated). virtual std::string spiresId() const { return info().spiresId(); } /// @brief Names & emails of paper/analysis authors. /// /// Names and email of authors in 'NAME \' format. The first /// name in the list should be the primary contact person. virtual std::vector authors() const { return info().authors(); } /// @brief Get a short description of the analysis. /// /// Short (one sentence) description used as an index entry. /// Use @a description() to provide full descriptive paragraphs /// of analysis details. virtual std::string summary() const { return info().summary(); } /// @brief Get a full description of the analysis. /// /// Full textual description of this analysis, what it is useful for, /// what experimental techniques are applied, etc. Should be treated /// as a chunk of restructuredText (http://docutils.sourceforge.net/rst.html), /// with equations to be rendered as LaTeX with amsmath operators. virtual std::string description() const { return info().description(); } /// @brief Information about the events needed as input for this analysis. /// /// Event types, energies, kinematic cuts, particles to be considered /// stable, etc. etc. Should be treated as a restructuredText bullet list /// (http://docutils.sourceforge.net/rst.html) virtual std::string runInfo() const { return info().runInfo(); } /// Experiment which performed and published this analysis. virtual std::string experiment() const { return info().experiment(); } /// Collider on which the experiment ran. virtual std::string collider() const { return info().collider(); } /// When the original experimental analysis was published. virtual std::string year() const { return info().year(); } /// The luminosity in inverse femtobarn virtual std::string luminosityfb() const { return info().luminosityfb(); } /// Journal, and preprint references. virtual std::vector references() const { return info().references(); } /// BibTeX citation key for this article. virtual std::string bibKey() const { return info().bibKey(); } /// BibTeX citation entry for this article. virtual std::string bibTeX() const { return info().bibTeX(); } /// Whether this analysis is trusted (in any way!) virtual std::string status() const { return (info().status().empty()) ? "UNVALIDATED" : info().status(); } /// Any work to be done on this analysis. virtual std::vector todos() const { return info().todos(); } /// Return the allowed pairs of incoming beams required by this analysis. virtual const std::vector& requiredBeams() const { return info().beams(); } /// Declare the allowed pairs of incoming beams required by this analysis. virtual Analysis& setRequiredBeams(const std::vector& requiredBeams) { info().setBeams(requiredBeams); return *this; } /// Sets of valid beam energy pairs, in GeV virtual const std::vector >& requiredEnergies() const { return info().energies(); } /// Get vector of analysis keywords virtual const std::vector & keywords() const { return info().keywords(); } /// Declare the list of valid beam energy pairs, in GeV virtual Analysis& setRequiredEnergies(const std::vector >& requiredEnergies) { info().setEnergies(requiredEnergies); return *this; } /// Return true if this analysis needs to know the process cross-section. /// @todo Remove this and require HepMC >= 2.06 bool needsCrossSection() const { return info().needsCrossSection(); } /// Declare whether this analysis needs to know the process cross-section from the generator. /// @todo Remove this and require HepMC >= 2.06 Analysis& setNeedsCrossSection(bool needed=true) { info().setNeedsCrossSection(needed); return *this; } //@} /// @name Internal metadata modifying methods //@{ /// Get the actual AnalysisInfo object in which all this metadata is stored (non-const). AnalysisInfo& info() { assert(_info && "No AnalysisInfo object :O"); return *_info; } //@} /// @name Run conditions //@{ /// Incoming beams for this run const ParticlePair& beams() const; /// Incoming beam IDs for this run const PdgIdPair beamIds() const; /// Centre of mass energy for this run double sqrtS() const; //@} /// @name Analysis / beam compatibility testing //@{ /// Check if analysis is compatible with the provided beam particle IDs and energies bool isCompatible(const ParticlePair& beams) const; /// Check if analysis is compatible with the provided beam particle IDs and energies bool isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const; /// Check if analysis is compatible with the provided beam particle IDs and energies bool isCompatible(const PdgIdPair& beams, const std::pair& energies) const; //@} /// Set the cross section from the generator Analysis& setCrossSection(double xs); /// Access the controlling AnalysisHandler object. AnalysisHandler& handler() const { return *_analysishandler; } protected: /// Get a Log object based on the name() property of the calling analysis object. Log& getLog() const; /// Get the process cross-section in pb. Throws if this hasn't been set. double crossSection() const; /// Get the process cross-section per generated event in pb. Throws if this /// hasn't been set. double crossSectionPerEvent() const; /// Get the number of events seen (via the analysis handler). Use in the /// finalize phase only. size_t numEvents() const; /// Get the sum of event weights seen (via the analysis handler). Use in the /// finalize phase only. double sumOfWeights() const; protected: /// @name Histogram paths //@{ /// Get the canonical histogram "directory" path for this analysis. const std::string histoDir() const; /// Get the canonical histogram path for the named histogram in this analysis. const std::string histoPath(const std::string& hname) const; /// Get the canonical histogram path for the numbered histogram in this analysis. const std::string histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const; /// Get the internal histogram name for given d, x and y (cf. HepData) const std::string makeAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const; //@} /// @name Histogram reference data //@{ /// Get reference data for a named histo /// @todo Move to the templated version when we have C++11 and can have a default fn template type const YODA::Scatter2D& refData(const string& hname) const; /// Get reference data for a numbered histo /// @todo Move to the templated version when we have C++11 and can have a default fn template type const YODA::Scatter2D& refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const; /// Get reference data for a named histo /// @todo Would be nice to just use these and ditch the S2D no-template version, /// but we need C++11 for default args in function templates // template /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject? template const T& refData(const string& hname) const { _cacheRefData(); MSG_TRACE("Using histo bin edges for " << name() << ":" << hname); if (!_refdata[hname]) { MSG_ERROR("Can't find reference histogram " << hname); throw Exception("Reference data " + hname + " not found."); } return dynamic_cast(*_refdata[hname]); } /// Get reference data for a numbered histo /// @todo Would be nice to just use these and ditch the S2D no-template version, /// but we need C++11 for default args in function templates // template /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject? template const T& refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { const string hname = makeAxisCode(datasetId, xAxisId, yAxisId); return refData(hname); } //@} /// @name Counter booking //@{ /// Book a counter. CounterPtr bookCounter(const std::string& name, const std::string& title=""); // const std::string& valtitle="" /// Book a counter, using a path generated from the dataset and axis ID codes /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. CounterPtr bookCounter(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title=""); // const std::string& valtitle="" //@} /// @name 1D histogram booking //@{ /// Book a 1D histogram with @a nbins uniformly distributed across the range @a lower - @a upper . Histo1DPtr bookHisto1D(const std::string& name, size_t nbins, double lower, double upper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram with non-uniform bins defined by the vector of bin edges @a binedges . Histo1DPtr bookHisto1D(const std::string& name, const std::vector& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram with binning from a reference scatter. Histo1DPtr bookHisto1D(const std::string& name, const Scatter2D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram, using the binnings in the reference data histogram. Histo1DPtr bookHisto1D(const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Histo1DPtr bookHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); //@} /// @name 2D histogram booking //@{ /// Book a 2D histogram with @a nxbins and @a nybins uniformly /// distributed across the ranges @a xlower - @a xupper and @a /// ylower - @a yupper respectively along the x- and y-axis. Histo2DPtr bookHisto2D(const std::string& name, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram with non-uniform bins defined by the /// vectorx of bin edges @a xbinedges and @a ybinedges. Histo2DPtr bookHisto2D(const std::string& name, const std::vector& xbinedges, const std::vector& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); // /// Book a 2D histogram with binning from a reference scatter. // Histo2DPtr bookHisto2D(const std::string& name, // const Scatter3D& refscatter, // const std::string& title="", // const std::string& xtitle="", // const std::string& ytitle="", // const std::string& ztitle=""); // /// Book a 2D histogram, using the binnings in the reference data histogram. // Histo2DPtr bookHisto2D(const std::string& name, // const std::string& title="", // const std::string& xtitle="", // const std::string& ytitle="", // const std::string& ztitle=""); // /// Book a 2D histogram, using the binnings in the reference data histogram. // /// // /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. // Histo2DPtr bookHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, // const std::string& title="", // const std::string& xtitle="", // const std::string& ytitle="", // const std::string& ztitle=""); //@} /// @name 1D profile histogram booking //@{ /// Book a 1D profile histogram with @a nbins uniformly distributed across the range @a lower - @a upper . Profile1DPtr bookProfile1D(const std::string& name, size_t nbins, double lower, double upper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges . Profile1DPtr bookProfile1D(const std::string& name, const std::vector& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram with binning from a reference scatter. Profile1DPtr bookProfile1D(const std::string& name, const Scatter2D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram, using the binnings in the reference data histogram. Profile1DPtr bookProfile1D(const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Profile1DPtr bookProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); //@} /// @name 2D profile histogram booking //@{ /// Book a 2D profile histogram with @a nxbins and @a nybins uniformly /// distributed across the ranges @a xlower - @a xupper and @a ylower - @a /// yupper respectively along the x- and y-axis. Profile2DPtr bookProfile2D(const std::string& name, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram with non-uniform bins defined by the vectorx /// of bin edges @a xbinedges and @a ybinedges. Profile2DPtr bookProfile2D(const std::string& name, const std::vector& xbinedges, const std::vector& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram with binning from a reference scatter. // Profile2DPtr bookProfile2D(const std::string& name, // const Scatter3D& refscatter, // const std::string& title="", // const std::string& xtitle="", // const std::string& ytitle="", // const std::string& ztitle=""); // /// Book a 2D profile histogram, using the binnings in the reference data histogram. // Profile2DPtr bookProfile2D(const std::string& name, // const std::string& title="", // const std::string& xtitle="", // const std::string& ytitle="", // const std::string& ztitle=""); // /// Book a 2D profile histogram, using the binnings in the reference data histogram. // /// // /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. // Profile2DPtr bookProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, // const std::string& title="", // const std::string& xtitle="", // const std::string& ytitle="", // const std::string& ztitle=""); //@} /// @name 2D scatter booking //@{ /// @brief Book a 2-dimensional data point set with the given name. /// /// @note Unlike histogram booking, scatter booking by default makes no /// attempt to use reference data to pre-fill the data object. If you want /// this, which is sometimes useful e.g. when the x-position is not really /// meaningful and can't be extracted from the data, then set the @a /// copy_pts parameter to true. This creates points to match the reference /// data's x values and errors, but with the y values and errors zeroed... /// assuming that there is a reference histo with the same name: if there /// isn't, an exception will be thrown. Scatter2DPtr bookScatter2D(const std::string& name, bool copy_pts=false, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// @brief Book a 2-dimensional data point set, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. /// /// @note Unlike histogram booking, scatter booking by default makes no /// attempt to use reference data to pre-fill the data object. If you want /// this, which is sometimes useful e.g. when the x-position is not really /// meaningful and can't be extracted from the data, then set the @a /// copy_pts parameter to true. This creates points to match the reference /// data's x values and errors, but with the y values and errors zeroed. Scatter2DPtr bookScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, bool copy_pts=false, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// @brief Book a 2-dimensional data point set with equally spaced x-points in a range. /// /// The y values and errors will be set to 0. Scatter2DPtr bookScatter2D(const std::string& name, size_t npts, double lower, double upper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// @brief Book a 2-dimensional data point set based on provided contiguous "bin edges". /// /// The y values and errors will be set to 0. Scatter2DPtr bookScatter2D(const std::string& hname, const std::vector& binedges, const std::string& title, const std::string& xtitle, const std::string& ytitle); //@} public: /// @name Analysis object manipulation /// @todo Should really be protected: only public to keep BinnedHistogram happy for now... //@{ /// Multiplicatively scale the given counter, @a cnt, by factor @s factor. void scale(CounterPtr cnt, double factor); /// Multiplicatively scale the given counters, @a cnts, by factor @s factor. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of CounterPtrs void scale(const std::vector& cnts, double factor) { for (auto& c : cnts) scale(c, factor); } /// @todo YUCK! template void scale(const CounterPtr (&cnts)[array_size], double factor) { // for (size_t i = 0; i < std::extent::value; ++i) scale(cnts[i], factor); for (auto& c : cnts) scale(c, factor); } /// Normalize the given histogram, @a histo, to area = @a norm. void normalize(Histo1DPtr histo, double norm=1.0, bool includeoverflows=true); /// Normalize the given histograms, @a histos, to area = @a norm. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo1DPtrs void normalize(const std::vector& histos, double norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// @todo YUCK! template void normalize(const Histo1DPtr (&histos)[array_size], double norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// Multiplicatively scale the given histogram, @a histo, by factor @s factor. void scale(Histo1DPtr histo, double factor); /// Multiplicatively scale the given histograms, @a histos, by factor @s factor. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo1DPtrs void scale(const std::vector& histos, double factor) { for (auto& h : histos) scale(h, factor); } /// @todo YUCK! template void scale(const Histo1DPtr (&histos)[array_size], double factor) { for (auto& h : histos) scale(h, factor); } /// Normalize the given histogram, @a histo, to area = @a norm. void normalize(Histo2DPtr histo, double norm=1.0, bool includeoverflows=true); /// Normalize the given histograms, @a histos, to area = @a norm. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo2DPtrs void normalize(const std::vector& histos, double norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// @todo YUCK! template void normalize(const Histo2DPtr (&histos)[array_size], double norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// Multiplicatively scale the given histogram, @a histo, by factor @s factor. void scale(Histo2DPtr histo, double factor); /// Multiplicatively scale the given histograms, @a histos, by factor @s factor. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo2DPtrs void scale(const std::vector& histos, double factor) { for (auto& h : histos) scale(h, factor); } /// @todo YUCK! template void scale(const Histo2DPtr (&histos)[array_size], double factor) { for (auto& h : histos) scale(h, factor); } /// Helper for counter division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const; /// Helper for histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Counter& c1, const YODA::Counter& c2, Scatter1DPtr s) const; /// Helper for histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const; /// Helper for histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const; /// Helper for profile histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const; /// Helper for profile histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Profile1D& p1, const YODA::Profile1D& p2, Scatter2DPtr s) const; /// Helper for 2D histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const; /// Helper for 2D histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Histo2D& h1, const YODA::Histo2D& h2, Scatter3DPtr s) const; /// Helper for 2D profile histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const; /// Helper for 2D profile histogram division with raw YODA objects /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Profile2D& p1, const YODA::Profile2D& p2, Scatter3DPtr s) const; /// Helper for histogram efficiency calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const; /// Helper for histogram efficiency calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void efficiency(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const; /// Helper for histogram asymmetry calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const; /// Helper for histogram asymmetry calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void asymm(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const; /// Helper for converting a differential histo to an integral one. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void integrate(Histo1DPtr h, Scatter2DPtr s) const; /// Helper for converting a differential histo to an integral one. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void integrate(const Histo1D& h, Scatter2DPtr s) const; //@} public: /// List of registered analysis data objects const vector& analysisObjects() const { return _analysisobjects; } protected: /// @name Data object registration, retrieval, and removal //@{ /// Register a data object in the histogram system void addAnalysisObject(AnalysisObjectPtr ao); /// Get a data object from the histogram system /// @todo Use this default function template arg in C++11 // template template const std::shared_ptr getAnalysisObject(const std::string& name) const { foreach (const AnalysisObjectPtr& ao, analysisObjects()) { if (ao->path() == histoPath(name)) return dynamic_pointer_cast(ao); } throw Exception("Data object " + histoPath(name) + " not found"); } /// Get a data object from the histogram system (non-const) /// @todo Use this default function template arg in C++11 // template template std::shared_ptr getAnalysisObject(const std::string& name) { foreach (const AnalysisObjectPtr& ao, analysisObjects()) { if (ao->path() == histoPath(name)) return dynamic_pointer_cast(ao); } throw Exception("Data object " + histoPath(name) + " not found"); } /// Unregister a data object from the histogram system (by name) void removeAnalysisObject(const std::string& path); /// Unregister a data object from the histogram system (by pointer) void removeAnalysisObject(AnalysisObjectPtr ao); /// Get a named Histo1D object from the histogram system const Histo1DPtr getHisto1D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Histo1D object from the histogram system (non-const) Histo1DPtr getHisto1D(const std::string& name) { return getAnalysisObject(name); } /// Get a Histo1D object from the histogram system by axis ID codes (non-const) const Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Histo1D object from the histogram system by axis ID codes (non-const) Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } // /// Get a named Histo2D object from the histogram system // const Histo2DPtr getHisto2D(const std::string& name) const { // return getAnalysisObject(name); // } // /// Get a named Histo2D object from the histogram system (non-const) // Histo2DPtr getHisto2D(const std::string& name) { // return getAnalysisObject(name); // } // /// Get a Histo2D object from the histogram system by axis ID codes (non-const) // const Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } // /// Get a Histo2D object from the histogram system by axis ID codes (non-const) // Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } /// Get a named Profile1D object from the histogram system const Profile1DPtr getProfile1D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Profile1D object from the histogram system (non-const) Profile1DPtr getProfile1D(const std::string& name) { return getAnalysisObject(name); } /// Get a Profile1D object from the histogram system by axis ID codes (non-const) const Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Profile1D object from the histogram system by axis ID codes (non-const) Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } // /// Get a named Profile2D object from the histogram system // const Profile2DPtr getProfile2D(const std::string& name) const { // return getAnalysisObject(name); // } // /// Get a named Profile2D object from the histogram system (non-const) // Profile2DPtr getProfile2D(const std::string& name) { // return getAnalysisObject(name); // } // /// Get a Profile2D object from the histogram system by axis ID codes (non-const) // const Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } // /// Get a Profile2D object from the histogram system by axis ID codes (non-const) // Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } /// Get a named Scatter2D object from the histogram system const Scatter2DPtr getScatter2D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Scatter2D object from the histogram system (non-const) Scatter2DPtr getScatter2D(const std::string& name) { return getAnalysisObject(name); } /// Get a Scatter2D object from the histogram system by axis ID codes (non-const) const Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Scatter2D object from the histogram system by axis ID codes (non-const) Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } //@} private: /// Name passed to constructor (used to find .info analysis data file, and as a fallback) string _defaultname; /// Pointer to analysis metadata object unique_ptr _info; /// Storage of all plot objects /// @todo Make this a map for fast lookup by path? vector _analysisobjects; - + bool _haveReadData; /// @name Cross-section variables //@{ double _crossSection; bool _gotCrossSection; //@} /// The controlling AnalysisHandler object. AnalysisHandler* _analysishandler; /// Collection of cached refdata to speed up many autobookings: the /// reference data file should only be read once. mutable std::map _refdata; private: /// @name Utility functions //@{ /// Get the reference data for this paper and cache it. void _cacheRefData() const; //@} /// The assignment operator is private and must never be called. /// In fact, it should not even be implemented. Analysis& operator=(const Analysis&); }; } // Include definition of analysis plugin system so that analyses automatically see it when including Analysis.hh #include "Rivet/AnalysisBuilder.hh" /// @def DECLARE_RIVET_PLUGIN /// Preprocessor define to prettify the global-object plugin hook mechanism. #define DECLARE_RIVET_PLUGIN(clsname) Rivet::AnalysisBuilder plugin_ ## clsname /// @def DECLARE_ALIASED_RIVET_PLUGIN /// Preprocessor define to prettify the global-object plugin hook mechanism, with an extra alias name for this analysis. // #define DECLARE_ALIASED_RIVET_PLUGIN(clsname, alias) Rivet::AnalysisBuilder plugin_ ## clsname ## ( ## #alias ## ) #define DECLARE_ALIASED_RIVET_PLUGIN(clsname, alias) DECLARE_RIVET_PLUGIN(clsname)( #alias ) /// @def DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR /// Preprocessor define to prettify the manky constructor with name string argument #define DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR(clsname) clsname() : Analysis(# clsname) {} /// @def DEFAULT_RIVET_ANALYSIS_CTOR /// Slight abbreviation for DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR #define DEFAULT_RIVET_ANALYSIS_CTOR(clsname) DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR(clsname) #endif diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh --- a/include/Rivet/AnalysisHandler.hh +++ b/include/Rivet/AnalysisHandler.hh @@ -1,250 +1,255 @@ // -*- 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) { 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; /// 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; /// Get the sum of the event weights seen - the weighted equivalent of the /// number of events. Should only really be used by external steering code /// or analyses in the finalize phase. double sumOfWeights() const; /// Set 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) void setSumOfWeights(const double& sum); /// Is cross-section information required by at least one child analysis? bool needCrossSection() const; /// Set the cross-section for the process being generated. AnalysisHandler& setCrossSection(double xs); /// Get the cross-section known to the handler. double crossSection() const { return _xs; } /// Whether the handler knows about a cross-section. bool hasCrossSection() const; /// 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 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(); + /// Function for postprocessing is introduced. This allows for further processing + /// e. g. additional analyses which were done in the past but should be included + /// in a certain analysis in one way or the other + void post(); + /// Method reads YODA files and adds the contained AnalysisObjects /// to the member map "_readObjects". void readData(const std::string& filename); //@} /// @name Histogram / data object access //@{ /// Get all analyses' plots as a vector of analysis objects. std::vector getData() const; /// Write all analyses' plots to the named file. void writeData(const std::string& filename) const; //@} private: - + /// Map containing read YODA objects. std::map< std::string, AnalysisObjectPtr > _readObjects; bool _haveReadData; /// The collection of Analysis objects to be used. set _analyses; /// @name Run properties //@{ /// Run name std::string _runname; /// Number of events seen. /// @todo Replace by a counter unsigned int _numEvents; /// Sum of event weights seen. /// @todo Replace by a counter double _sumOfWeights, _sumOfWeightsSq; /// 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; //@} 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,818 +1,885 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Analysis.hh" #include "Rivet/AnalysisHandler.hh" #include "Rivet/AnalysisInfo.hh" #include "Rivet/Tools/BeamConstraint.hh" namespace Rivet { Analysis::Analysis(const string& name) - : _crossSection(-1.0), + : _haveReadData(false), + _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); } + + + /// xyz This method replaces the content of all those pointers in "_analysisobjects" which have + /// a partner in "objectsRead". The AnalysisObjectPtr are matched by their path. + void Analysis::replaceByData( std::map< std::string, AnalysisObjectPtr > readObjects) { + if( _haveReadData ) { + MSG_WARNING("Analysis " << name() << " already got data. This is the second time and the older objects will be overwritten!"); + } + + for( AnalysisObjectPtr ao : _analysisobjects ) { + + if( !readObjects[ao->path()] ) { + continue; + } + + /// xyz type cast from AnalysisObjectPtr to respective pointer in order to check for number of + // entries and points, respectively. Further, the operator "=" for YODA::AnalysisObject does + // only copy the path/names and not the entries. However, this is done by the operator respective + // "=" of YODA::Histo1D etc. + // @note may not be the most efficient solution + Histo1DPtr tmpHisto1D, tmpHisto1DRead; + Histo2DPtr tmpHisto2D, tmpHisto2DRead; + Profile1DPtr tmpProfile1D, tmpProfile1DRead; + Profile2DPtr tmpProfile2D, tmpProfile2DRead; + Scatter1DPtr tmpScatter1D, tmpScatter1DRead; + Scatter2DPtr tmpScatter2D, tmpScatter2DRead; + Scatter3DPtr tmpScatter3D, tmpScatter3DRead; + + // assume that type casting in the following is successful. Will be set to false if none of the + // type castings was successful + bool replaced = true; + + try { + + if( ( tmpHisto1D = dynamic_pointer_cast< YODA::Histo1D >( ao ) ) && ( tmpHisto1DRead = dynamic_pointer_cast< YODA::Histo1D >( readObjects[ao->path()] ) ) ) { + *tmpHisto1D = *tmpHisto1DRead; + } else if ( ( tmpHisto2D = dynamic_pointer_cast< YODA::Histo2D >( ao ) ) && ( tmpHisto2DRead = dynamic_pointer_cast< YODA::Histo2D >( readObjects[ao->path()] ) ) ) { + *tmpHisto2D = *tmpHisto2DRead; + } else if ( ( tmpProfile1D = dynamic_pointer_cast< YODA::Profile1D >( ao ) ) && ( tmpProfile1DRead = dynamic_pointer_cast< YODA::Profile1D >( readObjects[ao->path()] ) ) ) { + *tmpProfile1D = *tmpProfile1DRead; + } else if ( ( tmpProfile2D = dynamic_pointer_cast< YODA::Profile2D >( ao ) ) && ( tmpProfile2DRead = dynamic_pointer_cast< YODA::Profile2D >( readObjects[ao->path()] ) ) ) { + *tmpProfile2D = *tmpProfile2DRead; + } else if ( ( tmpScatter1D = dynamic_pointer_cast< YODA::Scatter1D >( ao ) ) && ( tmpScatter1DRead = dynamic_pointer_cast< YODA::Scatter1D >( readObjects[ao->path()] ) ) ) { + *tmpScatter1D = *tmpScatter1DRead; + } else if ( ( tmpScatter2D = dynamic_pointer_cast< YODA::Scatter2D >( ao ) ) && ( tmpScatter2DRead = dynamic_pointer_cast< YODA::Scatter2D >( readObjects[ao->path()] ) ) ) { + *tmpScatter2D = *tmpScatter2DRead; + } else if ( ( tmpScatter3D = dynamic_pointer_cast< YODA::Scatter3D >( ao ) ) && ( tmpScatter3DRead = dynamic_pointer_cast< YODA::Scatter3D >( readObjects[ao->path()] ) ) ) { + *tmpScatter3D = *tmpScatter3DRead; + } else { + MSG_ERROR( "Unknown analysis object type " << ao->type() << " with path " << ao->path() ); + replaced = false; + } + + } catch(...) { + throw UserError("Replacing histograms by given data failed in analysis"); + } + + if( replaced ) { + MSG_INFO("Histogram with path " << ao->path() << " replaced"); + _haveReadData = true; // set to true because at least one object has been read + } + + } + + } + 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() + "/" + makeAxisCode(datasetId, xAxisId, yAxisId); } const string Analysis::makeAxisCode(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::sumOfWeights() const { return handler().sumOfWeights(); } /////////////////////////////////////////// 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::isCompatible(const PdgIdPair& beams, const pair& energies) const { // First check the beam IDs bool beamIdsOk = false; foreach (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; foreach (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; /// @todo Need to also check internal consistency of the analysis' /// beam requirements with those of the projections it uses. } /////////////////////////////////////////// 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(name()); } } const Scatter2D& Analysis::refData(const string& hname) const { _cacheRefData(); MSG_TRACE("Using histo bin edges for " << name() << ":" << hname); if (!_refdata[hname]) { MSG_ERROR("Can't find reference histogram " << hname); throw Exception("Reference data " + hname + " not found."); } return dynamic_cast(*_refdata[hname]); } const Scatter2D& Analysis::refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { const string hname = makeAxisCode(datasetId, xAxisId, yAxisId); return refData(hname); } CounterPtr Analysis::bookCounter(const string& cname, const string& title) { // const string& xtitle, // const string& ytitle) { const string path = histoPath(cname); CounterPtr ctr = make_shared(path, title); addAnalysisObject(ctr); MSG_TRACE("Made counter " << cname << " for " << name()); // hist->setAnnotation("XLabel", xtitle); // hist->setAnnotation("YLabel", ytitle); return ctr; } CounterPtr Analysis::bookCounter(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title) { // const string& xtitle, // const string& ytitle) { const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId); return bookCounter(axisCode, title); } Histo1DPtr Analysis::bookHisto1D(const string& hname, size_t nbins, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Histo1DPtr hist = make_shared(nbins, lower, upper, path, title); addAnalysisObject(hist); MSG_TRACE("Made histogram " << hname << " for " << name()); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); return hist; } Histo1DPtr Analysis::bookHisto1D(const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Histo1DPtr hist = make_shared(binedges, path, title); addAnalysisObject(hist); MSG_TRACE("Made histogram " << hname << " for " << name()); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); return hist; } Histo1DPtr Analysis::bookHisto1D(const string& hname, const Scatter2D& refscatter, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Histo1DPtr hist = make_shared(refscatter, path); addAnalysisObject(hist); MSG_TRACE("Made histogram " << hname << " for " << name()); hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); return 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 = makeAxisCode(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); addAnalysisObject(hist); MSG_TRACE("Made 2D histogram " << hname << " for " << name()); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); hist->setAnnotation("ZLabel", ztitle); return 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); addAnalysisObject(hist); MSG_TRACE("Made 2D histogram " << hname << " for " << name()); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); hist->setAnnotation("ZLabel", ztitle); return hist; } // 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) ); // addAnalysisObject(hist); // MSG_TRACE("Made 2D histogram " << hname << " for " << name()); // hist->setTitle(title); // hist->setAnnotation("XLabel", xtitle); // hist->setAnnotation("YLabel", ytitle); // hist->setAnnotation("ZLabel", ztitle); // return 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 = makeAxisCode(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); addAnalysisObject(prof); MSG_TRACE("Made profile histogram " << hname << " for " << name()); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); return 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); addAnalysisObject(prof); MSG_TRACE("Made profile histogram " << hname << " for " << name()); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); return prof; } 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); addAnalysisObject(prof); MSG_TRACE("Made profile histogram " << hname << " for " << name()); prof->setTitle(title); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); return 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 = makeAxisCode(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); addAnalysisObject(prof); MSG_TRACE("Made 2D profile histogram " << hname << " for " << name()); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); prof->setAnnotation("ZLabel", ztitle); return 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); addAnalysisObject(prof); MSG_TRACE("Made 2D profile histogram " << hname << " for " << name()); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); prof->setAnnotation("ZLabel", ztitle); return prof; } // 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) ); // addAnalysisObject(prof); // MSG_TRACE("Made 2D profile histogram " << hname << " for " << name()); // prof->setTitle(title); // prof->setAnnotation("XLabel", xtitle); // prof->setAnnotation("YLabel", ytitle); // prof->setAnnotation("ZLabel", ztitle); // return 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 = makeAxisCode(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 = makeAxisCode(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); foreach (Point2D& p, s->points()) p.setY(0, 0); } else { s = make_shared(path); } addAnalysisObject(s); MSG_TRACE("Made scatter " << hname << " for " << name()); s->setTitle(title); s->setAnnotation("XLabel", xtitle); s->setAnnotation("YLabel", ytitle); return 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); } addAnalysisObject(s); MSG_TRACE("Made scatter " << hname << " for " << name()); s->setTitle(title); s->setAnnotation("XLabel", xtitle); s->setAnnotation("YLabel", ytitle); return 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); } addAnalysisObject(s); MSG_TRACE("Made scatter " << hname << " for " << name()); s->setTitle(title); s->setAnnotation("XLabel", xtitle); s->setAnnotation("YLabel", ytitle); return 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 { 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 { 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; } } } } diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc --- a/src/Core/AnalysisHandler.cc +++ b/src/Core/AnalysisHandler.cc @@ -1,393 +1,431 @@ // -*- 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/WriterYODA.h" #include "YODA/ReaderYODA.h" namespace Rivet { AnalysisHandler::AnalysisHandler(const string& runname) // xyz Initialize introduced member "_haveReadData" : _haveReadData(false), _runname(runname), _numEvents(0), _sumOfWeights(0.0), _xs(NAN), _initialised(false), _ignoreBeams(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"); _numEvents = 0; _sumOfWeights = 0.0; _sumOfWeightsSq = 0.0; // 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 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 _numEvents += 1; _sumOfWeights += event.weight(); _sumOfWeightsSq += sqr(event.weight()); MSG_DEBUG("Event #" << _numEvents << " 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()); } } void AnalysisHandler::analyze(const GenEvent* ge) { if (ge == NULL) { 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"); for (AnaHandle a : _analyses) { a->setCrossSection(_xs); try { a->finalize(); } catch (const Error& err) { cerr << "Error in " << a->name() << "::finalize method: " << err.what() << endl; exit(1); } } // Print out number of events processed MSG_INFO("Processed " << _numEvents << " event" << (_numEvents == 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; } + /// This method does three things: (1) In case, the hanlder is not initialized, + /// Analysis::init() is called for each analysis, (2) if YODA objects have been + /// read before, those are passed to the analysis which replaces the content of + /// the booked objects and after all (3) Analysis::post() is called. + void AnalysisHandler::post() { + MSG_INFO("Post-processing analyses"); + if( _haveReadData ) { + MSG_INFO("Replace by data in case paths are matching"); + } + for (AnaHandle a : _analyses) { + if( !_initialised ) { + MSG_INFO( "No MC running: The post-processing is based on YODA files only." ); + try { + a->_allowProjReg = true; + a->init(); + } catch( const Error& err ) { + MSG_ERROR( "Unexpected error in " << a->name() << "::init() during the initialisation for the post-processing: " << err.what() ); + exit(1); + } + } + if( _haveReadData ) { + try { + MSG_INFO("Replacing data ..."); + a->replaceByData( _readObjects ); + } catch (const UserError& err) { + MSG_ERROR("Error in analysis " << a->name() << ": " << err.what() ) ; + } + } + try { + a->post(); + } catch (const Error& err) { + MSG_ERROR("Unexpected error in post() in analysis " << a->name() ); + } + } + MSG_INFO("Post-processing finished"); + } + 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. for (const AnaHandle& a : _analyses) { if (a->name() == analysisname) { MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate"); return *this; } } AnaHandle analysis( AnalysisLoader::getAnalysis(analysisname) ); if (analysis.get() != 0) { // < Check for null analysis. MSG_DEBUG("Adding analysis '" << analysisname << "'"); 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; } vector AnalysisHandler::getData() const { vector rtn; rtn.push_back( make_shared(YODA::Dbn0D(_numEvents, _sumOfWeights, _sumOfWeightsSq), "/_EVTCOUNT") ); YODA::Scatter1D::Points pts; pts.insert(YODA::Point1D(_xs, _xserr)); rtn.push_back( make_shared(pts, "/_XSEC") ); for (const AnaHandle a : analyses()) { vector aos = a->analysisObjects(); // MSG_WARNING(a->name() << " " << aos.size()); for (const AnalysisObjectPtr ao : aos) { // Exclude paths starting with /TMP/ from final write-out /// @todo This needs to be much more nuanced for re-entrant histogramming if (ao->path().find("/TMP/") != string::npos) continue; rtn.push_back(ao); } } sort(rtn.begin(), rtn.end(), [](AnalysisObjectPtr a, AnalysisObjectPtr b) { return a->path() < b->path(); } ); return rtn; } void AnalysisHandler::writeData(const string& filename) const { const vector aos = getData(); try { YODA::WriterYODA::write(filename, aos.begin(), aos.end()); } catch (...) { /// @todo Move to specific YODA::WriteError type when YODA >= 1.5.0 is well-established throw UserError("Unexpected error in writing file to: " + filename); } } /// xyz Reads Objects from a given YODA file and stores those in a map. No duplicates possible /// since they would be overwritten. @note Method has to be called before AnalysisHandler::init /// since there potential read objects are given to the added analyses. void AnalysisHandler::readData( const std::string& filename ) { YODA::Reader &reader = YODA::mkReader( filename ); std::vector< YODA::AnalysisObject* > aos; try { reader.read( filename, aos ); } catch(const YODA::ReadError& err) { MSG_ERROR( "Error in AnalysisHandler::readData: " << err.what() ); exit(1); } MSG_INFO( "Read data from " << filename ); for( const auto ao : aos ) { //MSG_INFO( "type of ao is " << typeid( ao ).name() ); AnalysisObjectPtr append( ao ); int numEntries = 0; try { if( Histo1DPtr tmpHisto1D = dynamic_pointer_cast< YODA::Histo1D >( append ) ) { numEntries = tmpHisto1D->numEntries(); } else if ( Histo2DPtr tmpHisto2D = dynamic_pointer_cast< YODA::Histo2D >( append ) ) { numEntries = tmpHisto2D->numEntries(); } else if ( Profile1DPtr tmpProfile1D = dynamic_pointer_cast< YODA::Profile1D >( append ) ) { numEntries = tmpProfile1D->numEntries(); } else if ( Profile2DPtr tmpProfile2D = dynamic_pointer_cast< YODA::Profile2D >( append ) ) { numEntries = tmpProfile2D->numEntries(); } else if ( Scatter1DPtr tmpScatter1D = dynamic_pointer_cast< YODA::Scatter1D >( append ) ) { numEntries = tmpScatter1D->numPoints(); } else if ( Scatter2DPtr tmpScatter2D = dynamic_pointer_cast< YODA::Scatter2D >( append ) ) { numEntries = tmpScatter2D->numPoints(); } else if ( Scatter3DPtr tmpScatter3D = dynamic_pointer_cast< YODA::Scatter3D >( append ) ) { numEntries = tmpScatter3D->numPoints(); } if( numEntries <= 0 ) { continue; } if ( !_readObjects[ao->path()] ) { _readObjects[ao->path()] = append; } // in case an object with the given path is already read try to merge the new and the old one and save this version // else // merge...; or fail } catch(...) { MSG_ERROR( "Unable to read data!" ); throw UserError( "Unexpected error while reading data!" ); } } _haveReadData = true; } + string AnalysisHandler::runName() const { return _runname; } size_t AnalysisHandler::numEvents() const { return _numEvents; } double AnalysisHandler::sumOfWeights() const { return _sumOfWeights; } void AnalysisHandler::setSumOfWeights(const double& sum) { _sumOfWeights=sum; } 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) { _xs = xs; 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; } }