diff --git a/include/Rivet/Analysis.hh b/include/Rivet/Analysis.hh --- a/include/Rivet/Analysis.hh +++ b/include/Rivet/Analysis.hh @@ -1,1440 +1,1441 @@ // -*- 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/Cuts.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Tools/ParticleUtils.hh" #include "Rivet/Tools/BinnedHistogram.hh" #include "Rivet/Tools/RivetMT2.hh" #include "Rivet/Tools/RivetYODA.hh" #include "Rivet/Tools/Percentile.hh" #include "Rivet/Projections/CentralityProjection.hh" #include /// @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 { // Convenience for analysis writers using std::cout; using std::cerr; using std::endl; using std::tuple; using std::stringstream; using std::swap; using std::numeric_limits; // 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() { } //@} 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. If options has been set for this instance, a /// corresponding string is appended at the end. virtual std::string name() const { return ( (info().name().empty()) ? _defaultname : info().name() ) + _optstring; } /// Get name of reference data file, which could be different from plugin name virtual std::string getRefDataName() const { return (info().getRefDataName().empty()) ? _defaultname : info().getRefDataName(); } /// Set name of reference data file, which could be different from plugin name virtual void setRefDataName(const std::string& ref_data="") { info().setRefDataName(!ref_data.empty() ? ref_data : 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(); } /// make-style commands for validating this analysis. virtual std::vector validation() const { return info().validation(); } /// 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; } //@} /// @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; /// Check if we are running rivet-merge. bool merging() const { return sqrtS() <= 0.0; } //@} /// @name Analysis / beam compatibility testing /// @todo Replace with beamsCompatible() with no args (calling beams() function internally) /// @todo Add beamsMatch() methods with same (shared-code?) tolerance as in beamsCompatible() //@{ /// 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; //@} /// 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; /// @brief Get the number of events seen (via the analysis handler). /// /// @note Use in the finalize phase only. size_t numEvents() const; /// @brief Get the sum of event weights seen (via the analysis handler). /// /// @note Use in the finalize phase only. double sumW() const; /// Alias double sumOfWeights() const { return sumW(); } /// @brief Get the sum of squared event weights seen (via the analysis handler). /// /// @note Use in the finalize phase only. double sumW2() 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 mkAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const; //@} /// @name Histogram reference data //@{ /// Get reference data for a named histo /// @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 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 = mkAxisCode(datasetId, xAxisId, yAxisId); return refData(hname); } //@} /// @name Counter booking //@{ /// Book a counter. CounterPtr & book(CounterPtr &, 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 & book(CounterPtr &, 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 & book(Histo1DPtr &,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 & book(Histo1DPtr &,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 non-uniform bins defined by the vector of bin edges @a binedges . Histo1DPtr & book(Histo1DPtr &,const std::string& name, const std::initializer_list& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram with binning from a reference scatter. Histo1DPtr & book(Histo1DPtr &,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 & book(Histo1DPtr &,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 & book(Histo1DPtr &,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 & book(Histo2DPtr &,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 /// vectors of bin edges @a xbinedges and @a ybinedges. Histo2DPtr & book(Histo2DPtr &,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 non-uniform bins defined by the /// vectors of bin edges @a xbinedges and @a ybinedges. Histo2DPtr & book(Histo2DPtr &,const std::string& name, const std::initializer_list& xbinedges, const std::initializer_list& 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 & book(Histo2DPtr &,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 & book(Histo2DPtr &,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 & book(Histo2DPtr &,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 & book(Profile1DPtr &, 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 & book(Profile1DPtr &, 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 non-uniform bins defined by the vector of bin edges @a binedges . Profile1DPtr & book(Profile1DPtr &, const std::string& name, const std::initializer_list& 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 & book(Profile1DPtr &, 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 & book(Profile1DPtr &, 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 & book(Profile1DPtr &, 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 & book(Profile2DPtr &, 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 & book(Profile2DPtr &, 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 non-uniform bins defined by the vectorx /// of bin edges @a xbinedges and @a ybinedges. Profile2DPtr & book(Profile2DPtr &, const std::string& name, const std::initializer_list& xbinedges, const std::initializer_list& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// @todo REINSTATE // /// Book a 2D profile histogram with binning from a reference scatter. // Profile2DPtr & book(const Profile2DPtr &, 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 & book(const Profile2DPtr &, 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 & book(const Profile2DPtr &, 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 & book(Scatter2DPtr & s2d, const string& hname, 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 & book(Scatter2DPtr & s2d, 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 & book(Scatter2DPtr & s2d, const string& hname, 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 & book(Scatter2DPtr & s2d, const string& hname, const std::vector& binedges, const std::string& title, const std::string& xtitle, const std::string& ytitle); /// Book a 2-dimensional data point set with x-points from an existing scatter and a new path. Scatter2DPtr book(Scatter2DPtr & s2d, const Scatter2DPtr & scPtr, const std::string& path, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); //@} public: /// @name Accessing options for this Analysis instance. //@{ /// Return the map of all options given to this analysis. const std::map & options() { return _options; } /// Get an option for this analysis instance as a string. std::string getOption(std::string optname) { if ( _options.find(optname) != _options.end() ) return _options.find(optname)->second; return ""; } /// Get an option for this analysis instance converted to a /// specific type (given by the specified @a def value). template T getOption(std::string optname, T def) { if (_options.find(optname) == _options.end()) return def; std::stringstream ss; ss << _options.find(optname)->second; T ret; ss >> ret; return ret; } //@} /// @brief Book a CentralityProjection /// /// Using a SingleValueProjection, @a proj, giving the value of an /// experimental observable to be used as a centrality estimator, /// book a CentralityProjection based on the experimentally /// measured pecentiles of this observable (as given by the /// reference data for the @a calHistName histogram in the @a /// calAnaName analysis. If a preloaded file with the output of a /// run using the @a calAnaName analysis contains a valid /// generated @a calHistName histogram, it will be used as an /// optional percentile binning. Also if this preloaded file /// contains a histogram with the name @a calHistName with an /// appended "_IMP" This histogram will be used to add an optional /// centrality percentile based on the generated impact /// parameter. If @increasing is true, a low (high) value of @proj /// is assumed to correspond to a more peripheral (central) event. const CentralityProjection& declareCentrality(const SingleValueProjection &proj, string calAnaName, string calHistName, const string projName, bool increasing = false); /// @brief Book a Percentile wrapper around AnalysisObjects. /// /// Based on a previously registered CentralityProjection named @a /// projName book one AnalysisObject for each @a centralityBin and /// name them according to the corresponding code in the @a ref /// vector. /// /// @todo Convert to just be called book() cf. others template Percentile bookPercentile(string projName, vector > centralityBins, vector > ref) { typedef typename ReferenceTraits::RefT RefT; typedef rivet_shared_ptr> WrapT; Percentile pctl(this, projName); const int nCent = centralityBins.size(); for (int iCent = 0; iCent < nCent; ++iCent) { const string axisCode = mkAxisCode(std::get<0>(ref[iCent]), std::get<1>(ref[iCent]), std::get<2>(ref[iCent])); const RefT & refscatter = refData(axisCode); WrapT wtf(_weightNames(), T(refscatter, histoPath(axisCode))); wtf = addAnalysisObject(wtf); CounterPtr cnt(_weightNames(), Counter(histoPath("TMP/COUNTER/" + axisCode))); cnt = addAnalysisObject(cnt); pctl.add(wtf, cnt, centralityBins[iCent]); } return pctl; } // /// @brief Book Percentile wrappers around AnalysisObjects. // /// // /// Based on a previously registered CentralityProjection named @a // /// projName book one (or several) AnalysisObject(s) named // /// according to @a ref where the x-axis will be filled according // /// to the percentile output(s) of the @projName. // /// // /// @todo Convert to just be called book() cf. others // template // PercentileXaxis bookPercentileXaxis(string projName, // tuple ref) { // typedef typename ReferenceTraits::RefT RefT; // typedef rivet_shared_ptr> WrapT; // PercentileXaxis pctl(this, projName); // const string axisCode = mkAxisCode(std::get<0>(ref), // std::get<1>(ref), // std::get<2>(ref)); // const RefT & refscatter = refData(axisCode); // WrapT wtf(_weightNames(), T(refscatter, histoPath(axisCode))); // wtf = addAnalysisObject(wtf); // CounterPtr cnt(_weightNames(), Counter()); // cnt = addAnalysisObject(cnt); // pctl.add(wtf, cnt); // return pctl; // } private: // Functions that have to be defined in the .cc file to avoid circular #includes /// Get the list of weight names from the handler vector _weightNames() const; /// Get the list of weight names from the handler YODA::AnalysisObjectPtr _getPreload(string name) const; /// Get the default/nominal weight index size_t _defaultWeightIndex() const; /// Get an AO from another analysis MultiweightAOPtr _getOtherAnalysisObject(const std::string & ananame, const std::string& name); /// Check that analysis objects aren't being booked/registered outside the init stage void _checkBookInit() const; /// Check if we are in the init stage. bool inInit() const; /// Check if we are in the finalize stage. bool inFinalize() const; private: /// To be used in finalize context only: class CounterAdapter { public: CounterAdapter(double x) : x_(x ) {} CounterAdapter(const YODA::Counter & c) : x_(c.val() ) {} // CounterAdapter(CounterPtr cp) : x_(cp->val() ) {} CounterAdapter(const YODA::Scatter1D & s) : x_(s.points()[0].x()) { assert( s.numPoints() == 1 || "Can only scale by a single value."); } // CounterAdapter(Scatter1DPtr sp) : x_(sp->points()[0].x()) { // assert( sp->numPoints() == 1 || "Can only scale by a single value."); // } operator double() const { return x_; } private: double x_; }; public: double dbl(double x) { return x; } double dbl(const YODA::Counter & c) { return c.val(); } double dbl(const YODA::Scatter1D & s) { assert( s.numPoints() == 1 ); return s.points()[0].x(); } /// @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, CounterAdapter 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, CounterAdapter factor) { for (auto& c : cnts) scale(c, factor); } /// @todo YUCK! template void scale(const CounterPtr (&cnts)[array_size], CounterAdapter 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, CounterAdapter 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, CounterAdapter norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// @todo YUCK! template void normalize(const Histo1DPtr (&histos)[array_size], CounterAdapter 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, CounterAdapter 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, CounterAdapter factor) { for (auto& h : histos) scale(h, factor); } /// @todo YUCK! template void scale(const Histo1DPtr (&histos)[array_size], CounterAdapter factor) { for (auto& h : histos) scale(h, factor); } /// Normalize the given histogram, @a histo, to area = @a norm. void normalize(Histo2DPtr histo, CounterAdapter 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, CounterAdapter norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// @todo YUCK! template void normalize(const Histo2DPtr (&histos)[array_size], CounterAdapter 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, CounterAdapter 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, CounterAdapter factor) { for (auto& h : histos) scale(h, factor); } /// @todo YUCK! template void scale(const Histo2DPtr (&histos)[array_size], CounterAdapter 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 //@{ /// Get a preloaded YODA object. template shared_ptr getPreload(string path) const { return dynamic_pointer_cast(_getPreload(path)); } /// Register a new data object, optionally read in preloaded data. template rivet_shared_ptr< Wrapper > registerAO(const YODAT & yao) { typedef Wrapper WrapperT; typedef shared_ptr YODAPtrT; typedef rivet_shared_ptr RAOT; if ( !inInit() && !inFinalize() ) { MSG_ERROR("Can't book objects outside of init()"); throw UserError(name() + ": Can't book objects outside of init() or finalize()."); } // First check that we haven't booked this before. This is // allowed when booking in finalze. for (auto & waold : analysisObjects()) { if ( yao.path() == waold.get()->basePath() ) { if ( !inInit() ) MSG_WARNING("Found double-booking of " << yao.path() << " in " << name() << ". Keeping previous booking"); return RAOT(dynamic_pointer_cast(waold.get())); } } shared_ptr wao = make_shared(); + wao->_basePath = yao.path(); YODAPtrT yaop = make_shared(yao); for (const string& weightname : _weightNames()) { // Create two YODA objects for each weight. Copy from // preloaded YODAs if present. First the finalized yoda: string finalpath = yao.path(); if ( weightname != "" ) finalpath += "[" + weightname + "]"; YODAPtrT preload = getPreload(finalpath); if ( preload ) { if ( !bookingCompatible(preload, yaop) ) { MSG_WARNING("Found incompatible pre-existing data object with same base path " << finalpath << " for " << name()); preload = nullptr; } else { MSG_TRACE("Using preloaded " << finalpath << " in " <_final.push_back(make_shared(*preload)); } } if ( !preload ) { wao->_final.push_back(make_shared(yao)); wao->_final.back()->setPath(finalpath); } // Then the raw filling yodas. string rawpath = "/RAW" + finalpath; preload = getPreload(rawpath); if ( preload ) { if ( !bookingCompatible(preload, yaop) ) { MSG_WARNING("Found incompatible pre-existing data object with same base path " << rawpath << " for " << name()); preload = nullptr; } else { MSG_TRACE("Using preloaded " << rawpath << " in " <_persistent.push_back(make_shared(*preload)); } } if ( !preload ) { wao->_persistent.push_back(make_shared(yao)); wao->_persistent.back()->setPath(rawpath); } } rivet_shared_ptr ret(wao); ret.get()->unsetActiveWeight(); if ( inFinalize() ) { // If booked in finalize() we assume it is the first time // finalize is run. ret.get()->pushToFinal(); ret.get()->setActiveFinalWeightIdx(0); } _analysisobjects.push_back(ret); return ret; } /// Register a data object in the histogram system template AO addAnalysisObject(const AO & aonew) { _checkBookInit(); for (const MultiweightAOPtr& ao : analysisObjects()) { // Check AO base-name first ao.get()->setActiveWeightIdx(_defaultWeightIndex()); aonew.get()->setActiveWeightIdx(_defaultWeightIndex()); if (ao->path() != aonew->path()) continue; // If base-name matches, check compatibility // NB. This evil is because dynamic_ptr_cast can't work on rivet_shared_ptr directly AO aoold = AO(dynamic_pointer_cast(ao.get())); //< OMG if ( !aoold || !bookingCompatible(aonew, aoold) ) { MSG_WARNING("Found incompatible pre-existing data object with same base path " << aonew->path() << " for " << name()); throw LookupError("Found incompatible pre-existing data object with same base path during AO booking"); } // Finally, check all weight variations for (size_t weightIdx = 0; weightIdx < _weightNames().size(); ++weightIdx) { aoold.get()->setActiveWeightIdx(weightIdx); aonew.get()->setActiveWeightIdx(weightIdx); if (aoold->path() != aonew->path()) { MSG_WARNING("Found incompatible pre-existing data object with different weight-path " << aonew->path() << " for " << name()); throw LookupError("Found incompatible pre-existing data object with same weight-path during AO booking"); } } // They're fully compatible: bind and return aoold.get()->unsetActiveWeight(); MSG_TRACE("Bound pre-existing data object " << aoold->path() << " for " << name()); return aoold; } // No equivalent found MSG_TRACE("Registered " << aonew->annotation("Type") << " " << aonew->path() << " for " << name()); aonew.get()->unsetActiveWeight(); _analysisobjects.push_back(aonew); return aonew; } /// 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(const MultiweightAOPtr & ao); // /// Get all data objects, for all analyses, from the AnalysisHandler // /// @todo Can we remove this? Why not call handler().getData()? // vector getAllData(bool includeorphans) const; /// Get a Rivet data object from the histogram system template const AO getAnalysisObject(const std::string& aoname) const { for (const MultiweightAOPtr& ao : analysisObjects()) { ao.get()->setActiveWeightIdx(_defaultWeightIndex()); if (ao->path() == histoPath(aoname)) { // return dynamic_pointer_cast(ao); return AO(dynamic_pointer_cast(ao.get())); } } throw LookupError("Data object " + histoPath(aoname) + " not found"); } // /// Get a data object from the histogram system // 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 LookupError("Data object " + histoPath(name) + " not found"); // } // /// Get a data object from the histogram system (non-const) // 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 LookupError("Data object " + histoPath(name) + " not found"); // } /// Get a data object from another analysis (e.g. preloaded /// calibration histogram). template AO getAnalysisObject(const std::string& ananame, const std::string& aoname) { MultiweightAOPtr ao = _getOtherAnalysisObject(ananame, aoname); // return dynamic_pointer_cast(ao); return AO(dynamic_pointer_cast(ao.get())); } // /// 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; /// @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; /// Options the (this instance of) the analysis map _options; /// The string of options. string _optstring; 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/src/Core/Analysis.cc b/src/Core/Analysis.cc --- a/src/Core/Analysis.cc +++ b/src/Core/Analysis.cc @@ -1,1073 +1,1073 @@ // -*- 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) : _analysishandler(nullptr) { 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 { std::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; } /////////////////////////////////////////// double Analysis::crossSection() const { const YODA::Scatter1D::Points& ps = handler().crossSection()->points(); if (ps.size() != 1) { string errMsg = "cross section missing for analysis " + name(); throw Error(errMsg); } return ps[0].x(); } double Analysis::crossSectionPerEvent() const { 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, false, false); // } CounterPtr & Analysis::book(CounterPtr & ctr, const string& cname, const string& title) { // const string path = histoPath(cname); // ctr = CounterPtr(handler().weightNames(), Counter(path, title)); // ctr = addAnalysisObject(ctr); // return ctr; return ctr = registerAO(Counter(histoPath(cname), title)); } CounterPtr & Analysis::book(CounterPtr & ctr, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return book(ctr, axisCode, title); } Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname, size_t nbins, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Histo1D hist = Histo1D(nbins, lower, upper, path, title); hist.setAnnotation("XLabel", xtitle); hist.setAnnotation("YLabel", ytitle); // histo = Histo1DPtr(handler().weightNames(), hist); // histo = addAnalysisObject(histo); // return histo; return histo = registerAO(hist); } Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname, const initializer_list& binedges, const string& title, const string& xtitle, const string& ytitle) { return book(histo, hname, vector{binedges}, title, xtitle, ytitle); } Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Histo1D hist = Histo1D(binedges, path, title); hist.setAnnotation("XLabel", xtitle); hist.setAnnotation("YLabel", ytitle); // histo = Histo1DPtr(handler().weightNames(), hist); // histo = addAnalysisObject(histo); // return histo; return histo = registerAO(hist); } Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname, const string& title, const string& xtitle, const string& ytitle) { const Scatter2D& refdata = refData(hname); return book(histo, hname, refdata, title, xtitle, ytitle); } Histo1DPtr & Analysis::book(Histo1DPtr & histo, 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 book(histo, axisCode, title, xtitle, ytitle); } Histo1DPtr & Analysis::book(Histo1DPtr& histo, const string& hname, const Scatter2D& refscatter, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Histo1D hist = Histo1D(refscatter, path); hist.setTitle(title); hist.setAnnotation("XLabel", xtitle); hist.setAnnotation("YLabel", ytitle); if (hist.hasAnnotation("IsRef")) hist.rmAnnotation("IsRef"); // histo = Histo1DPtr(handler().weightNames(), hist); // histo = addAnalysisObject(histo); // return histo; return histo = registerAO(hist); } ///////////////// Histo2DPtr & Analysis::book(Histo2DPtr & h2d,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); Histo2D hist(nxbins, xlower, xupper, nybins, ylower, yupper, path, title); hist.setAnnotation("XLabel", xtitle); hist.setAnnotation("YLabel", ytitle); hist.setAnnotation("ZLabel", ztitle); // h2d = Histo2DPtr(handler().weightNames(), hist); // h2d = addAnalysisObject(h2d); // return h2d; return h2d = registerAO(hist); } Histo2DPtr & Analysis::book(Histo2DPtr & h2d,const string& hname, const initializer_list& xbinedges, const initializer_list& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { return book(h2d, hname, vector{xbinedges}, vector{ybinedges}, title, xtitle, ytitle, ztitle); } Histo2DPtr & Analysis::book(Histo2DPtr & h2d,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); Histo2D hist(xbinedges, ybinedges, path, title); hist.setAnnotation("XLabel", xtitle); hist.setAnnotation("YLabel", ytitle); hist.setAnnotation("ZLabel", ztitle); // h2d = Histo2DPtr(handler().weightNames(), hist); // h2d = addAnalysisObject(h2d); // return h2d; return h2d = registerAO(hist); } Histo2DPtr & Analysis::book(Histo2DPtr & histo, const string& hname, const Scatter3D& refscatter, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Histo2D hist = Histo2D(refscatter, path); hist.setTitle(title); hist.setAnnotation("XLabel", xtitle); hist.setAnnotation("YLabel", ytitle); hist.setAnnotation("ZLabel", ztitle); if (hist.hasAnnotation("IsRef")) hist.rmAnnotation("IsRef"); // histo = Histo2DPtr(handler().weightNames(), hist); // histo = addAnalysisObject(histo); // return histo; return histo = registerAO(hist); } Histo2DPtr & Analysis::book(Histo2DPtr & histo, const string& hname, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const Scatter3D& refdata = refData(hname); return book(histo, hname, refdata, title, xtitle, ytitle, ztitle); } Histo2DPtr & Analysis::book(Histo2DPtr & histo, 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 book(histo, axisCode, title, xtitle, ytitle, ztitle); } ///////////////// Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname, size_t nbins, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Profile1D prof(nbins, lower, upper, path, title); prof.setAnnotation("XLabel", xtitle); prof.setAnnotation("YLabel", ytitle); // p1d = Profile1DPtr(handler().weightNames(), prof); // p1d = addAnalysisObject(p1d); // return p1d; return p1d = registerAO(prof); } Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname, const initializer_list& binedges, const string& title, const string& xtitle, const string& ytitle) { return book(p1d, hname, vector{binedges}, title, xtitle, ytitle); } Profile1DPtr & Analysis::book(Profile1DPtr & p1d, const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Profile1D prof(binedges, path, title); prof.setAnnotation("XLabel", xtitle); prof.setAnnotation("YLabel", ytitle); // p1d = Profile1DPtr(handler().weightNames(), prof); // p1d = addAnalysisObject(p1d); // return p1d; return p1d = registerAO(prof); } Profile1DPtr & Analysis::book(Profile1DPtr & p1d, const string& hname, const Scatter2D& refscatter, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Profile1D prof(refscatter, path); prof.setTitle(title); prof.setAnnotation("XLabel", xtitle); prof.setAnnotation("YLabel", ytitle); if (prof.hasAnnotation("IsRef")) prof.rmAnnotation("IsRef"); // p1d = Profile1DPtr(handler().weightNames(), prof); // p1d = addAnalysisObject(p1d); // return p1d; return p1d = registerAO(prof); } Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname, const string& title, const string& xtitle, const string& ytitle) { const Scatter2D& refdata = refData(hname); book(p1d, hname, refdata, title, xtitle, ytitle); return p1d; } Profile1DPtr & Analysis::book(Profile1DPtr & p1d,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 book(p1d, axisCode, title, xtitle, ytitle); } /////////////////// Profile2DPtr & Analysis::book(Profile2DPtr & p2d, 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); Profile2D prof(nxbins, xlower, xupper, nybins, ylower, yupper, path, title); prof.setAnnotation("XLabel", xtitle); prof.setAnnotation("YLabel", ytitle); prof.setAnnotation("ZLabel", ztitle); // p2d = Profile2DPtr(handler().weightNames(), prof); // p2d = addAnalysisObject(p2d); // return p2d; return p2d = registerAO(prof); } Profile2DPtr & Analysis::book(Profile2DPtr & p2d, const string& hname, const initializer_list& xbinedges, const initializer_list& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { return book(p2d, hname, vector{xbinedges}, vector{ybinedges}, title, xtitle, ytitle, ztitle); } Profile2DPtr & Analysis::book(Profile2DPtr & p2d, 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); Profile2D prof(xbinedges, ybinedges, path, title); prof.setAnnotation("XLabel", xtitle); prof.setAnnotation("YLabel", ytitle); prof.setAnnotation("ZLabel", ztitle); // p2d = Profile2DPtr(handler().weightNames(), prof); // p2d = addAnalysisObject(p2d); // return p2d; return p2d = registerAO(prof); } /// @todo REINSTATE // Profile2DPtr Analysis::book(Profile2DPtr& prof,const string& hname, // const Scatter3D& refscatter, // const string& title, // const string& xtitle, // const string& ytitle, // const string& ztitle) { // const string path = histoPath(hname); // /// @todo Add no-metadata argument to YODA copy constructors // Profile2D prof(refscatter, path); // prof.setTitle(title); // prof.setAnnotation("XLabel", xtitle); // prof.setAnnotation("YLabel", ytitle); // prof.setAnnotation("ZLabel", ztitle); // if (prof.hasAnnotation("IsRef")) prof.rmAnnotation("IsRef"); // p2d = Profile2DPtr(handler().weightNames(), prof); // p2d = addAnalysisObject(p2d); // return p2d; // } // Profile2DPtr Analysis::book(Profile2DPtr& prof, const string& hname, // const string& title, // const string& xtitle, // const string& ytitle, // const string& ztitle) { // const Scatter3D& refdata = refData(hname); // return book(prof, hname, refdata, title, xtitle, ytitle, ztitle); // } /////////////// /// @todo Should be able to book Scatter1Ds /////////////// Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, 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 book(s2d, axisCode, copy_pts, title, xtitle, ytitle); } Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, const string& hname, bool copy_pts, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Scatter2D scat; if (copy_pts) { const Scatter2D& refdata = refData(hname); scat = Scatter2D(refdata, path); for (Point2D& p : scat.points()) p.setY(0, 0); } else { scat = Scatter2D(path); } scat.setTitle(title); scat.setAnnotation("XLabel", xtitle); scat.setAnnotation("YLabel", ytitle); if (scat.hasAnnotation("IsRef")) scat.rmAnnotation("IsRef"); // s2d = Scatter2DPtr(handler().weightNames(), scat); // s2d = addAnalysisObject(s2d); // return s2d; return s2d = registerAO(scat); } Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, const string& hname, size_t npts, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Scatter2D scat; const double binwidth = (upper-lower)/npts; for (size_t pt = 0; pt < npts; ++pt) { const double bincentre = lower + (pt + 0.5) * binwidth; scat.addPoint(bincentre, 0, binwidth/2.0, 0); } scat.setTitle(title); scat.setAnnotation("XLabel", xtitle); scat.setAnnotation("YLabel", ytitle); // s2d = Scatter2DPtr(handler().weightNames(), scat); // s2d = addAnalysisObject(s2d); // return s2d; return s2d = registerAO(scat); } Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Scatter2D scat; 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]; scat.addPoint(bincentre, 0, binwidth/2.0, 0); } scat.setTitle(title); scat.setAnnotation("XLabel", xtitle); scat.setAnnotation("YLabel", ytitle); // s2d = Scatter2DPtr(handler().weightNames(), scat); // s2d = addAnalysisObject(s2d); // return s2d; return s2d = registerAO(scat); } /// @todo Book Scatter3Ds? ///////////////////// 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, Analysis::CounterAdapter factor) { if (!cnt) { MSG_WARNING("Failed to scale counter=NULL in analysis " << name() << " (scale=" << double(factor) << ")"); return; } if (std::isnan(double(factor)) || std::isinf(double(factor))) { MSG_WARNING("Failed to scale counter=" << cnt->path() << " in analysis: " << name() << " (invalid scale factor = " << double(factor) << ")"); factor = 0; } MSG_TRACE("Scaling counter " << cnt->path() << " by factor " << double(factor)); try { cnt->scaleW(factor); } catch (YODA::Exception& we) { MSG_WARNING("Could not scale counter " << cnt->path()); return; } } void Analysis::normalize(Histo1DPtr histo, Analysis::CounterAdapter norm, bool includeoverflows) { if (!histo) { MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")"); return; } MSG_TRACE("Normalizing histo " << histo->path() << " to " << double(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, Analysis::CounterAdapter factor) { if (!histo) { MSG_WARNING("Failed to scale histo=NULL in analysis " << name() << " (scale=" << double(factor) << ")"); return; } if (std::isnan(double(factor)) || std::isinf(double(factor))) { MSG_WARNING("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << double(factor) << ")"); factor = 0; } MSG_TRACE("Scaling histo " << histo->path() << " by factor " << double(factor)); try { histo->scaleW(factor); } catch (YODA::Exception& we) { MSG_WARNING("Could not scale histo " << histo->path()); return; } } void Analysis::normalize(Histo2DPtr histo, Analysis::CounterAdapter norm, bool includeoverflows) { if (!histo) { MSG_ERROR("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")"); return; } MSG_TRACE("Normalizing histo " << histo->path() << " to " << double(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, Analysis::CounterAdapter factor) { if (!histo) { MSG_ERROR("Failed to scale histo=NULL in analysis " << name() << " (scale=" << double(factor) << ")"); return; } if (std::isnan(double(factor)) || std::isinf(double(factor))) { MSG_ERROR("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << double(factor) << ")"); factor = 0; } MSG_TRACE("Scaling histo " << histo->path() << " by factor " << double(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?!? ////////////////////////////////// // namespace { // void errormsg(std::string name) { // // #ifdef HAVE_BACKTRACE // // void * buffer[4]; // // backtrace(buffer, 4); // // backtrace_symbols_fd(buffer, 4 , 1); // // #endif // std::cerr << name << ": Can't book objects outside of init().\n"; // assert(false); // } // } namespace Rivet { // void Analysis::addAnalysisObject(const MultiweightAOPtr & ao) { // if (handler().stage() == AnalysisHandler::Stage::INIT) { // _analysisobjects.push_back(ao); // } // else { // errormsg(name()); // } // } void Analysis::removeAnalysisObject(const string& path) { for (auto it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) { if ((*it)->path() == path) { _analysisobjects.erase(it); break; } } } void Analysis::removeAnalysisObject(const MultiweightAOPtr & ao) { for (auto 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" ) { YODA::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" ) { YODA::Histo1DPtr genhists = getPreload("/" + calAnaName + "/" + calHistName); // for ( YODA::AnalysisObjectPtr ao : handler().getData(true) ) { // if ( ao->path() == histpath ) // genhist = dynamic_pointer_cast(ao); // } if ( !genhists || genhists->numEntries() <= 1 ) { MSG_WARNING("No generated calibration histogram for " << "CentralityProjection " << projName << " found " << "(requested histogram " << calHistName << " in " << calAnaName << ")"); } else { MSG_INFO("Found calibration histogram " << sel << " " << genhists->path()); cproj.add(PercentileProjection(proj, *genhists, increasing), sel); } } else if ( sel == "IMP" ) { YODA::Histo1DPtr imphists = getPreload("/" + calAnaName + "/" + calHistName + "_IMP"); if ( !imphists || imphists->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 << " " << imphists->path()); cproj.add(PercentileProjection(ImpactParameterProjection(), *imphists, true), sel); } } else if ( sel == "USR" ) { #if HEPMC_VERSION_CODE >= 3000000 YODA::Histo1DPtr usrhists = getPreload("/" + calAnaName + "/" + calHistName + "_USR"); if ( !usrhists || usrhists->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 << " " << usrhists->path()); cproj.add((UserCentEstimate(), usrhists*, 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); } vector Analysis::_weightNames() const { return handler().weightNames(); } YODA::AnalysisObjectPtr Analysis::_getPreload(string path) const { return handler().getPreload(path); } size_t Analysis::_defaultWeightIndex() const { return handler().defaultWeightIndex(); } MultiweightAOPtr Analysis::_getOtherAnalysisObject(const std::string & ananame, const std::string& name) { std::string path = "/" + ananame + "/" + name; const auto& ana = handler().analysis(ananame); return ana->getAnalysisObject(name); //< @todo includeorphans check?? } void Analysis::_checkBookInit() const { if (handler().stage() != AnalysisHandler::Stage::INIT) { MSG_ERROR("Can't book objects outside of init()"); throw UserError(name() + ": Can't book objects outside of init()."); } } bool Analysis::inInit() const { - return handler().stage() != AnalysisHandler::Stage::INIT; + return handler().stage() == AnalysisHandler::Stage::INIT; } bool Analysis::inFinalize() const { - return handler().stage() != AnalysisHandler::Stage::FINALIZE; + return handler().stage() == AnalysisHandler::Stage::FINALIZE; } }