Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-04-19 09:06:56

0001 // -*- C++ -*-
0002 #ifndef RIVET_Analysis_HH
0003 #define RIVET_Analysis_HH
0004 
0005 #include "Rivet/Config/RivetCommon.hh"
0006 #include "Rivet/AnalysisInfo.hh"
0007 #include "Rivet/Event.hh"
0008 #include "Rivet/Projection.hh"
0009 #include "Rivet/ProjectionApplier.hh"
0010 #include "Rivet/ProjectionHandler.hh"
0011 #include "Rivet/AnalysisLoader.hh"
0012 #include "Rivet/Tools/Cuts.hh"
0013 #include "Rivet/Tools/Logging.hh"
0014 #include "Rivet/Tools/ParticleUtils.hh"
0015 #ifdef HAVE_H5
0016 #include "Rivet/Tools/RivetHDF5.hh"
0017 #endif
0018 #include "Rivet/Tools/HistoGroup.hh"
0019 #include "Rivet/Tools/RivetMT2.hh"
0020 #include "Rivet/Tools/RivetPaths.hh"
0021 #include "Rivet/Tools/RivetYODA.hh"
0022 #include "Rivet/Tools/Percentile.hh"
0023 #include "Rivet/Tools/Cutflow.hh"
0024 #include "Rivet/Projections/CentralityProjection.hh"
0025 #include <tuple>
0026 
0027 
0028 /// @def vetoEvent
0029 /// Preprocessor define for vetoing events, including the log message and return.
0030 #define vetoEvent                                                       \
0031   do { MSG_DEBUG("Vetoing event on line " << __LINE__ << " of " << __FILE__); return; } while(0)
0032 
0033 
0034 namespace Rivet {
0035 
0036 
0037   // Convenience for analysis writers
0038   using std::cout;
0039   using std::cerr;
0040   using std::endl;
0041   using std::tuple;
0042   using std::stringstream;
0043   using std::swap;
0044   using std::numeric_limits;
0045 
0046 
0047   // Forward declaration
0048   class AnalysisHandler;
0049 
0050 
0051   /// @brief This is the base class of all analysis classes in Rivet.
0052   ///
0053   /// There are
0054   /// three virtual functions which should be implemented in base classes:
0055   ///
0056   /// void init() is called by Rivet before a run is started. Here the
0057   /// analysis class should book necessary histograms. The needed
0058   /// projections should probably rather be constructed in the
0059   /// constructor.
0060   ///
0061   /// void analyze(const Event&) is called once for each event. Here the
0062   /// analysis class should apply the necessary Projections and fill the
0063   /// histograms.
0064   ///
0065   /// void finalize() is called after a run is finished. Here the analysis
0066   /// class should do whatever manipulations are necessary on the
0067   /// histograms. Writing the histograms to a file is, however, done by
0068   /// the AnalysisHandler class.
0069   class Analysis : public ProjectionApplier {
0070   public:
0071 
0072     /// The AnalysisHandler is a friend.
0073     friend class AnalysisHandler;
0074 
0075 
0076     /// Constructor
0077     Analysis(const std::string& name);
0078 
0079     /// The destructor
0080     virtual ~Analysis() {}
0081 
0082     /// The assignment operator is private and must be deleted, so it can never be called.
0083     Analysis& operator = (const Analysis&) = delete;
0084 
0085 
0086   public:
0087 
0088     /// @name Main analysis methods
0089     /// @{
0090 
0091     /// Initialize this analysis object. A concrete class should here
0092     /// book all necessary histograms. An overridden function must make
0093     /// sure it first calls the base class function.
0094     virtual void init() { }
0095 
0096     /// Analyze one event. A concrete class should here apply the
0097     /// necessary projections on the \a event and fill the relevant
0098     /// histograms. An overridden function must make sure it first calls
0099     /// the base class function.
0100     virtual void analyze(const Event& event) = 0;
0101 
0102     /// Finalize this analysis object. A concrete class should here make
0103     /// all necessary operations on the histograms. Writing the
0104     /// histograms to a file is, however, done by the Rivet class. An
0105     /// overridden function must make sure it first calls the base class
0106     /// function.
0107     virtual void finalize() { }
0108 
0109     ///Call the projection applier _syncDeclQueue() method.
0110     ///(It should be hidden for all projection appliers other than analyses)
0111     ///TODO for TP: Is this the right block for this method to be in?
0112     void syncDeclQueue() {
0113       this->_syncDeclQueue();
0114       this->markAsOwned();
0115     }
0116 
0117     /// @}
0118 
0119 
0120   public:
0121 
0122     /// @name Metadata
0123     ///
0124     /// Metadata is used for querying from the command line and also for
0125     /// building web pages and the analysis pages in the Rivet manual.
0126     /// @{
0127 
0128     /// Get the AnalysisInfo object to parse its info file in which the metadata is stored.
0129     void loadInfo() { info().parseInfoFile(); }
0130 
0131     /// Get the actual AnalysisInfo object in which all this metadata is stored.
0132     const AnalysisInfo& info() const {
0133       if (!_info) throw Error("No AnalysisInfo object :-O");
0134       return *_info;
0135     }
0136 
0137     /// @brief Get the name of the analysis.
0138     ///
0139     /// By default this is computed by combining the results of the
0140     /// experiment, year and Spires ID metadata methods and you should
0141     /// only override it if there's a good reason why those won't
0142     /// work. If options has been set for this instance, a
0143     /// corresponding string is appended at the end.
0144     virtual std::string name() const {
0145       return ( (info().name().empty()) ? _defaultname : info().name() ) + _optstring;
0146     }
0147 
0148     /// @brief Get the path to a data file associated with this analysis
0149     ///
0150     /// The searched-for filename will be \<ANANAME\>.\<extn\> of suffix is empty/unspecified,
0151     /// or \<ANANAME\>-\<suffix\>.\<extn\> if a non-zero suffix is specified.
0152     std::string analysisDataPath(const std::string& extn, const std::string& suffix="") {
0153       string filename = name() + (suffix.empty() ? "" : "-") + suffix + "." + extn;
0154       return findAnalysisDataFile(filename);
0155     }
0156 
0157     /// Get the Inspire ID code for this analysis.
0158     virtual std::string inspireID() const {
0159       return info().inspireID();
0160     }
0161 
0162     /// Get the SPIRES ID code for this analysis (~deprecated).
0163     virtual std::string spiresID() const {
0164       return info().spiresID();
0165     }
0166 
0167     /// @brief Names & emails of paper/analysis authors.
0168     ///
0169     /// Names and email of authors in 'NAME \<EMAIL\>' format. The first
0170     /// name in the list should be the primary contact person.
0171     virtual std::vector<std::string> authors() const {
0172       return info().authors();
0173     }
0174 
0175     /// @brief Get a short description of the analysis.
0176     ///
0177     /// Short (one sentence) description used as an index entry.
0178     /// Use @a description() to provide full descriptive paragraphs
0179     /// of analysis details.
0180     virtual std::string summary() const {
0181       return info().summary();
0182     }
0183 
0184     /// @brief Get a full description of the analysis.
0185     ///
0186     /// Full textual description of this analysis, what it is useful for,
0187     /// what experimental techniques are applied, etc. Should be treated
0188     /// as a chunk of restructuredText (http://docutils.sourceforge.net/rst.html),
0189     /// with equations to be rendered as LaTeX with amsmath operators.
0190     virtual std::string description() const {
0191       return info().description();
0192     }
0193 
0194     /// @brief Information about the events needed as input for this analysis.
0195     ///
0196     /// Event types, energies, kinematic cuts, particles to be considered
0197     /// stable, etc. etc. Should be treated as a restructuredText bullet list
0198     /// (http://docutils.sourceforge.net/rst.html)
0199     virtual std::string runInfo() const {
0200       return info().runInfo();
0201     }
0202 
0203     /// Experiment which performed and published this analysis.
0204     virtual std::string experiment() const {
0205       return info().experiment();
0206     }
0207 
0208     /// Collider on which the experiment ran.
0209     virtual std::string collider() const {
0210       return info().collider();
0211     }
0212 
0213     /// When the original experimental analysis was published.
0214     virtual std::string year() const {
0215       return info().year();
0216     }
0217 
0218     /// The integrated luminosity in inverse femtobarn
0219     virtual double luminosityfb() const {
0220       return info().luminosityfb();
0221     }
0222     /// The integrated luminosity in inverse picobarn
0223     virtual double luminosity() const {
0224       return info().luminosity();
0225     }
0226     /// The integrated luminosity in inverse picobarn
0227     double luminositypb() const { return luminosity(); }
0228 
0229     /// Journal, and preprint references.
0230     virtual std::vector<std::string> references() const {
0231       return info().references();
0232     }
0233 
0234     /// BibTeX citation key for this article.
0235     virtual std::string bibKey() const {
0236       return info().bibKey();
0237     }
0238 
0239     /// BibTeX citation entry for this article.
0240     virtual std::string bibTeX() const {
0241       return info().bibTeX();
0242     }
0243 
0244     /// Whether this analysis is trusted (in any way!)
0245     virtual std::string status() const {
0246       return (info().status().empty()) ? "UNVALIDATED" : info().status();
0247     }
0248 
0249     /// A warning message from the info file, if there is one
0250     virtual std::string warning() const {
0251       return info().warning();
0252     }
0253 
0254     /// Any work to be done on this analysis.
0255     virtual std::vector<std::string> todos() const {
0256       return info().todos();
0257     }
0258 
0259     /// make-style commands for validating this analysis.
0260     virtual std::vector<std::string> validation() const {
0261       return info().validation();
0262     }
0263 
0264     /// Does this analysis have a reentrant finalize()?
0265     virtual bool reentrant() const {
0266       return info().reentrant();
0267     }
0268 
0269     /// Get vector of analysis keywords
0270     virtual const std::vector<std::string>& keywords() const {
0271       return info().keywords();
0272     }
0273 
0274     /// Positive filtering regex for ref-data HepData sync
0275     virtual std::string refMatch() const {
0276       return info().refMatch();
0277     }
0278 
0279     /// Negative filtering regex for ref-data HepData sync
0280     virtual std::string refUnmatch() const {
0281       return info().refUnmatch();
0282     }
0283 
0284     /// Positive filtering regex for setting double precision in Writer
0285     virtual std::string writerDoublePrecision() const {
0286       return info().writerDoublePrecision();
0287     }
0288 
0289     /// Return the allowed pairs of incoming beams required by this analysis.
0290     virtual const std::vector<PdgIdPair>& requiredBeamIDs() const {
0291       return info().beamIDs();
0292     }
0293     /// Declare the allowed pairs of incoming beams required by this analysis.
0294     virtual Analysis& setRequiredBeamIDs(const std::vector<PdgIdPair>& beamids) {
0295       info().setBeamIDs(beamids);
0296       return *this;
0297     }
0298 
0299     /// Sets of valid beam energy pairs, in GeV
0300     virtual const std::vector<std::pair<double, double> >& requiredBeamEnergies() const {
0301       return info().energies();
0302     }
0303     /// Declare the list of valid beam energy pairs, in GeV
0304     virtual Analysis& setRequiredBeamEnergies(const std::vector<std::pair<double, double> >& energies) {
0305       info().setEnergies(energies);
0306       return *this;
0307     }
0308 
0309 
0310     /// Location of reference data YODA file
0311     virtual std::string refFile() const {
0312       return info().refFile();
0313     }
0314     /// Get name of reference data file, which could be different from plugin name
0315     virtual std::string refDataName() const {
0316       return (info().getRefDataName().empty()) ? _defaultname : info().getRefDataName();
0317     }
0318     /// Set name of reference data file, which could be different from plugin name
0319     virtual void setRefDataName(const std::string& ref_data="") {
0320       info().setRefDataName(!ref_data.empty() ? ref_data : name());
0321     }
0322 
0323 
0324     /// @brief Get the actual AnalysisInfo object in which all this metadata is stored (non-const).
0325     ///
0326     /// @note For *internal* use!
0327     AnalysisInfo& info() {
0328       if (!_info) throw Error("No AnalysisInfo object :-O");
0329       return *_info;
0330     }
0331 
0332     /// @}
0333 
0334 
0335     /// @name Run conditions
0336     /// @{
0337 
0338     /// Check if we are running rivet-merge
0339     bool merging() const {
0340       return sqrtS() <= 0.0;
0341     }
0342 
0343     /// Check if the given conditions are compatible with this analysis' declared constraints
0344     bool compatibleWithRun() const;
0345 
0346     /// @}
0347 
0348 
0349     /// @name Analysis / beam compatibility testing
0350     ///
0351     /// @{
0352 
0353     /// Incoming beams for this run
0354     const ParticlePair& beams() const;
0355 
0356     /// Incoming beam IDs for this run
0357     PdgIdPair beamIDs() const;
0358 
0359     /// Incoming beam energies for this run
0360     pair<double,double> beamEnergies() const;
0361 
0362     /// Allowed centre-of-mass energies (in GeV) for this routine
0363     vector<double> allowedEnergies() const;
0364 
0365     /// Centre of mass energy for this run
0366     double sqrtS() const;
0367 
0368     /// Check if analysis is compatible with the provided beam particle IDs and energies
0369     bool beamsMatch(const ParticlePair& beams) const;
0370 
0371     /// Check if analysis is compatible with the provided beam particle IDs and energies
0372     bool beamsMatch(PdgId beam1, PdgId beam2, double e1, double e2) const;
0373 
0374     /// Check if analysis is compatible with the provided beam particle IDs and energies in GeV
0375     bool beamsMatch(const PdgIdPair& beams, const std::pair<double,double>& energies) const;
0376 
0377     /// Check if analysis is compatible with the provided beam particle IDs
0378     bool beamIDsMatch(PdgId beam1, PdgId beam2) const;
0379 
0380     /// Check if analysis is compatible with the provided beam particle IDs
0381     bool beamIDsMatch(const PdgIdPair& beamids) const;
0382 
0383     /// Check if analysis is compatible with the provided beam energies
0384     bool beamEnergiesMatch(double e1, double e2) const;
0385 
0386     /// Check if analysis is compatible with the provided beam energies
0387     bool beamEnergiesMatch(const std::pair<double,double>& energies) const;
0388 
0389     /// Check if analysis is compatible with the provided CoM energy
0390     bool beamEnergyMatch(const std::pair<double,double>& energies) const;
0391 
0392     /// Check if analysis is compatible with the provided CoM energy
0393     bool beamEnergyMatch(double sqrts) const;
0394 
0395     /// Check if sqrtS is compatible with provided value
0396     bool isCompatibleWithSqrtS(double energy, double tolerance=1e-5) const;
0397 
0398     /// @}
0399 
0400 
0401     /// Access the controlling AnalysisHandler object.
0402     AnalysisHandler& handler() const { return *_analysishandler; }
0403 
0404 
0405   protected:
0406 
0407     /// Get a Log object based on the name() property of the calling analysis object.
0408     Log& getLog() const;
0409 
0410     /// Get the process cross-section in pb. Throws if this hasn't been set.
0411     double crossSection() const;
0412 
0413     /// Get the process cross-section per generated event in pb. Throws if this
0414     /// hasn't been set.
0415     double crossSectionPerEvent() const;
0416 
0417     /// Get the process cross-section error in pb. Throws if this hasn't been set.
0418     double crossSectionError() const;
0419 
0420     /// Get the process cross-section error per generated event in
0421     /// pb. Throws if this hasn't been set.
0422     double crossSectionErrorPerEvent() const;
0423 
0424     /// @brief Get the number of events seen (via the analysis handler).
0425     ///
0426     /// @note Use in the finalize phase only.
0427     size_t numEvents() const;
0428 
0429     /// @brief Get the sum of event weights seen (via the analysis handler).
0430     ///
0431     /// @note Use in the finalize phase only.
0432     double sumW() const;
0433     /// Alias
0434     double sumOfWeights() const { return sumW(); }
0435 
0436     /// @brief Get the sum of squared event weights seen (via the analysis handler).
0437     ///
0438     /// @note Use in the finalize phase only.
0439     double sumW2() const;
0440 
0441 
0442   protected:
0443 
0444     /// @defgroup analysis_histopaths Histogram paths
0445     ///
0446     /// @todo Add "tmp" flags to return /ANA/TMP/foo/bar paths
0447     ///
0448     /// @{
0449 
0450     /// Get the canonical histogram "directory" path for this analysis.
0451     const std::string histoDir() const;
0452 
0453     /// Get the canonical histogram path for the named histogram in this analysis.
0454     const std::string histoPath(const std::string& hname) const;
0455 
0456     /// Get the canonical histogram path for the numbered histogram in this analysis.
0457     const std::string histoPath(unsigned int datasetID, unsigned int xAxisID, unsigned int yAxisID) const;
0458 
0459     /// Get the internal histogram name for given d, x and y (cf. HepData)
0460     const std::string mkAxisCode(unsigned int datasetID, unsigned int xAxisID, unsigned int yAxisID) const;
0461 
0462     /// @}
0463 
0464 
0465     #ifdef HAVE_H5
0466     /// @name Auxiliary HDF5 reference data
0467     /// @{
0468 
0469     /// Read in an aux data HDF5 file
0470     H5::File auxFile() const {
0471       return H5::readFile(name()+".h5");
0472     }
0473 
0474     /// Read HDF5 file @a filename
0475     template <typename T>
0476     bool auxData(const string& dsname, T& rtndata) {
0477       /// @todo Cache loading the H5 files? On MPI?
0478       return H5::readData(name()+".h5", dsname, rtndata);
0479     }
0480 
0481     template <typename T>
0482     T auxData(const string& dsname) {
0483       /// @todo Cache loading the H5 files? On MPI?
0484       return H5::readData<T>(name()+".h5", dsname);
0485     }
0486 
0487     /// @}
0488     #endif
0489 
0490 
0491     /// @name Histogram reference data
0492     /// @{
0493 
0494     /// Get all reference data objects for this analysis
0495     const std::map<std::string, YODA::AnalysisObjectPtr>& refData() const {
0496       _cacheRefData();
0497       return _refdata;
0498     }
0499 
0500 
0501     /// Get reference data for a named histo
0502     /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject?
0503     template <typename T=YODA::Estimate1D>
0504     const T& refData(const string& hname) const {
0505       _cacheRefData();
0506       MSG_TRACE("Using histo bin edges for " << name() << ":" << hname);
0507       if (!_refdata[hname]) {
0508         MSG_ERROR("Can't find reference histogram " << hname);
0509         throw Exception("Reference data " + hname + " not found.");
0510       }
0511       return dynamic_cast<T&>(*_refdata[hname]);
0512     }
0513 
0514 
0515     /// Get reference data for a numbered histo
0516     /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject?
0517     template <typename T=YODA::Estimate1D>
0518     const T& refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
0519       const string hname = mkAxisCode(datasetId, xAxisId, yAxisId);
0520       return refData<T>(hname);
0521     }
0522 
0523     /// @}
0524 
0525 
0526     /// @defgroup analysis_cbook Counter booking
0527     ///
0528     /// @todo Add "tmp" flags to book in standard temporary paths
0529     ///
0530     /// @{
0531 
0532     /// Book a counter.
0533     CounterPtr& book(CounterPtr&, const std::string& name);
0534 
0535     /// Book a counter, using a path generated from the dataset and axis ID codes
0536     ///
0537     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
0538     CounterPtr& book(CounterPtr&, unsigned int datasetID, unsigned int xAxisID, unsigned int yAxisID);
0539 
0540     /// @}
0541 
0542 
0543     /// @name Estimate booking
0544     /// @{
0545 
0546     /// Book an estimate.
0547     Estimate0DPtr& book(Estimate0DPtr&, const std::string& name);
0548 
0549     /// Book an estimate, using a path generated from the dataset and axis ID codes
0550     ///
0551     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
0552     Estimate0DPtr& book(Estimate0DPtr&, unsigned int datasetID, unsigned int xAxisID, unsigned int yAxisID);
0553 
0554 
0555     /// @}
0556 
0557     /// @name BinnedDbn booking
0558     /// @{
0559 
0560     /// Book a ND histogram with @a nbins uniformly distributed across the range @a lower - @a upper .
0561     template<size_t DbnN, typename... AxisT, typename = YODA::enable_if_all_CAxisT<AxisT...>>
0562     BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao,
0563                                        const std::string& name, const std::vector<size_t>& nbins,
0564                                        const std::vector<std::pair<double,double>>& loUpPairs) {
0565       if (nbins.size() != loUpPairs.size()) {
0566         throw RangeError("Vectors should have the same size!");
0567       }
0568       const string path = histoPath(name);
0569 
0570       YODA::BinnedDbn<DbnN, AxisT...> yao(nbins, loUpPairs, path);
0571       _setWriterPrecision(path, yao);
0572 
0573       return ao = registerAO(yao);
0574     }
0575 
0576     // Specialiation for 1D
0577     Histo1DPtr& book(Histo1DPtr& ao, const std::string& name,
0578                                      const size_t nbins, const double lower, const double upper) {
0579       return book(ao, name, vector<size_t>{nbins},
0580                   vector<pair<double,double>>{{lower,upper}});
0581     }
0582     //
0583     Profile1DPtr& book(Profile1DPtr& ao, const std::string& name,
0584                                          const size_t nbins, const double lower, const double upper) {
0585       return book(ao, name, vector<size_t>{nbins},
0586                   vector<pair<double,double>>{{lower,upper}});
0587     }
0588 
0589     // Specialiation for 2D
0590     Histo2DPtr& book(Histo2DPtr& ao, const std::string& name,
0591                                      const size_t nbinsX, const double lowerX, const double upperX,
0592                                      const size_t nbinsY, const double lowerY, const double upperY) {
0593       return book(ao, name, vector<size_t>{nbinsX,nbinsY},
0594                   vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}});
0595     }
0596     //
0597     Profile2DPtr& book(Profile2DPtr& ao, const std::string& name,
0598                                          const size_t nbinsX, const double lowerX, const double upperX,
0599                                          const size_t nbinsY, const double lowerY, const double upperY) {
0600       return book(ao, name, vector<size_t>{nbinsX,nbinsY},
0601                   vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}});
0602     }
0603 
0604     // Specialiation for 3D
0605     Histo3DPtr& book(Histo3DPtr& ao, const std::string& name,
0606                                      const size_t nbinsX, const double lowerX, const double upperX,
0607                                      const size_t nbinsY, const double lowerY, const double upperY,
0608                                      const size_t nbinsZ, const double lowerZ, const double upperZ) {
0609       return book(ao, name, vector<size_t>{nbinsX,nbinsY,nbinsZ},
0610                   vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}, {lowerZ,upperZ}});
0611     }
0612     //
0613     Profile3DPtr& book(Profile3DPtr& ao, const std::string& name,
0614                                          const size_t nbinsX, const double lowerX, const double upperX,
0615                                          const size_t nbinsY, const double lowerY, const double upperY,
0616                                          const size_t nbinsZ, const double lowerZ, const double upperZ) {
0617       return book(ao, name, vector<size_t>{nbinsX,nbinsY,nbinsZ},
0618                   vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}, {lowerZ,upperZ}});
0619     }
0620 
0621     /// Book a ND histogram with non-uniform bins defined by the vector of bin edges @a binedges .
0622     template<size_t DbnN, typename... AxisT>
0623     BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const std::string& name,
0624                                        const std::vector<AxisT>&... binedges) {
0625       const string path = histoPath(name);
0626       YODA::BinnedDbn<DbnN, AxisT...> yao(binedges..., path);
0627       _setWriterPrecision(path, yao);
0628 
0629       return ao = registerAO(yao);
0630     }
0631 
0632     /// Book a ND histogram with non-uniform bins defined by the vector of bin edges @a binedges .
0633     template<size_t DbnN, typename... AxisT>
0634     BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const std::string& name,
0635                                        const std::initializer_list<AxisT>&... binedges) {
0636       return book(ao, name, vector<AxisT>{binedges} ...);
0637     }
0638 
0639     /// Book a ND histogram with binning from a reference scatter.
0640     template<size_t DbnN, typename... AxisT>
0641     BinnedDbnPtr<DbnN,  AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const std::string& name,
0642                                         const YODA::BinnedEstimate<AxisT...>& refest) {
0643       const string path = histoPath(name);
0644 
0645       YODA::BinnedDbn<DbnN, AxisT...> yao(refest.binning(), path);
0646       for (const string& a : yao.annotations()) {
0647         if (a != "Path")  yao.rmAnnotation(a);
0648       }
0649       _setWriterPrecision(path, yao);
0650       return ao = registerAO(yao);
0651     }
0652 
0653     /// Book a ND histogram, using the binnings in the reference data histogram.
0654     template<size_t DbnN, typename... AxisT>
0655     BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const std::string& name) {
0656       return book(ao, name, refData<YODA::BinnedEstimate<AxisT...>>(name));
0657     }
0658 
0659     /// Book a ND histogram, using the binnings in the reference data histogram.
0660     ///
0661     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
0662     template<size_t DbnN, typename... AxisT>
0663     BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const unsigned int datasetID,
0664                                        const unsigned int xAxisID, const unsigned int yAxisID) {
0665       const string name = mkAxisCode(datasetID, xAxisID, yAxisID);
0666       return book(ao, name);
0667     }
0668 
0669     /// @}
0670 
0671     /// @name HistoGroup booking
0672     /// @{
0673 
0674     template <typename GroupAxisT, typename... AxisT>
0675     HistoGroupPtr<GroupAxisT, AxisT...>& book(HistoGroupPtr<GroupAxisT, AxisT...>& ao,
0676                                               const std::vector<GroupAxisT>& edges,
0677                                               const std::vector<std::string>& names) {
0678       ao = make_shared<HistoGroup<GroupAxisT, AxisT...>>(edges);
0679       if (ao->numBins() != names.size()) {
0680         throw RangeError("Binning and reference-data names don't match!");
0681       }
0682       for (auto& b : ao->bins()) {
0683         const string& refname = names[b.index()-1];
0684         book(b, refname, refData<YODA::BinnedEstimate<AxisT...>>(refname));
0685       }
0686       return ao;
0687     }
0688 
0689     template <typename GroupAxisT, typename... AxisT>
0690     HistoGroupPtr<GroupAxisT, AxisT...>& book(HistoGroupPtr<GroupAxisT, AxisT...>& ao,
0691                                               const std::vector<GroupAxisT>& edges) {
0692       return ao = make_shared<HistoGroup<GroupAxisT, AxisT...>>(edges);
0693     }
0694 
0695     template <typename GroupAxisT, typename... AxisT>
0696     HistoGroupPtr<GroupAxisT, AxisT...>& book(HistoGroupPtr<GroupAxisT, AxisT...>& ao,
0697                                               std::initializer_list<GroupAxisT>&& edges) {
0698       return ao = make_shared<HistoGroup<GroupAxisT, AxisT...>>(std::move(edges));
0699     }
0700 
0701     /// @}
0702 
0703     /// @name BinnedEstimate booking
0704     /// @{
0705 
0706     /// Book a ND estimate with @a nbins uniformly distributed across the range @a lower - @a upper .
0707     template<typename... AxisT, typename = YODA::enable_if_all_CAxisT<AxisT...>>
0708     BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao,
0709                                       const std::string& name, const std::vector<size_t>& nbins,
0710                                       const std::vector<std::pair<double,double>>& loUpPairs) {
0711       if (nbins.size() != loUpPairs.size()) {
0712         throw RangeError("Vectors should have the same size!");
0713       }
0714       const string path = histoPath(name);
0715 
0716       YODA::BinnedEstimate<AxisT...> yao(nbins, loUpPairs, path);
0717       _setWriterPrecision(path, yao);
0718 
0719       return ao = registerAO(yao);
0720     }
0721 
0722 
0723     // Specialiation for 1D
0724     Estimate1DPtr& book(Estimate1DPtr& ao, const std::string& name,
0725                                            const size_t nbins, const double lower, const double upper) {
0726       return book(ao, name, vector<size_t>{nbins},
0727                   vector<pair<double,double>>{{lower,upper}});
0728     }
0729 
0730     // Specialiation for 2D
0731     Estimate2DPtr& book(Estimate2DPtr& ao, const std::string& name,
0732                                            const size_t nbinsX, const double lowerX, const double upperX,
0733                                            const size_t nbinsY, const double lowerY, const double upperY) {
0734       return book(ao, name, vector<size_t>{nbinsX,nbinsY},
0735                   vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}});
0736     }
0737 
0738     // Specialiation for 3D
0739     Estimate3DPtr& book(Estimate3DPtr& ao, const std::string& name,
0740                                            const size_t nbinsX, const double lowerX, const double upperX,
0741                                            const size_t nbinsY, const double lowerY, const double upperY,
0742                                            const size_t nbinsZ, const double lowerZ, const double upperZ) {
0743       return book(ao, name, vector<size_t>{nbinsX,nbinsY,nbinsZ},
0744                   vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}, {lowerZ,upperZ}});
0745     }
0746 
0747     /// Book a ND estimate with non-uniform bins defined by the vector of bin edges @a binedges .
0748     template<typename... AxisT>
0749     BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao, const std::string& name,
0750                                       const std::vector<AxisT>&... binedges) {
0751       const string path = histoPath(name);
0752       YODA::BinnedEstimate<AxisT...> yao(binedges..., path);
0753       _setWriterPrecision(path, yao);
0754 
0755       return ao = registerAO(yao);
0756     }
0757 
0758     /// Book a ND estimate with non-uniform bins defined by the vector of bin edges @a binedges .
0759     template<typename... AxisT>
0760     BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao, const std::string& name,
0761                                       const std::initializer_list<AxisT>&... binedges) {
0762       return book(ao, name, vector<AxisT>{binedges} ...);
0763     }
0764 
0765     /// Book a ND estimate, using the binnings in the reference data histogram.
0766     template<typename... AxisT>
0767     BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao, const std::string& name) {
0768 
0769       const string path = histoPath(name);
0770       YODA::BinnedEstimate<AxisT...> yao;
0771       try {
0772         yao = YODA::BinnedEstimate<AxisT...>(refData<YODA::BinnedEstimate<AxisT...>>(name).binning());
0773       } catch (...) {
0774         MSG_DEBUG("Couldn't retrieve reference binning, continue with nullary AO constructor.");
0775       }
0776       yao.setPath(path);
0777       _setWriterPrecision(path, yao);
0778       return ao = registerAO(yao);
0779     }
0780 
0781     /// Book a ND estimate, using the binnings in the reference data histogram.
0782     ///
0783     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
0784     template<typename... AxisT>
0785     BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao, const unsigned int datasetID,
0786                                       const unsigned int xAxisID, const unsigned int yAxisID) {
0787       const string name = mkAxisCode(datasetID, xAxisID, yAxisID);
0788       return book(ao, name);
0789     }
0790 
0791     /// @}
0792 
0793 
0794     /// @name Scatter booking
0795     /// @{
0796 
0797     /// @brief Book a N-dimensional data point set with the given name.
0798     ///
0799     /// @note Unlike histogram booking, scatter booking by default makes no
0800     /// attempt to use reference data to pre-fill the data object. If you want
0801     /// this, which is sometimes useful e.g. when the x-position is not really
0802     /// meaningful and can't be extracted from the data, then set the @a
0803     /// copy_pts parameter to true. This creates points to match the reference
0804     /// data's x values and errors, but with the y values and errors zeroed...
0805     /// assuming that there is a reference histo with the same name: if there
0806     /// isn't, an exception will be thrown.
0807     template<size_t N>
0808     ScatterNDPtr<N>& book(ScatterNDPtr<N>& snd, const string& name, const bool copy_pts = false) {
0809       const string path = histoPath(name);
0810       YODA::ScatterND<N> scat(path);
0811       if (copy_pts) {
0812         const YODA::ScatterND<N> ref = refData<YODA::EstimateND<N-1>>(name).mkScatter();
0813         for (YODA::PointND<N> p : ref.points()) {
0814           p.setVal(N-1, 0.0);
0815           p.setErr(N-1, 0.0);
0816           scat.addPoint(std::move(p));
0817         }
0818       }
0819       _setWriterPrecision(path, scat);
0820       return snd = registerAO(scat);
0821     }
0822 
0823     /// @brief Book a N-dimensional data point set, using the binnings in the reference data histogram.
0824     ///
0825     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
0826     ///
0827     /// @note Unlike histogram booking, scatter booking by default makes no
0828     /// attempt to use reference data to pre-fill the data object. If you want
0829     /// this, which is sometimes useful e.g. when the x-position is not really
0830     /// meaningful and can't be extracted from the data, then set the @a
0831     /// copy_pts parameter to true. This creates points to match the reference
0832     /// data's x values and errors, but with the y values and errors zeroed.
0833     template<size_t N>
0834     ScatterNDPtr<N>& book(ScatterNDPtr<N>& snd, const unsigned int datasetID, const unsigned int xAxisID,
0835                                                 const unsigned int yAxisID, const bool copy_pts = false) {
0836       const string axisCode = mkAxisCode(datasetID, xAxisID, yAxisID);
0837       return book(snd, axisCode, copy_pts);
0838     }
0839 
0840     /// @brief Book a N-dimensional data point set with equally spaced x-points in a range.
0841     ///
0842     /// The y values and errors will be set to 0.
0843     ///
0844     /// @todo Remove this when we switch to BinnedEstimates
0845     Scatter2DPtr& book(Scatter2DPtr& snd, const string& name,
0846                        const size_t npts, const double lower, const double upper) {
0847       const string path = histoPath(name);
0848 
0849       Scatter2D scat(path);
0850       const double binwidth = (upper-lower)/npts;
0851       for (size_t pt = 0; pt < npts; ++pt) {
0852         const double bincentre = lower + (pt + 0.5) * binwidth;
0853         scat.addPoint(bincentre, 0, binwidth/2.0, 0);
0854       }
0855       _setWriterPrecision(path, scat);
0856 
0857       return snd = registerAO(scat);
0858     }
0859     //
0860     Scatter3DPtr& book(Scatter3DPtr& snd, const string& name,
0861                        const size_t nptsX, const double lowerX, const double upperX,
0862                        const size_t nptsY, const double lowerY, const double upperY) {
0863       const string path = histoPath(name);
0864 
0865       Scatter3D scat(path);
0866       const double xbinwidth = (upperX-lowerX)/nptsX;
0867       const double ybinwidth = (upperY-lowerY)/nptsY;
0868       for (size_t xpt = 0; xpt < nptsX; ++xpt) {
0869         const double xbincentre = lowerX + (xpt + 0.5) * xbinwidth;
0870         for (size_t ypt = 0; ypt < nptsY; ++ypt) {
0871           const double ybincentre = lowerY + (ypt + 0.5) * ybinwidth;
0872           scat.addPoint(xbincentre, ybincentre, 0, 0.5*xbinwidth, 0.5*ybinwidth, 0);
0873         }
0874       }
0875       _setWriterPrecision(path, scat);
0876 
0877       return snd = registerAO(scat);
0878     }
0879 
0880     /// @brief Book a 2-dimensional data point set based on provided contiguous "bin edges".
0881     ///
0882     /// The y values and errors will be set to 0.
0883     ///
0884     /// @todo Remove this when we switch to BinnedEstimates
0885     Scatter2DPtr& book(Scatter2DPtr& snd, const string& name,
0886                        const std::vector<double>& binedges) {
0887       const string path = histoPath(name);
0888 
0889       Scatter2D scat(path);
0890       for (size_t pt = 0; pt < binedges.size()-1; ++pt) {
0891         const double bincentre = (binedges[pt] + binedges[pt+1]) / 2.0;
0892         const double binwidth = binedges[pt+1] - binedges[pt];
0893         scat.addPoint(bincentre, 0, binwidth/2.0, 0);
0894       }
0895       _setWriterPrecision(path, scat);
0896 
0897       return snd = registerAO(scat);
0898     }
0899     //
0900     Scatter3DPtr& book(Scatter3DPtr& snd, const string& name,
0901                        const std::vector<double>& binedgesX,
0902                        const std::vector<double>& binedgesY) {
0903       const string path = histoPath(name);
0904 
0905       Scatter3D scat(path);
0906       for (size_t xpt = 0; xpt < binedgesX.size()-1; ++xpt) {
0907         const double xbincentre = (binedgesX[xpt] + binedgesX[xpt+1]) / 2.0;
0908         const double xbinwidth = binedgesX[xpt+1] - binedgesX[xpt];
0909         for (size_t ypt = 0; ypt < binedgesY.size()-1; ++ypt) {
0910           const double ybincentre = (binedgesY[ypt] + binedgesY[ypt+1]) / 2.0;
0911           const double ybinwidth = binedgesY[ypt+1] - binedgesY[ypt];
0912           scat.addPoint(xbincentre, ybincentre, 0, 0.5*xbinwidth, 0.5*ybinwidth, 0);
0913         }
0914       }
0915       _setWriterPrecision(path, scat);
0916 
0917       return snd = registerAO(scat);
0918     }
0919 
0920     /// Book a 2-dimensional data point set with x-points from an existing scatter and a new path.
0921     template<size_t N>
0922     ScatterNDPtr<N>& book(ScatterNDPtr<N>& snd, const string& name, const YODA::ScatterND<N>& refscatter) {
0923       const string path = histoPath(name);
0924 
0925       YODA::ScatterND<N> scat(refscatter, path);
0926       for (const string& a : scat.annotations()) {
0927         if (a != "Path")  scat.rmAnnotation(a);
0928       }
0929       _setWriterPrecision(path, scat);
0930 
0931       return snd = registerAO(scat);
0932     }
0933 
0934     /// @}
0935 
0936     /// @defgroup Cutflow booking
0937     /// @{
0938 
0939     /// Book a Cutflow object defined by the vector of @a edges
0940     CutflowPtr& book(CutflowPtr& ao, const string& name, const std::vector<std::string>& edges) {
0941       const string path = histoPath(name);
0942       Cutflow yao(edges, path);
0943       _setWriterPrecision(path, yao);
0944       return ao = registerAO(yao);
0945     }
0946 
0947     /// Book a Cutflow object defined by the vector of @a edges
0948     CutflowPtr& book(CutflowPtr& ao, const string& name, const std::initializer_list<std::string>& edges) {
0949       return book(ao, name, vector<std::string>{edges});
0950     }
0951 
0952     /// @}
0953 
0954     /// @name Cutflows booking
0955     /// @{
0956 
0957     CutflowsPtr& book(CutflowsPtr& ao, const std::vector<std::string>& edges,
0958                                            const std::vector<std::vector<std::string>>& innerEdges) {
0959       ao = make_shared<Cutflows>(edges);
0960       if (ao->numBins() !=innerEdges.size()) {
0961         throw RangeError("Outer and Inner edges don't match");
0962       }
0963       for (auto& b : ao->bins()) {
0964         book(b, b.xEdge(), innerEdges[b.index()-1]);
0965       }
0966       return ao;
0967     }
0968 
0969     CutflowsPtr& book(CutflowsPtr& ao, const std::vector<std::string>& edges) {
0970       return ao = make_shared<Cutflows>(edges);
0971     }
0972 
0973     CutflowsPtr& book(CutflowsPtr& ao, std::initializer_list<std::string>&& edges) {
0974       return ao = make_shared<Cutflows>(std::move(edges));
0975     }
0976 
0977     /// @}
0978 
0979 
0980     /// @name Virtual helper function to allow classes deriving
0981     /// from Analysis (e.g. CumulantAnalysis) to load external
0982     /// raw AOs into their local AOs (needed in heavy-ion land).
0983     virtual void rawHookIn(YODA::AnalysisObjectPtr yao) {
0984       (void) yao; // suppress unused variable warning
0985     }
0986 
0987     /// @name Virtual helper function to allow classes deriving from
0988     /// Analysis (e.g. CumulantAnalysis) to fiddle with raw AOs
0989     /// post-finalize/before writing them out (needed in heavy-ion land).
0990     virtual void rawHookOut(const vector<MultiplexAOPtr>& raos, size_t iW) {
0991       (void) raos; // suppress unused variable warning
0992       (void) iW; // suppress unused variable warning
0993     }
0994 
0995 
0996   public:
0997 
0998     /// @name Accessing options for this Analysis instance.
0999     /// @{
1000 
1001     /// Return the map of all options given to this analysis.
1002     const std::map<std::string,std::string>& options() const {
1003       return _options;
1004     }
1005 
1006     /// Get an option for this analysis instance as a string.
1007     std::string getOption(std::string optname, string def="") const {
1008       if ( _options.find(optname) != _options.end() )
1009         return _options.find(optname)->second;
1010       return def;
1011     }
1012 
1013     /// @brief Sane overload for literal character strings (which don't play well with stringstream)
1014     ///
1015     /// Note this isn't a template specialisation, because we can't return a non-static
1016     /// char*, and T-as-return-type is built into the template function definition.
1017     std::string getOption(std::string optname, const char* def) {
1018       return getOption<std::string>(optname, def);
1019     }
1020 
1021     /// @brief Get an option for this analysis instance converted to a specific type
1022     ///
1023     /// The return type is given by the specified @a def value, or by an explicit template
1024     /// type-argument, e.g. getOption<double>("FOO", 3).
1025     ///
1026     /// @warning To avoid accidents, strings not convertible to the requested
1027     /// type will throw a Rivet::ReadError exception.
1028     template<typename T>
1029     T getOption(std::string optname, T def) const {
1030       if (_options.find(optname) == _options.end()) return def;
1031       std::stringstream ss;
1032       ss.exceptions(std::ios::failbit);
1033       T ret;
1034       ss << _options.find(optname)->second;
1035       try {
1036         ss >> ret;
1037       } catch (...) {
1038         throw ReadError("Could not read user-provided option into requested type");
1039       }
1040       return ret;
1041     }
1042 
1043     /// @brief Get an option for this analysis instance converted to a bool
1044     ///
1045     /// Specialisation for bool, to allow use of "yes/no", "true/false"
1046     /// and "on/off" strings, with fallback casting to bool based on int
1047     /// value. An empty value will be treated as false.
1048     ///
1049     /// @warning To avoid accidents, strings not matching one of the above
1050     /// patterns will throw a Rivet::ReadError exception.
1051     ///
1052     /// @todo Make this a template-specialisation... needs to be outside the class body?
1053     // template<>
1054     // bool getOption<bool>(std::string optname, bool def) const {
1055     bool getOption(std::string optname, bool def) const {
1056       if (_options.find(optname) == _options.end()) return def;
1057       const std::string val = getOption(optname);
1058       const std::string lval = toLower(val);
1059       if (lval.empty()) return false;
1060       if (lval == "true" || lval == "yes" || lval == "on") return true;
1061       if (lval == "false" || lval == "no" || lval == "off") return false;
1062       return bool(getOption<int>(optname, 0));
1063     }
1064 
1065     /// @}
1066 
1067 
1068     /// @name Booking heavy ion features
1069     /// @{
1070 
1071     /// @brief Book a CentralityProjection
1072     ///
1073     /// Using a SingleValueProjection, @a proj, giving the value of an
1074     /// experimental observable to be used as a centrality estimator,
1075     /// book a CentralityProjection based on the experimentally
1076     /// measured pecentiles of this observable (as given by the
1077     /// reference data for the @a calHistName histogram in the @a
1078     /// calAnaName analysis. If a preloaded file with the output of a
1079     /// run using the @a calAnaName analysis contains a valid
1080     /// generated @a calHistName histogram, it will be used as an
1081     /// optional percentile binning. Also if this preloaded file
1082     /// contains a histogram with the name @a calHistName with an
1083     /// appended "_IMP" This histogram will be used to add an optional
1084     /// centrality percentile based on the generated impact
1085     /// parameter. If @a increasing is true, a low (high) value of @a proj
1086     /// is assumed to correspond to a more peripheral (central) event.
1087     const CentralityProjection&
1088     declareCentrality(const SingleValueProjection &proj,
1089                       string calAnaName, string calHistName,
1090                       const string projName,
1091               PercentileOrder pctorder=PercentileOrder::DECREASING);
1092 
1093 
1094     /// @brief Book a Percentile Multiplexer around AnalysisObjects.
1095     ///
1096     /// Based on a previously registered CentralityProjection named @a
1097     /// projName book one AnalysisObject for each @a centralityBin and
1098     /// name them according to the corresponding code in the @a ref
1099     /// vector.
1100     template<typename T>
1101     Percentile<T> book(const string& projName,
1102                        const vector<pair<double, double>>& centralityBins,
1103                        const vector<tuple<size_t, size_t, size_t>>& ref) {
1104 
1105       using RefT = typename ReferenceTraits<T>::RefT;
1106       using WrapT = MultiplexPtr<Multiplexer<T>>;
1107 
1108       Percentile<T> pctl(this, projName);
1109 
1110       const size_t nCent = centralityBins.size();
1111       for (size_t iCent = 0; iCent < nCent; ++iCent) {
1112         const string axisCode = mkAxisCode(std::get<0>(ref[iCent]),
1113                                            std::get<1>(ref[iCent]),
1114                                            std::get<2>(ref[iCent]));
1115         const RefT& refscatter = refData<RefT>(axisCode);
1116 
1117         WrapT wtf(_weightNames(), T(refscatter, histoPath(axisCode)));
1118         wtf = addAnalysisObject(wtf);
1119 
1120         CounterPtr cnt(_weightNames(), Counter(histoPath("TMP/COUNTER/" + axisCode)));
1121         cnt = addAnalysisObject(cnt);
1122 
1123         pctl.add(wtf, cnt, centralityBins[iCent]);
1124       }
1125       return pctl;
1126     }
1127 
1128 
1129     // /// @brief Book Percentile Multiplexers around AnalysisObjects.
1130     // ///
1131     // /// Based on a previously registered CentralityProjection named @a
1132     // /// projName book one (or several) AnalysisObject(s) named
1133     // /// according to @a ref where the x-axis will be filled according
1134     // /// to the percentile output(s) of the @projName.
1135     // ///
1136     // /// @todo Convert to just be called book() cf. others
1137     // template <class T>
1138     // PercentileXaxis<T> bookPercentileXaxis(string projName,
1139     //                                        tuple<int, int, int> ref) {
1140 
1141     //   typedef typename ReferenceTraits<T>::RefT RefT;
1142     //   typedef MultiplexPtr<Multiplexer<T>> WrapT;
1143 
1144     //   PercentileXaxis<T> pctl(this, projName);
1145 
1146     //   const string axisCode = mkAxisCode(std::get<0>(ref),
1147     //                                      std::get<1>(ref),
1148     //                                      std::get<2>(ref));
1149     //   const RefT & refscatter = refData<RefT>(axisCode);
1150 
1151     //   WrapT wtf(_weightNames(), T(refscatter, histoPath(axisCode)));
1152     //   wtf = addAnalysisObject(wtf);
1153 
1154     //   CounterPtr cnt(_weightNames(), Counter());
1155     //   cnt = addAnalysisObject(cnt);
1156 
1157     //   pctl.add(wtf, cnt);
1158     //   return pctl;
1159     // }
1160 
1161     /// @}
1162 
1163 
1164   private:
1165 
1166     // Functions that have to be defined in the .cc file to avoid circular #includes
1167 
1168     /// Get the list of weight names from the handler
1169     vector<string> _weightNames() const;
1170 
1171     /// Get the list of weight names from the handler
1172     YODA::AnalysisObjectPtr _getPreload (const string& name) const;
1173 
1174     /// Get an AO from another analysis
1175     MultiplexAOPtr _getOtherAnalysisObject(const std::string & ananame, const std::string& name);
1176 
1177     /// Check that analysis objects aren't being booked/registered outside the init stage
1178     void _checkBookInit() const;
1179 
1180     /// Check if we are in the init stage.
1181     bool _inInit() const;
1182 
1183     /// Check if we are in the finalize stage.
1184     bool _inFinalize() const;
1185 
1186     /// Set DP annotation
1187     template <typename YODAT>
1188     void _setWriterPrecision(const string& path, YODAT& yao) {
1189       const string re = _info->writerDoublePrecision();
1190       if (re != "") {
1191         std::smatch match;
1192         const bool needsDP = std::regex_search(path, match, std::regex(re));
1193         if (needsDP)  yao.template setAnnotation("WriterDoublePrecision", "1");
1194       }
1195     }
1196 
1197 
1198   private:
1199 
1200     /// To be used in finalize context only
1201     class CounterAdapter {
1202     public:
1203 
1204       CounterAdapter(double x) : x_(x) {}
1205 
1206       CounterAdapter(const YODA::Counter& c) : x_(c.val()) {}
1207 
1208       CounterAdapter(const YODA::Estimate& e) : x_(e.val()) {}
1209 
1210       CounterAdapter(const YODA::Scatter1D& s) : x_(s.points()[0].x()) {
1211         if (s.numPoints() != 1)  throw RangeError("Can only scale by a single value.");
1212       }
1213 
1214       operator double() const { return x_; }
1215 
1216     private:
1217       double x_;
1218 
1219     };
1220 
1221 
1222   public:
1223 
1224     double dbl(double x) { return x; }
1225     double dbl(const YODA::Counter& c) { return c.val(); }
1226     double dbl(const YODA::Estimate0D& e) { return e.val(); }
1227     double dbl(const YODA::Scatter1D& s) {
1228       if ( s.numPoints() != 1 ) throw RangeError("Only scatter with single value supported.");
1229       return s.points()[0].x();
1230     }
1231 
1232 
1233   protected:
1234 
1235     /// @name Analysis object manipulation
1236     ///
1237     /// @{
1238 
1239     /// Multiplicatively scale the given AnalysisObject, @a ao, by factor @a factor.
1240     template<typename T>
1241     void scale(MultiplexPtr<Multiplexer<T>>& ao, CounterAdapter factor) {
1242       if (!ao) {
1243         MSG_WARNING("Failed to scale AnalysisObject=NULL in analysis "
1244                     << name() << " (scale=" << double(factor) << ")");
1245         return;
1246       }
1247       if (std::isnan(double(factor)) || std::isinf(double(factor))) {
1248         MSG_WARNING("Failed to scale AnalysisObject=" << ao->path() << " in analysis: "
1249                     << name() << " (invalid scale factor = " << double(factor) << ")");
1250         factor = 0;
1251       }
1252       MSG_TRACE("Scaling AnalysisObject " << ao->path() << " by factor " << double(factor));
1253       try {
1254         if constexpr( isFillable<T>::value ) {
1255           ao->scaleW(factor);
1256         }
1257         else {
1258           ao->scale(factor);
1259         }
1260       }
1261       catch (YODA::Exception& we) {
1262         MSG_WARNING("Could not scale AnalysisObject " << ao->path());
1263         return;
1264       }
1265     }
1266 
1267     /// Multiplicatively scale the given histogram group, @a group, by factor @a factor.
1268     template<typename GroupAxisT, typename... AxisT>
1269     void scale(HistoGroupPtr<GroupAxisT, AxisT...>& group, CounterAdapter factor) {
1270       if (!group) {
1271         MSG_WARNING("Failed to scale AnalysisObject=NULL in analysis "
1272                     << name() << " (scale=" << double(factor) << ")");
1273         return;
1274       }
1275       if (std::isnan(double(factor)) || std::isinf(double(factor))) {
1276         MSG_WARNING("Failed to scale histo group in analysis: "
1277                     << name() << " (invalid scale factor = " << double(factor) << ")");
1278         factor = 0;
1279       }
1280       MSG_TRACE("Scaling histo group by factor " << double(factor));
1281       try {
1282         group->scaleW(factor);
1283       }
1284       catch (YODA::Exception& we) {
1285         MSG_WARNING("Could not scale histo group.");
1286         return;
1287       }
1288     }
1289 
1290     /// Multiplicatively scale the given histogram group, @a group, by factors @a factors.
1291     template<typename GroupAxisT, typename... AxisT>
1292     void scale(HistoGroupPtr<GroupAxisT, AxisT...>& group, const vector<double>& factors) {
1293       if (!group) {
1294         MSG_WARNING("Failed to scale AnalysisObject=NULL in analysis " << name());
1295         return;
1296       }
1297       if (group->numBins(true) != factors.size()) {
1298         throw RangeError(name() + ": Number of scale factors does not match group binning");
1299         return;
1300       }
1301       for (auto& b : group->bins(true)) {
1302         if (!b.get())  continue;
1303         double factor = factors[b.index()];
1304         if (std::isnan(factor) || std::isinf(factor)) {
1305           MSG_WARNING("Failed to scale componment of histo group in analysis: "
1306                       << name() << " (invalid scale factor = " << factor << ")");
1307           factor = 0;
1308         }
1309         MSG_TRACE("Scaling histo group element by factor " << factor);
1310         try {
1311           b->scaleW(factor);
1312         }
1313         catch (YODA::Exception& we) {
1314           MSG_WARNING("Could not scale component of histo group.");
1315         }
1316       }
1317     }
1318 
1319     /// Multiplicatively scale the cutflow group, @a group, by factor @a factor.
1320     void scale(CutflowsPtr& group, CounterAdapter factor) {
1321       if (!group) {
1322         MSG_WARNING("Failed to scale AnalysisObject=NULL in analysis "
1323                     << name() << " (scale=" << double(factor) << ")");
1324         return;
1325       }
1326       if (std::isnan(double(factor)) || std::isinf(double(factor))) {
1327         MSG_WARNING("Failed to scale histo group in analysis: "
1328                     << name() << " (invalid scale factor = " << double(factor) << ")");
1329         factor = 0;
1330       }
1331       MSG_TRACE("Scaling histo group by factor " << double(factor));
1332       try {
1333         group->scale(factor);
1334       }
1335       catch (YODA::Exception& we) {
1336         MSG_WARNING("Could not scale histo group.");
1337         return;
1338       }
1339     }
1340 
1341     /// Iteratively scale the AOs in the map @a aos, by factor @a factor.
1342     template<typename T, typename U>
1343     void scale(std::map<T, U>& aos, CounterAdapter factor) {
1344       for (auto& item : aos)  scale(item.second, factor);
1345     }
1346 
1347     /// Iteratively scale the AOs in the iterable @a aos, by factor @a factor.
1348     template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
1349     void scale(AORange& aos, CounterAdapter factor) {
1350       for (auto& ao : aos)  scale(ao, factor);
1351     }
1352 
1353     /// Iteratively scale the AOs in the initialiser list @a aos, by factor @a factor.
1354     template <typename T>
1355     void scale(std::initializer_list<T> aos, CounterAdapter factor) {
1356       for (auto& ao : std::vector<T>{aos})  scale(ao, factor);
1357     }
1358 
1359     /// Iteratively scale the AOs in the map @a aos, by factors @a factors.
1360     template<typename T, typename U>
1361     void scale(std::map<T, U>& aos, const vector<double>& factors) {
1362       for (auto& item : aos)  scale(item.second, factors);
1363     }
1364 
1365     /// Iteratively scale the AOs in the iterable @a aos, by factors @a factors.
1366     template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
1367     void scale(AORange& aos, const vector<double>& factors) {
1368       for (auto& ao : aos)  scale(ao, factors);
1369     }
1370 
1371     /// Iteratively scale the AOs in the initialiser list @a aos, by factors @a factors.
1372     template <typename T>
1373     void scale(std::initializer_list<T> aos, const vector<double>& factors) {
1374       for (auto& ao : std::vector<T>{aos})  scale(ao, factors);
1375     }
1376 
1377     /// Scale the given histogram group, @a group, by the group axis width
1378     template<typename GroupAxisT, typename... AxisT>
1379     void divByGroupWidth(HistoGroupPtr<GroupAxisT, AxisT...>& group) {
1380       if (!group) {
1381         MSG_WARNING("Failed to scale HistoGroup=NULL in analysis "
1382                     << name() << " by group axis width");
1383         return;
1384       }
1385       group->divByGroupWidth();
1386     }
1387 
1388     /// Iteratively scale the HistoGroups in the map @a aos, by the group axis width
1389     template<typename T, typename U>
1390     void divByGroupWidth(std::map<T, U>& aos) {
1391       for (auto& item : aos)  divByGroupWidth(item.second);
1392     }
1393 
1394     /// Iteratively scale the HistoGroups in the iterable @a aos, by the group axis width
1395     template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
1396     void divByGroupWidth(AORange& aos) {
1397       for (auto& ao : aos)  divByGroupWidth(ao);
1398     }
1399 
1400     /// Iteratively scale the HistoGroups in the initialiser list @a aos, by the group axis width
1401     template <typename T>
1402     void divByGroupWidth(std::initializer_list<T> aos) {
1403       for (auto& ao : std::vector<T>{aos})  divByGroupWidth(ao);
1404     }
1405 
1406 
1407     /// Normalize the given analysis object, @a ao to a target @a norm.
1408     template <size_t DbnN, typename... AxisT>
1409     void normalize(BinnedDbnPtr<DbnN, AxisT...> ao, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1410       if (!ao) {
1411         MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")");
1412         return;
1413       }
1414       MSG_TRACE("Normalizing histo " << ao->path() << " to " << double(norm));
1415       try {
1416         const double hint = ao->integral(includeoverflows);
1417         if (hint == 0)  MSG_DEBUG("Skipping histo with null area " << ao->path());
1418         else            ao->normalize(norm, includeoverflows);
1419       }
1420       catch (YODA::Exception& we) {
1421         MSG_WARNING("Could not normalize histo " << ao->path());
1422         return;
1423       }
1424     }
1425 
1426     /// Normalize each AO in the given histogram group, @a group to a target @a norm.
1427     template <typename GroupAxisT, typename... AxisT>
1428     void normalize(HistoGroupPtr<GroupAxisT, AxisT...> group, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1429       if (!group) {
1430         MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")");
1431         return;
1432       }
1433       MSG_TRACE("Normalizing histo group  to " << double(norm));
1434       try {
1435         const double hint = group->integral(includeoverflows);
1436         if (hint == 0)  MSG_DEBUG("Skipping histo group with null area.");
1437         else            group->normalize(norm, includeoverflows);
1438       }
1439       catch (YODA::Exception& we) {
1440         MSG_WARNING("Could not normalize histo group.");
1441         return;
1442       }
1443     }
1444 
1445 
1446     /// Iteratively normalise the AOs in the iterable @a iter, by factor @a factor.
1447     template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
1448     void normalize(AORange& aos, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1449       for (auto& ao : aos)  normalize(ao, norm, includeoverflows);
1450     }
1451 
1452     /// Iteratively normalise the AOs in the initialiser list @a iter to a target @a norm.
1453     template<typename T>
1454     void normalize(std::initializer_list<T>&& aos, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1455       for (auto& ao : aos)  normalize(ao, norm, includeoverflows);
1456     }
1457 
1458     /// Iteratively normalise the AOs in the map @a aos to a target @a norm.
1459     template<typename T, typename U>
1460     void normalize(std::map<T, U>& aos, //BinnedDbnPtr<DbnN, AxisT...>>& aos,
1461                    const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1462       for (auto& item : aos)  normalize(item.second, norm, includeoverflows);
1463     }
1464 
1465     /// Normalize the given histogram group, @a group to a target @a norm.
1466     template <typename GroupAxisT, typename... AxisT>
1467     void normalizeGroup(HistoGroupPtr<GroupAxisT, AxisT...> group, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1468       if (!group) {
1469         MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")");
1470         return;
1471       }
1472       MSG_TRACE("Normalizing histo group  to " << double(norm));
1473       try {
1474         const double hint = group->integral(includeoverflows);
1475         if (hint == 0)  MSG_DEBUG("Skipping histo group with null area.");
1476         else            group->normalizeGroup(norm, includeoverflows);
1477       }
1478       catch (YODA::Exception& we) {
1479         MSG_WARNING("Could not normalize histo group.");
1480         return;
1481       }
1482     }
1483 
1484     /// Iteratively normalise the HistoGroups in the iterable @a iter, by factor @a factor.
1485     template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
1486     void normalizeGroup(AORange& aos, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1487       for (auto& ao : aos)  normalizeGroup(ao, norm, includeoverflows);
1488     }
1489 
1490     /// Iteratively normalise the HistoGroups in the initialiser list @a iter to a target @a norm.
1491     template<typename T>
1492     void normalizeGroup(std::initializer_list<T>&& aos, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1493       for (auto& ao : aos)  normalizeGroup(ao, norm, includeoverflows);
1494     }
1495 
1496     /// Iteratively normalise the HistoGroups in the map @a aos to a target @a norm.
1497     template<typename T, typename U>
1498     void normalizeGroup(std::map<T, U>& aos, //BinnedDbnPtr<DbnN, AxisT...>>& aos,
1499                    const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1500       for (auto& item : aos)  normalizeGroup(item.second, norm, includeoverflows);
1501     }
1502 
1503 
1504     /// Helper for histogram conversion to an inert estimate type
1505     ///
1506     /// @note Assigns to the (already registered) output estimate, @a est.
1507     /// Preserves the path information of the target.
1508     template<size_t DbnN, typename... AxisT>
1509     void barchart(BinnedDbnPtr<DbnN, AxisT...> ao, BinnedEstimatePtr<AxisT...> est) const {
1510       const string path = est->path();
1511       *est = ao->mkEstimate(path, "stats", false); //< do NOT divide by bin area cf. a differential dsigma/dX histogram
1512     }
1513 
1514     /// Helper for counter division.
1515     ///
1516     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1517     void divide(CounterPtr c1, CounterPtr c2, Estimate0DPtr est) const;
1518 
1519     /// Helper for histogram division with raw YODA objects.
1520     ///
1521     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1522     void divide(const YODA::Counter& c1, const YODA::Counter& c2, Estimate0DPtr est) const;
1523 
1524     /// Helper for counter division.
1525     ///
1526     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1527     void divide(Estimate0DPtr e1, Estimate0DPtr e2, Estimate0DPtr est) const;
1528 
1529     /// Helper for estimate division with raw YODA objects.
1530     ///
1531     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1532     void divide(const YODA::Estimate0D& e1, const YODA::Estimate0D& e2, Estimate0DPtr est) const;
1533 
1534 
1535     /// Helper for histogram division.
1536     ///
1537     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1538     template<size_t DbnN, typename... AxisT>
1539     void divide(const YODA::BinnedDbn<DbnN, AxisT...>& h1, const YODA::BinnedDbn<DbnN, AxisT...>& h2,
1540                 BinnedEstimatePtr<AxisT...> est) const {
1541       const string path = est->path();
1542       *est = h1 / h2;
1543       est->setPath(path);
1544     }
1545     //
1546     template<size_t DbnN, typename... AxisT>
1547     void divide(BinnedDbnPtr<DbnN, AxisT...> h1, BinnedDbnPtr<DbnN, AxisT...> h2,
1548                              BinnedEstimatePtr<AxisT...> est) const {
1549       return divide(*h1, *h2, est);
1550     }
1551 
1552     /// Helper for binned estimate division.
1553     ///
1554     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1555     template<typename... AxisT>
1556     void divide(const YODA::BinnedEstimate<AxisT...>& e1, const YODA::BinnedEstimate<AxisT...>& e2,
1557                 BinnedEstimatePtr<AxisT...> est) const {
1558       const string path = est->path();
1559       *est = e1 / e2;
1560       est->setPath(path);
1561     }
1562     //
1563     template<typename... AxisT>
1564     void divide(BinnedEstimatePtr<AxisT...> e1, BinnedEstimatePtr<AxisT...> e2,
1565                              BinnedEstimatePtr<AxisT...> est) const {
1566       return divide(*e1, *e2, est);
1567     }
1568 
1569 
1570 
1571     /// Helper for counter efficiency calculation.
1572     ///
1573     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1574     void efficiency(CounterPtr c1, CounterPtr c2, Estimate0DPtr est) const {
1575       efficiency(*c1, *c2, est);
1576     }
1577 
1578     /// Helper for counter efficiency calculation.
1579     ///
1580     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1581     void efficiency(const YODA::Counter& c1, const YODA::Counter& c2, Estimate0DPtr est) const {
1582       const string path = est->path();
1583       *est = YODA::efficiency(c1, c2);
1584       est->setPath(path);
1585     }
1586 
1587 
1588 
1589     /// Helper for histogram efficiency calculation.
1590     ///
1591     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1592     template<size_t DbnN, typename... AxisT>
1593     void efficiency(const YODA::BinnedDbn<DbnN, AxisT...>& h1, const YODA::BinnedDbn<DbnN, AxisT...>& h2,
1594                     BinnedEstimatePtr<AxisT...> est) const {
1595       const string path = est->path();
1596       *est = YODA::efficiency(h1, h2);
1597       est->setPath(path);
1598     }
1599     //
1600     template<size_t DbnN, typename... AxisT>
1601     void efficiency(BinnedDbnPtr<DbnN, AxisT...> h1, BinnedDbnPtr<DbnN, AxisT...> h2,
1602                     BinnedEstimatePtr<AxisT...> est) const {
1603       efficiency(*h1, *h2, est);
1604     }
1605 
1606 
1607     /// Helper for estimate efficiency calculation.
1608     ///
1609     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1610     template<typename... AxisT>
1611     void efficiency(const YODA::BinnedEstimate<AxisT...>& e1, const YODA::BinnedEstimate<AxisT...>& e2,
1612                     BinnedEstimatePtr<AxisT...> est) const {
1613       const string path = est->path();
1614       *est = YODA::efficiency(e1, e2);
1615       est->setPath(path);
1616     }
1617     //
1618     template<typename... AxisT>
1619     void efficiency(BinnedEstimatePtr<AxisT...> e1, BinnedEstimatePtr<AxisT...> e2,
1620                     BinnedEstimatePtr<AxisT...> est) const {
1621       efficiency(*e1, *e2, est);
1622     }
1623 
1624 
1625     /// Helper for histogram asymmetry calculation.
1626     ///
1627     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1628     template<size_t DbnN, typename... AxisT>
1629     void asymm(const YODA::BinnedDbn<DbnN, AxisT...>& h1, const YODA::BinnedDbn<DbnN, AxisT...>& h2,
1630                BinnedEstimatePtr<AxisT...> est) const {
1631       const string path = est->path();
1632       *est = YODA::asymm(h1, h2);
1633       est->setPath(path);
1634     }
1635     //
1636     template<size_t DbnN, typename... AxisT>
1637     void asymm(BinnedDbnPtr<DbnN, AxisT...> h1, BinnedDbnPtr<DbnN, AxisT...> h2,
1638                BinnedEstimatePtr<AxisT...> est) const {
1639       asymm(*h1, *h2, est);
1640     }
1641 
1642     /// Helper for estimate asymmetry calculation.
1643     ///
1644     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1645     template<typename... AxisT>
1646     void asymm(const YODA::BinnedEstimate<AxisT...>& e1, const YODA::BinnedEstimate<AxisT...>& e2,
1647                BinnedEstimatePtr<AxisT...> est) const {
1648       const string path = est->path();
1649       *est = YODA::asymm(e1, e2);
1650       est->setPath(path);
1651     }
1652     //
1653     template<typename... AxisT>
1654     void asymm(BinnedEstimatePtr<AxisT...> e1, BinnedEstimatePtr<AxisT...> e2,
1655                BinnedEstimatePtr<AxisT...> est) const {
1656       asymm(*e1, *e2, est);
1657     }
1658 
1659     /// Helper for converting a differential histo to an integral one.
1660     ///
1661     /// @note Assigns to the (already registered) output estimate, @a est. Preserves the path information of the target.
1662     template<size_t DbnN, typename... AxisT>
1663     void integrate(const YODA::BinnedDbn<DbnN, AxisT...>& h, BinnedEstimatePtr<AxisT...> est) const {
1664       const string path = est->path();
1665       *est = mkIntegral(h);
1666       est->setPath(path);
1667     }
1668     //
1669     template<size_t DbnN, typename... AxisT>
1670     void integrate(BinnedDbnPtr<DbnN, AxisT...>& h, BinnedEstimatePtr<AxisT...> est) const {
1671       integrate(*h, est);
1672     }
1673 
1674     /// @}
1675 
1676 
1677   public:
1678 
1679     /// List of registered analysis data objects
1680     const vector<MultiplexAOPtr>& analysisObjects() const {
1681       return _analysisobjects;
1682     }
1683 
1684 
1685   protected:
1686 
1687     /// @name Data object registration, retrieval, and removal
1688     /// @{
1689 
1690     /// Get the default/nominal weight index
1691     size_t defaultWeightIndex() const;
1692 
1693     /// Get a preloaded YODA object.
1694     template <typename YODAT>
1695     shared_ptr<YODAT> getPreload(const string& path) const {
1696       return dynamic_pointer_cast<YODAT>(_getPreload(path));
1697     }
1698 
1699 
1700     /// Register a new data object, optionally read in preloaded data.
1701     template <typename YODAT>
1702     MultiplexPtr< Multiplexer<YODAT> > registerAO(const YODAT& yao) {
1703       using MultiplexerT = Multiplexer<YODAT>;
1704       using YODAPtrT = shared_ptr<YODAT>;
1705       using RAOT = MultiplexPtr<MultiplexerT>;
1706 
1707       if ( !_inInit() && !_inFinalize() ) {
1708         MSG_ERROR("Can't book objects outside of init() or finalize()");
1709         throw UserError(name() + ": Can't book objects outside of init() or finalize().");
1710       }
1711 
1712       // First check that we haven't booked this before.
1713       // This is allowed when booking in finalize: just warn in that case.
1714       // If in init(), throw an exception: it's 99.9% never going to be intentional.
1715       for (auto& waold : analysisObjects()) {
1716         if ( yao.path() == waold.get()->basePath() ) {
1717           const string msg = "Found double-booking of " + yao.path() + " in " + name();
1718           if ( _inInit() ) {
1719             MSG_ERROR(msg);
1720             throw LookupError(msg);
1721           } else {
1722             MSG_WARNING(msg + ". Keeping previous booking");
1723           }
1724           return RAOT(dynamic_pointer_cast<MultiplexerT>(waold.get()));
1725         }
1726       }
1727 
1728       shared_ptr<MultiplexerT> wao = make_shared<MultiplexerT>();
1729       wao->_basePath = yao.path();
1730       YODAPtrT yaop = make_shared<YODAT>(yao);
1731 
1732       for (const string& weightname : _weightNames()) {
1733         // Create two YODA objects for each weight. Copy from
1734         // preloaded YODAs if present. First the finalized yoda:
1735         string finalpath = yao.path();
1736         if ( weightname != "" ) finalpath +=  "[" + weightname + "]";
1737         YODAPtrT preload = getPreload<YODAT>(finalpath);
1738         if ( preload ) {
1739           if ( !bookingCompatible(preload, yaop) ) {
1740             /// @todo What about if/when we want to make the final objects the Scatter or binned persistent type?
1741             MSG_WARNING("Found incompatible pre-existing data object with same base path "
1742                         << finalpath <<  " for " << name());
1743             preload = nullptr;
1744           } else {
1745             MSG_TRACE("Using preloaded " << finalpath << " in " <<name());
1746             wao->_final.push_back(make_shared<YODAT>(*preload));
1747           }
1748         }
1749         else {
1750           wao->_final.push_back(make_shared<YODAT>(yao));
1751           wao->_final.back()->setPath(finalpath);
1752         }
1753 
1754         // Then the raw filling yodas.
1755         string rawpath = "/RAW" + finalpath;
1756         preload = getPreload<YODAT>(rawpath);
1757         if ( preload ) {
1758           if ( !bookingCompatible(preload, yaop) ) {
1759             MSG_WARNING("Found incompatible pre-existing data object with same base path "
1760                         << rawpath <<  " for " << name());
1761             preload = nullptr;
1762           } else {
1763             MSG_TRACE("Using preloaded " << rawpath << " in " <<name());
1764             wao->_persistent.push_back(make_shared<YODAT>(*preload));
1765           }
1766         }
1767         else {
1768           wao->_persistent.push_back(make_shared<YODAT>(yao));
1769           wao->_persistent.back()->setPath(rawpath);
1770         }
1771       }
1772       MultiplexPtr<MultiplexerT> ret(wao);
1773 
1774       ret.get()->unsetActiveWeight();
1775       if ( _inFinalize() ) {
1776         // If booked in finalize() we assume it is the first time
1777         // finalize is run.
1778         ret.get()->pushToFinal();
1779         ret.get()->setActiveFinalWeightIdx(0);
1780       }
1781       _analysisobjects.push_back(ret);
1782 
1783       return ret;
1784     }
1785 
1786 
1787     /// Register a data object in the histogram system
1788     template <typename AO=MultiplexAOPtr>
1789     AO addAnalysisObject(const AO& aonew) {
1790       _checkBookInit();
1791 
1792       for (const MultiplexAOPtr& ao : analysisObjects()) {
1793 
1794         // Check AO base-name first
1795         ao.get()->setActiveWeightIdx(defaultWeightIndex());
1796         aonew.get()->setActiveWeightIdx(defaultWeightIndex());
1797         if (ao->path() != aonew->path()) continue;
1798 
1799         // If base-name matches, check compatibility
1800         // NB. This evil is because dynamic_ptr_cast can't work on MultiplexPtr directly
1801         AO aoold = AO(dynamic_pointer_cast<typename AO::value_type>(ao.get())); //< OMG
1802         if ( !aoold || !bookingCompatible(aonew, aoold) ) {
1803           MSG_WARNING("Found incompatible pre-existing data object with same base path "
1804                       << aonew->path() <<  " for " << name());
1805           throw LookupError("Found incompatible pre-existing data object with same base path during AO booking");
1806         }
1807 
1808         // Finally, check all weight variations
1809         for (size_t weightIdx = 0; weightIdx < _weightNames().size(); ++weightIdx) {
1810           aoold.get()->setActiveWeightIdx(weightIdx);
1811           aonew.get()->setActiveWeightIdx(weightIdx);
1812           if (aoold->path() != aonew->path()) {
1813             MSG_WARNING("Found incompatible pre-existing data object with different weight-path "
1814                         << aonew->path() <<  " for " << name());
1815             throw LookupError("Found incompatible pre-existing data object with same weight-path during AO booking");
1816           }
1817         }
1818 
1819         // They're fully compatible: bind and return
1820         aoold.get()->unsetActiveWeight();
1821         MSG_TRACE("Bound pre-existing data object " << aoold->path() <<  " for " << name());
1822         return aoold;
1823       }
1824 
1825       // No equivalent found
1826       MSG_TRACE("Registered " << aonew->annotation("Type") << " " << aonew->path() <<  " for " << name());
1827       aonew.get()->unsetActiveWeight();
1828 
1829       _analysisobjects.push_back(aonew);
1830       return aonew;
1831     }
1832 
1833     /// Unregister a data object from the histogram system (by name)
1834     void removeAnalysisObject(const std::string& path);
1835 
1836     /// Unregister a data object from the histogram system (by pointer)
1837     void removeAnalysisObject(const MultiplexAOPtr& ao);
1838 
1839     /// Get a Rivet data object from the histogram system
1840     template <typename AO=MultiplexAOPtr>
1841     const AO getAnalysisObject(const std::string& aoname) const {
1842       for (const MultiplexAOPtr& ao : analysisObjects()) {
1843         ao.get()->setActiveWeightIdx(defaultWeightIndex());
1844         if (ao->path() == histoPath(aoname)) {
1845           // return dynamic_pointer_cast<AO>(ao);
1846           return AO(dynamic_pointer_cast<typename AO::value_type>(ao.get()));
1847         }
1848       }
1849       throw LookupError("Data object " + histoPath(aoname) + " not found");
1850     }
1851 
1852 
1853     // /// Get a data object from the histogram system
1854     // template <typename AO=YODA::AnalysisObject>
1855     // const std::shared_ptr<AO> getAnalysisObject(const std::string& name) const {
1856     //   foreach (const AnalysisObjectPtr& ao, analysisObjects()) {
1857     //     if (ao->path() == histoPath(name)) return dynamic_pointer_cast<AO>(ao);
1858     //   }
1859     //   throw LookupError("Data object " + histoPath(name) + " not found");
1860     // }
1861 
1862     // /// Get a data object from the histogram system (non-const)
1863     // template <typename AO=YODA::AnalysisObject>
1864     // std::shared_ptr<AO> getAnalysisObject(const std::string& name) {
1865     //   foreach (const AnalysisObjectPtr& ao, analysisObjects()) {
1866     //     if (ao->path() == histoPath(name)) return dynamic_pointer_cast<AO>(ao);
1867     //   }
1868     //   throw LookupError("Data object " + histoPath(name) + " not found");
1869     // }
1870 
1871 
1872     /// Get a data object from another analysis (e.g. preloaded
1873     /// calibration histogram).
1874     template <typename AO=MultiplexAOPtr>
1875     AO getAnalysisObject(const std::string& ananame,
1876                          const std::string& aoname) {
1877       MultiplexAOPtr ao = _getOtherAnalysisObject(ananame, aoname);
1878       // return dynamic_pointer_cast<AO>(ao);
1879       return AO(dynamic_pointer_cast<typename AO::value_type>(ao.get()));
1880     }
1881 
1882     /// @}
1883 
1884 
1885   private:
1886 
1887     /// Name passed to constructor (used to find .info analysis data file, and as a fallback)
1888     string _defaultname;
1889 
1890     /// Pointer to analysis metadata object
1891     unique_ptr<AnalysisInfo> _info;
1892 
1893     /// Storage of all plot objects
1894     /// @todo Make this a map for fast lookup by path?
1895     vector<MultiplexAOPtr> _analysisobjects;
1896 
1897     /// @name Cross-section variables
1898     /// @{
1899     double _crossSection;
1900     bool _gotCrossSection;
1901     /// @}
1902 
1903     /// The controlling AnalysisHandler object.
1904     AnalysisHandler* _analysishandler;
1905 
1906     /// Collection of cached refdata to speed up many autobookings: the
1907     /// reference data file should only be read once.
1908     mutable std::map<std::string, YODA::AnalysisObjectPtr> _refdata;
1909 
1910     /// Options the (this instance of) the analysis
1911     map<string, string> _options;
1912 
1913     /// The string of options.
1914     string _optstring;
1915 
1916 
1917   private:
1918 
1919     /// @name Utility functions
1920     /// @{
1921 
1922     /// Get the reference data for this paper and cache it.
1923     void _cacheRefData() const;
1924 
1925     /// @}
1926 
1927   };
1928 
1929 
1930   // // Template specialisation for literal character strings (which don't play well with stringstream)
1931   // template<>
1932   // inline std::string Analysis::getOption(std::string optname, const char* def) {
1933   //   return getOption<std::string>(optname, def); //.c_str();
1934   // }
1935 
1936 
1937 }
1938 
1939 
1940 // Include definition of analysis plugin system so that analyses automatically see it when including Analysis.hh
1941 #include "Rivet/AnalysisBuilder.hh"
1942 
1943 
1944 /// @defgroup anamacros Analysis macros
1945 /// @{
1946 
1947 /// @def RIVET_DECLARE_PLUGIN
1948 /// Preprocessor define to prettify the global-object plugin hook mechanism
1949 #define RIVET_DECLARE_PLUGIN(clsname) ::Rivet::AnalysisBuilder<clsname> plugin_ ## clsname
1950 
1951 /// @def RIVET_DECLARE_ALIASED_PLUGIN
1952 /// Preprocessor define to prettify the global-object plugin hook mechanism, with an extra alias name for this analysis
1953 #define RIVET_DECLARE_ALIASED_PLUGIN(clsname, alias) RIVET_DECLARE_PLUGIN(clsname)( #alias )
1954 
1955 /// @def RIVET_DEFAULT_ANALYSIS_CTOR
1956 /// Preprocessor define to prettify the awkward constructor with name string argument
1957 #define RIVET_DEFAULT_ANALYSIS_CTOR(clsname) clsname() : Analysis(# clsname) {}
1958 
1959 /// @def RIVET_REGISTER_TYPE
1960 /// Preprocessor define to prettify on-the-fly type registration
1961 #define RIVET_REGISTER_TYPE(...) handler().registerType<__VA_ARGS__>()
1962 
1963 /// @def RIVET_REGISTER_BINNED_SET
1964 /// Preprocessor define to prettify on-the-fly type registration
1965 #define RIVET_REGISTER_BINNED_SET(...) { \
1966     RIVET_REGISTER_TYPE(YODA::BinnedHisto<__VA_ARGS__>); \
1967     RIVET_REGISTER_TYPE(YODA::BinnedProfile<__VA_ARGS__>); \
1968     RIVET_REGISTER_TYPE(YODA::BinnedEstimate<__VA_ARGS__>); }
1969 
1970 /// @}
1971 
1972 
1973 #endif