Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-17 09:13:44

0001 // -*- C++ -*-
0002 #ifndef RIVET_RivetHandler_HH
0003 #define RIVET_RivetHandler_HH
0004 
0005 #include "Rivet/Config/RivetCommon.hh"
0006 #include "Rivet/Particle.hh"
0007 #include "Rivet/AnalysisLoader.hh"
0008 #include "Rivet/Tools/RivetYODA.hh"
0009 #include "Rivet/Tools/Utils.hh"
0010 #include "Rivet/ProjectionHandler.hh"
0011 #include "YODA/ReaderYODA.h"
0012 
0013 #include <fstream>
0014 #include <unordered_map>
0015 
0016 namespace Rivet {
0017 
0018 
0019   // Forward declaration and smart pointer for Analysis
0020   class Analysis;
0021   using AnaHandle = std::shared_ptr<Analysis>;
0022 
0023 
0024   /// @brief The key class for coordination of Analysis objects and the event loop
0025   ///
0026   /// A class which handles a number of analysis objects to be applied to
0027   /// generated events. An {@link Analysis}' AnalysisHandler is also responsible
0028   /// for handling the final writing-out of histograms.
0029   class AnalysisHandler {
0030 
0031     using TypeHandlePtr = std::shared_ptr<TypeBaseHandle>;
0032     using TypeRegister = std::unordered_map<string, TypeHandlePtr>;
0033     using TypeRegisterItr = typename TypeRegister::const_iterator;
0034 
0035   public:
0036 
0037     using Annotations = std::map<std::string, std::string>;
0038 
0039     /// Preferred / default constructor.
0040     AnalysisHandler();
0041 
0042     /// The copy constructor is deleted, so it can never be called.
0043     AnalysisHandler(const AnalysisHandler&) = delete;
0044 
0045     /// The assignment operator is deleted, so it can never be called.
0046     AnalysisHandler& operator=(const AnalysisHandler&) = delete;
0047 
0048     /// The destructor is not virtual, as this class should not be inherited from.
0049     ~AnalysisHandler();
0050 
0051 
0052     /// @name Run properties
0053     /// @{
0054 
0055     /// Get the number of events seen. Should only really be used by external
0056     /// steering code or analyses in the finalize phase.
0057     ///
0058     /// N.B. This only reports the count for the last collapsed event group
0059     /// and hence ignores any additional sub-events seen so far.
0060     size_t numEvents() const {
0061       const double N = _eventCounter.get()->persistent(defaultWeightIndex())->numEntries();
0062       return  size_t(N + 0.5 - (N<0)); // round to nearest integer
0063     }
0064 
0065     /// Get the effective number of events seen. Should only really be used by external
0066     /// steering code or analyses in the finalize phase.
0067     ///
0068     /// N.B. This only reports the count for the last collapsed event group
0069     /// and hence ignores any additional sub-events seen so far.
0070     double effNumEvents() const {
0071       if ((bool)_eventCounter) { return _eventCounter->effNumEntries(); }
0072       return _eventCounter.get()->persistent(defaultWeightIndex())->effNumEntries();
0073     }
0074 
0075     /// @brief Access the sum of the event weights seen
0076     ///
0077     /// This is the weighted equivalent of the number of events. It should only
0078     /// be used by external steering code or analyses in the finalize phase.
0079     double sumW() const {
0080       if ((bool)_eventCounter) { return _eventCounter->sumW(); }
0081       return _eventCounter.get()->persistent(defaultWeightIndex())->sumW();
0082     }
0083     /// Access to the sum of squared-weights
0084     double sumW2() const {
0085       if ((bool)_eventCounter) { return _eventCounter->sumW2(); }
0086       return _eventCounter.get()->persistent(defaultWeightIndex())->sumW2();
0087     }
0088 
0089     /// @}
0090 
0091 
0092     /// @name Event weights
0093     /// @{
0094 
0095     /// Names of event weight categories
0096     const vector<string>& weightNames() const { return _weightNames; }
0097 
0098     /// Are any of the weights non-numeric?
0099     size_t numWeights() const { return _weightNames.size(); }
0100 
0101     /// Are any of the weights non-numeric?
0102     bool haveNamedWeights() const;
0103 
0104     /// Set the weight names from a GenEvent
0105     void setWeightNames(const GenEvent& ge);
0106 
0107     /// Set the weight names from a vector<string>
0108     void setWeightNames(const vector<string>& weightNames);
0109 
0110     /// Get the index of the nominal weight-stream
0111     size_t defaultWeightIndex() const { return _rivetDefaultWeightIdx; }
0112 
0113     /// @brief Access the array of sum of the event weights seen
0114     vector<double> weightSumWs() const;
0115 
0116     /// Set the weight cap
0117     void setWeightCap(const double maxWeight) { _weightCap = maxWeight; }
0118 
0119     /// Set the name of the nominal weight stream
0120     void setNominalWeightName(const std::string& name) { _nominalWeightName = name; }
0121 
0122     /// Ignore all weight streams other than the nominal
0123     void skipMultiWeights(bool skip=false) { _skipMultiWeights = skip; }
0124 
0125     /// Specify weight-name patterns to accept
0126     void matchWeightNames(const std::string& patterns) { _matchWeightNames = patterns; }
0127 
0128     /// Specify weight-name patterns to reject
0129     void unmatchWeightNames(const std::string& patterns) { _unmatchWeightNames = patterns; }
0130 
0131     /// Set the relative width of the NLO smearing window.
0132     void setNLOSmearing(double frac) { _NLOSmearing = frac; }
0133 
0134     /// @}
0135 
0136 
0137     /// @name Cross-sections
0138     /// @{
0139 
0140     /// Get the cross-section known to the handler
0141     Estimate0DPtr crossSection() const { return _xs; }
0142 
0143     /// Set all cross-sections for the process being generated specifically (preferred)
0144     void setCrossSection(const vector<pair<double,double>>& xsecs, bool isUserSupplied = false);
0145 
0146     /// Set all cross-sections for the process being generated, based on nominal weight
0147     void setCrossSection(const pair<double, double>& xsec, bool isUserSupplied=false);
0148 
0149     /// Set the cross-section for the process being generated (alternative signature)
0150     void setCrossSection(double xsec, double xsecerr, bool isUserSupplied=false) {
0151       setCrossSection({xsec, xsecerr}, isUserSupplied);
0152     }
0153 
0154     /// Update the internal cross-section average when running over multiple files
0155     ///
0156     /// @note This method should only be called when switching HepMC file mid-run.
0157     void updateCrossSection();
0158 
0159     /// Toggle to signal a change in HepMC input file
0160     void notifyEndOfFile() { _isEndOfFile = true; }
0161 
0162     /// Get the nominal cross-section
0163     double nominalCrossSection() const;
0164 
0165     /// Get the nominal cross-section
0166     double nominalCrossSectionError() const;
0167 
0168     /// @}
0169 
0170 
0171     /// @name Beams
0172     /// @{
0173 
0174     /// Set the beam particles for this run
0175     AnalysisHandler& setRunBeams(const ParticlePair& beams);
0176 
0177     /// Get the beam particles for this run, usually determined from the first event
0178     const ParticlePair& runBeams() const { return _beams; }
0179 
0180     /// Get beam IDs for this run, usually determined from the first event
0181     PdgIdPair runBeamIDs() const;
0182 
0183     /// Get beam IDs for this run, usually determined from the first event
0184     pair<double,double> runBeamEnergies() const;
0185 
0186     /// Get energy for this run, usually determined from the first event
0187     double runSqrtS() const;
0188 
0189     /// Option to disable analysis-compatibility checks
0190     void setCheckBeams(bool check=true) { _checkBeams = check; }
0191 
0192     /// Option to disable run-consistency checks
0193     // void setCheckConsistency(bool check=true) { _checkConsistency = check; }
0194     // Check event consistency with the run, usually determined from the first event
0195     // bool consistentWithRun(Event& event) {
0196 
0197     /// @}
0198 
0199 
0200     /// @name Run-based annotations
0201     /// @{
0202 
0203     // Get all the annotation names
0204     std::vector<std::string> annotations() const {
0205       return _beaminfo->annotations();
0206     }
0207 
0208     /// Check if an annotation is defined
0209     bool hasAnnotation(const std::string& name) const {
0210       return _beaminfo->hasAnnotation(name);
0211     }
0212 
0213     /// Get an annotation by name (as a string)
0214     const std::string& annotation(const std::string& name) const {
0215       return _beaminfo->annotation(name);
0216     }
0217 
0218     /// Get an annotation by name (as a string) with a default in case the annotation is not found
0219     const std::string& annotation(const std::string& name, const std::string& defaultreturn) const {
0220       return _beaminfo->annotation(name, defaultreturn);
0221     }
0222 
0223     /// @brief Get an annotation by name (copied to another type)
0224     ///
0225     /// @note Templated on return type
0226     template <typename T>
0227     const T annotation(const std::string& name) const {
0228       return _beaminfo->annotation<T>(name);
0229     }
0230 
0231     /// @brief Get an annotation by name (copied to another type) with a default in case the annotation is not found
0232     ///
0233     /// @note Templated on return type
0234     template <typename T>
0235     const T annotation(const std::string& name, T&& defaultreturn) const {
0236       return _beaminfo->annotation<T>(name, std::forward<T>(defaultreturn));
0237     }
0238 
0239     /// @brief Add or set an annotation by name (templated for remaining types)
0240     ///
0241     /// @note Templated on arg type, but stored as a string.
0242     template <typename T>
0243     void setAnnotation(const std::string& name, T&& value) {
0244       _beaminfo->annotation<T>(name, std::forward<T>(value));
0245     }
0246 
0247 
0248     /// Set all annotations at once
0249     void setAnnotations(const Annotations& anns) {
0250       _beaminfo->setAnnotations(anns);
0251     }
0252 
0253     /// Delete an annotation by name
0254     void rmAnnotation(const std::string& name) {
0255       _beaminfo->rmAnnotation(name);
0256     }
0257 
0258 
0259     /// Delete an annotation by name
0260     void clearAnnotations() {
0261       _beaminfo->clearAnnotations();
0262     }
0263 
0264     /// @}
0265 
0266 
0267     /// @name AO type handling
0268     /// @{
0269 
0270     /// Register an AO type handle into type map and YODA reader
0271     template<typename T>
0272     void registerType() {
0273       const std::string name = T().type();
0274       const TypeRegisterItr& res = _register.find(name);
0275       if (res == _register.end()) {
0276         _register[name] = make_shared<TypeHandle<T>>();
0277       }
0278       _reader.registerType<T>(); // also let YODA know
0279     }
0280 
0281     /// If @a dst is the same subclass as @a src, copy the contents of
0282     /// @a src into @a dst and return true. Otherwise return false.
0283     bool copyAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst, const double scale=1.0);
0284 
0285     /// If @a dst is the same subclass as @a src, scale the contents of @a src with
0286     /// @a scale and add it to @a dst and return true. Otherwise return false.
0287     bool addAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr& dst, const double scale);
0288 
0289     /// @}
0290 
0291 
0292     /// @name Analysis handling
0293     /// @{
0294 
0295     /// Get a list of the currently registered analyses' names.
0296     std::vector<std::string> analysisNames() const;
0297 
0298     /// Get a list of the official analysis names for this release.
0299     std::vector<std::string> stdAnalysisNames() const;
0300 
0301     /// Get the collection of currently registered analyses.
0302     const std::map<std::string, AnaHandle>& analysesMap() const {
0303       return _analyses;
0304     }
0305 
0306     /// Get the collection of currently registered analyses.
0307     std::vector<AnaHandle> analyses() const {
0308       std::vector<AnaHandle> rtn;
0309       rtn.reserve(_analyses.size());
0310       for (const auto& apair : _analyses) rtn.push_back(apair.second);
0311       return rtn;
0312     }
0313 
0314     /// Get a registered analysis by name.
0315     AnaHandle analysis(const std::string& analysisname) {
0316       if ( _analyses.find(analysisname) == _analyses.end() )
0317         throw LookupError("No analysis named '" + analysisname + "' registered in AnalysisHandler");
0318       try {
0319         return _analyses[analysisname];
0320       } catch (...) {
0321         throw LookupError("No analysis named '" + analysisname + "' registered in AnalysisHandler");
0322       }
0323     }
0324 
0325     /// Add an analysis to the run list by object
0326     AnalysisHandler& addAnalysis(Analysis* analysis);
0327 
0328     /// @brief Add an analysis to the run list using its name.
0329     ///
0330     /// The actual Analysis to be used will be obtained via
0331     /// AnalysisLoader::getAnalysis(string).  If no matching analysis is found,
0332     /// no analysis is added (i.e. the null pointer is checked and discarded.
0333     AnalysisHandler& addAnalysis(const std::string& analysisname);
0334 
0335     /// @brief Add an analysis with a map of analysis options.
0336     AnalysisHandler& addAnalysis(const std::string& analysisname, std::map<string, string> pars);
0337 
0338     /// @brief Add analyses to the run list using their names.
0339     ///
0340     /// The actual {@link Analysis}' to be used will be obtained via
0341     /// AnalysisHandler::addAnalysis(string), which in turn uses
0342     /// AnalysisLoader::getAnalysis(string). If no matching analysis is found
0343     /// for a given name, no analysis is added, but also no error is thrown.
0344     AnalysisHandler& addAnalyses(const std::vector<std::string>& analysisnames);
0345 
0346 
0347     /// Remove an analysis from the run list using its name.
0348     AnalysisHandler& removeAnalysis(const std::string& analysisname);
0349 
0350     /// Remove analyses from the run list using their names.
0351     AnalysisHandler& removeAnalyses(const std::vector<std::string>& analysisnames);
0352 
0353     /// @}
0354 
0355 
0356     /// @name Main init/execute/finalise
0357     /// @{
0358 
0359     /// Initialize a run, with the run beams taken from the example event
0360     void init(const GenEvent& event);
0361 
0362     /// @brief Analyze the given @a event by reference
0363     ///
0364     /// This function will call the AnalysisBase::analyze() function of all
0365     /// included analysis objects.
0366     ///
0367     /// @note Despite the event being passed as const, its units etc. may be changed, hence non-const.
0368     void analyze(GenEvent& event);
0369 
0370     /// @brief Analyze the given @a event by pointer
0371     ///
0372     /// This function will call the AnalysisBase::analyze() function of all
0373     /// included analysis objects, after checking the event pointer validity.
0374     void analyze(GenEvent* event);
0375 
0376     /// Finalize a run. This function calls the AnalysisBase::finalize()
0377     /// functions of all included analysis objects.
0378     void finalize();
0379 
0380     /// @}
0381 
0382 
0383     /// @name Histogram / data object access
0384     /// @{
0385 
0386     /// After all subevents in an event group have been processed, push
0387     /// all histo fills to the relevant histograms.
0388     void collapseEventGroup();
0389 
0390     /// @brief Read analysis plots into the histo collection from the given stream
0391     ///
0392     /// Use the @a fmt flag to specify the YODA output format (yoda, yoda.gz, yoda.h5, ...)
0393     void readData(std::istream& istr, const string& fmt, bool preload = true);
0394 
0395     /// Read analysis plots into the histo collection (via addData) from the named file.
0396     void readData(const std::string& filename, bool preload = true);
0397 
0398     /// Get all YODA analysis objects (across all weights, optionally including RAW)
0399     ///
0400     /// @note We'll live with the mixed-case "Yoda" here, since the consistent all-caps would be worse!
0401     vector<YODA::AnalysisObjectPtr> getYodaAOs(const bool includeraw=false, const bool mkinert=true) const;
0402 
0403     /// Get all raw YODA analysis objects (across all weights)
0404     vector<YODA::AnalysisObjectPtr> getRawAOs() const;
0405 
0406     /// Get all raw YODA analysis object paths (across all weights)
0407     vector<std::string> getRawAOPaths() const;
0408 
0409     /// Get a pointer to a preloaded yoda object with the given path,
0410     /// or null if path is not found.
0411     const YODA::AnalysisObjectPtr getPreload(const string& path) const {
0412       auto it = _preloads.find(path);
0413       if ( it == _preloads.end() ) return nullptr;
0414       return it->second;
0415     }
0416 
0417     /// @brief Write all analyses' plots (via getData) to the given stream.
0418     ///
0419     /// Use the @a fmt flag to specify the YODA output format (yoda, yoda.gz, yoda.h5, ...)
0420     void writeData(std::ostream& ostr, const string& fmt) const;
0421 
0422     /// Write all analyses' plots (via getData) to the named file.
0423     void writeData(const string& filename) const;
0424 
0425     /// @brief Configure the AnalysisObject dump rate and destination.
0426     ///
0427     /// Tell Rivet to dump intermediate result to a file named @a
0428     /// dumpfile every @a period'th event. If @a period is not positive,
0429     /// no periodic finalization will be done.
0430     void setFinalizePeriod(const string& dumpfile, int period) {
0431       _dumpPeriod = period;
0432       _dumpFile = dumpfile;
0433     }
0434     /// @brief Configure the AnalysisObject dump rate and destination.
0435     void setNoFinalizePeriod() {
0436       setFinalizePeriod("DUMMY", -1);
0437     }
0438 
0439     /// @brief Set filename of the bootstrap file
0440     void setBootstrapFilename(const string& filename) {
0441       _bootstrapfilename = filename;
0442     }
0443 
0444     /// Return a vector of (AO path, AO numBins) pairs to decode the fills layout
0445     vector<pair<string,size_t>> fillLayout() const;
0446 
0447     /// Return a vector of the binary fill outcome (was/wasn't filled) at each fill position
0448     vector<bool> fillOutcomes() const;
0449 
0450     /// Return a vector of the fill fraction at each fill position
0451     vector<double> fillFractions() const;
0452 
0453     /// @brief Merge the vector of YODA files, using the cross-section and weight information provided in each.
0454     ///
0455     /// Each file in @a aofiles is assumed to have been produced by Rivet. By
0456     /// default the files are assumed to contain different processes (or the
0457     /// same processs but mutually exclusive cuts), but if @a equiv if true, the
0458     /// files are assumed to contain output of completely equivalent (but
0459     /// statistically independent) Rivet runs. The corresponding analyses will
0460     /// be loaded and their analysis objects will be filled with the merged
0461     /// result. finalize() will be run on each relevant analysis. The resulting
0462     /// YODA file can then be written out by writeData().
0463     ///
0464     /// If @a delopts is non-empty, it is assumed to contain names of different
0465     /// options to be merged into the same analysis objects.
0466     ///
0467 
0468     void mergeYODAs(const vector<string>& aofiles,
0469                     const vector<string>& delopts=vector<string>(),
0470                     const vector<string>& addopts=vector<string>(),
0471                     const vector<string>& matches=vector<string>(),
0472                     const vector<string>& unmatches=vector<string>(),
0473                     const bool equiv=false, const bool reentrantOnly = true);
0474 
0475     /// A method to merge another AnalysisHandler into the current one
0476     void merge(AnalysisHandler &other);
0477 
0478     /// @brief A method to prepare a re-entrant run for a given set
0479     /// of AO paths and serialized AO data
0480     ///
0481     /// The @a unscale parameter multiplies fillable objects with sumW/xsec to counteract
0482     /// the cross-section scaling in finalize() when merging different processes (non-equiv)
0483     void loadAOs(const vector<string>& aoPaths, const vector<double>& aoData);
0484 
0485     /// @}
0486 
0487     /// @name MPI (de-)serialisation
0488     ///@{
0489 
0490     vector<double> serializeContent(bool fixed_length = false) {
0491       if (!_initialised)
0492         throw Error("AnalysisHandler has not been initialised!");
0493 
0494       collapseEventGroup();
0495 
0496       // Loop over raw AOs and work out the size of the content data
0497       const vector<YODA::AnalysisObjectPtr> raos = getRawAOs();
0498       size_t total = 0;
0499       for (size_t i = 0; i < raos.size(); ++i) {
0500         total += raos[i]->lengthContent(fixed_length)+1;
0501       }
0502       total += _beaminfo->numBins()+1;
0503 
0504       // Loop over raw AOs and retrieve the content data
0505       std::vector<double> data; // serialized data vector
0506       data.reserve(total); // pre-allocate enough memory
0507       // Add beam IDs
0508       data.push_back(_beaminfo->numBins());
0509       for (const string& beamID : _beaminfo->xEdges()) {
0510         data.push_back(_beamInfoLabelToID(beamID));
0511       }
0512       // Add raw YODA AO content
0513       for (size_t i = 0; i < raos.size(); ++i) {
0514         vector<double> tmp = raos[i]->serializeContent(fixed_length);
0515         data.push_back(tmp.size()); // length of the AO
0516         data.insert(std::end(data),
0517                     std::make_move_iterator(std::begin(tmp)),
0518                     std::make_move_iterator(std::end(tmp)));
0519       }
0520       return data;
0521     }
0522 
0523     void deserializeContent(const vector<double>& data, size_t nprocs = 0) {
0524       if (!_initialised)
0525         throw Error("AnalysisHandler has not been initialised!");
0526 
0527       collapseEventGroup();
0528 
0529       // get Rivet AOs for access to raw AO pointers
0530       vector<MultiplexAOPtr> raos = getRivetAOs();
0531 
0532 
0533       // beam info first
0534       size_t iAO = 0, iW = 0, nBeams = data[0], offset = 1;
0535       if (nprocs)  nBeams /= nprocs;
0536       const auto itr = data.cbegin();
0537       // set beam IDs
0538       vector<int> edges{itr+offset, itr+offset+nBeams};
0539       vector<string> labels; labels.reserve(edges.size());
0540       size_t id = 0;
0541       for (int edge : edges) {
0542         if (nprocs >= 2)  edge /= nprocs;
0543         labels.push_back(_mkBeamInfoLabel(++id, edge));
0544       }
0545 
0546       _beaminfo = make_shared<YODA::BinnedEstimate<string>>(labels, "/TMP/_BEAMPZ");
0547       offset += nBeams;
0548       // set beam momenta
0549       size_t beamLen = *(itr + offset); ++offset;
0550       if (nprocs)  beamLen /= nprocs;
0551       std::vector<double> energies{itr+offset, itr+offset+beamLen};
0552       if (nprocs >= 2) {
0553         for (double& e : energies) { e /= nprocs; }
0554       }
0555       _beaminfo->deserializeContent(energies);
0556       for (auto& b : _beaminfo->bins(true))  b.rmErrs();
0557       offset += beamLen;
0558 
0559       // then the multiweighted AOs
0560       while (offset < data.size()) {
0561         if (iW < numWeights())  raos[iAO].get()->setActiveWeightIdx(iW);
0562         else {
0563           raos[iAO].get()->unsetActiveWeight();
0564           iW = 0; ++iAO; // move on to next AO
0565           raos[iAO].get()->setActiveWeightIdx(iW);
0566         }
0567 
0568         // obtain content length and set content iterators
0569         size_t aoLen = *(itr + offset); ++offset;
0570         if (nprocs)  aoLen /= nprocs;
0571         auto first = itr + offset;
0572         auto last = first + aoLen;
0573         // load data into AO
0574         raos[iAO].get()->activeAO()->deserializeContent(std::vector<double>{first, last});
0575 
0576         ++iW; offset += aoLen; // increment offset
0577       }
0578       raos[iAO].get()->unsetActiveWeight();
0579       // Reset cross-section bookkeeping
0580       _ntrials = 0.0;
0581       _fileCounter = CounterPtr(weightNames(), Counter("_FILECOUNT"));
0582       _xserr = CounterPtr(weightNames(), Counter("XSECERR"));
0583       if (nprocs >= 2) {
0584         for (size_t iW = 0; iW < numWeights(); ++iW) {
0585           *_fileCounter.get()->persistent(iW) = *_eventCounter.get()->persistent(iW);
0586           _xs.get()->persistent(iW)->scale(1.0/nprocs);
0587         }
0588       }
0589     }
0590 
0591     /// @}
0592 
0593 
0594     /// @name Processing stage
0595     /// @{
0596 
0597     /// Indicate which Rivet stage we're in.
0598     /// At the moment, only INIT is used to enable booking.
0599     enum class Stage { OTHER, INIT, FINALIZE };
0600 
0601     /// Return the current processing stage.
0602     Stage stage() const { return _stage; }
0603 
0604     /// @}
0605 
0606   private:
0607 
0608     /// @name Internal helper functions
0609     /// @{
0610 
0611     /// Get a logger object.
0612     Log& getLog() const;
0613 
0614     /// Get all multi-weight Rivet analysis object wrappers.
0615     vector<MultiplexAOPtr> getRivetAOs() const;
0616 
0617     /// Helper function to strip specific options from data object paths.
0618     void stripOptions(YODA::AnalysisObjectPtr ao, const vector<string>& delopts) const;
0619 
0620     /// @brief Merge the AO map @a newaos into @a allaos
0621     void mergeAOS(map<string, YODA::AnalysisObjectPtr> &allaos,
0622                   const map<string, YODA::AnalysisObjectPtr> &newaos,
0623                   map<string, std::array<double,4>> &allxsecs,
0624                   const vector<string>& delopts=vector<string>(),
0625                   const vector<string>& optAnas=vector<string>(),
0626                   const vector<string>& optKeys=vector<string>(),
0627                   const vector<string>& optVals=vector<string>(),
0628                   const bool equiv=false,
0629                   const bool overwrite_xsec = false,
0630                   const double user_xsec = 1.0);
0631 
0632 
0633     /// @brief A method to prepare a re-entrant run for a given set of analysis objects
0634     ///
0635     /// The @a unscale parameter multiplies fillable objects with sumW/xsec to counteract
0636     /// the cross-section scaling in finalize() when merging different processes (non-equiv)
0637     void loadAOs(const map<string, YODA::AnalysisObjectPtr>& allAOs,
0638                  const bool unscale = false, const bool reentrantOnly = true);
0639 
0640     /// @brief A method to set the internal _beaminfo object from a pair of beams
0641     void _setRunBeamInfo(const ParticlePair& beams);
0642 
0643     /// @brief A method to set the internal _beaminfo object from an existing YODA AO
0644     void _setRunBeamInfo(YODA::AnalysisObjectPtr ao);
0645 
0646     /// Construct beaminfo label from ID
0647     string _mkBeamInfoLabel(size_t n, PdgId id) {
0648       return "BEAM"+std::to_string(n)+"("+std::to_string(id)+")";
0649     }
0650 
0651     /// Retrieve ID from beaminfo label
0652     PdgId _beamInfoLabelToID(const string& label) {
0653       size_t pos = label.find("(");
0654       string beamID = label.substr(pos+1, label.size()-pos-2);
0655       return std::stoi(beamID);
0656     }
0657 
0658 
0659     /// @}
0660 
0661 
0662   private:
0663 
0664     /// Current handler processing stage.
0665     Stage _stage = Stage::OTHER;
0666 
0667     /// The collection of Analysis objects to be used.
0668     std::map<std::string, AnaHandle> _analyses;
0669 
0670     /// A vector of pre-loaded object which do not have a valid Analysis plugged in.
0671     ///
0672     /// @todo Rename to _preloadedAOs for consistency
0673     map<string,YODA::AnalysisObjectPtr> _preloads;
0674 
0675     /// A vector containing copies of analysis objects after finalize() has been run.
0676     vector<YODA::AnalysisObjectPtr> _finalizedAOs;
0677 
0678     /// Loads a default set of YODA AO types into the TypeHandleMap
0679     template<typename... Args>
0680     void registerDefaultTypes();
0681 
0682     /// The TypeRegister of all registered types
0683     TypeRegister _register;
0684 
0685     /// A reference to the ReaderYODA singleton
0686     YODA::Reader& _reader = YODA::ReaderYODA::create();
0687 
0688 
0689     /// @name Run properties
0690     /// @{
0691 
0692     /// Weight names
0693     std::vector<std::string> _weightNames;
0694     std::vector<std::valarray<double> > _subEventWeights;
0695     //size_t _numWeightTypes; // always == WeightVector.size()
0696 
0697     /// Weight indices
0698     std::vector<size_t> _weightIndices;
0699 
0700     /// Event counter
0701     CounterPtr _eventCounter;
0702 
0703     /// Cross-section known to AH
0704     Estimate0DPtr _xs;
0705 
0706     /// Running average of the cross-section uncertainties
0707     CounterPtr _xserr;
0708 
0709     /// Beam info known to AH
0710     YODA::BinnedEstimatePtr<string> _beaminfo;
0711 
0712     /// Total number of trials
0713     double _ntrials;
0714 
0715     /// Total number of entries of the current HepMC file
0716     CounterPtr _fileCounter;
0717 
0718     /// Toggle for multi-file runs
0719     bool _isEndOfFile;
0720 
0721     /// Nominal cross-section
0722     std::pair<double,double> _userxs;
0723 
0724     /// Beams used by this run.
0725     ParticlePair _beams;
0726 
0727     /// Flag to check if init has been called
0728     bool _initialised;
0729 
0730     /// Flag whether input event beams should be ignored in compatibility check
0731     bool _checkBeams;
0732 
0733     /// Flag to check if multiweights should be ignored
0734     bool _skipMultiWeights;
0735 
0736     /// String of weight names (or regex) to select multiweights
0737     std::string _matchWeightNames;
0738 
0739     /// String of weight names (or regex) to veto multiweights
0740     std::string _unmatchWeightNames;
0741 
0742     /// String giving the nominal weight name
0743     std::string _nominalWeightName;
0744 
0745     /// weight cap value
0746     double _weightCap;
0747 
0748     /// The relative width of the NLO smearing window.
0749     ///
0750     /// @todo Improve & standardise name
0751     double _NLOSmearing;
0752 
0753     /// Current event number
0754     int _eventNumber;
0755 
0756     /// The index in the (original) weight vector for the nominal weight stream
0757     size_t _defaultWeightIdx;
0758 
0759     /// The index in the (possibly pruned) weight vector for the nominal weight stream
0760     size_t _rivetDefaultWeightIdx;
0761 
0762     /// The index of the (possibly user-specified) intended default stream
0763     int _customDefaultWeightIdx;
0764 
0765     /// How often Rivet runs finalize() and writes the result to a YODA file.
0766     int _dumpPeriod;
0767 
0768     /// The name of a YODA file to which Rivet periodically dumps results.
0769     string _dumpFile;
0770 
0771     /// Flag to indicate periodic AO dumping is in progress
0772     bool _dumping;
0773 
0774     /// Flag to indicate bootstrap dumping is in progress
0775     ofstream _fbootstrap;
0776 
0777     /// Name of the file that the bootstrap dumping should be written to
0778     std::string _bootstrapfilename;
0779 
0780     /// Projection Handler
0781     ProjectionHandler _projHandler;
0782 
0783     /// @}
0784 
0785   };
0786 
0787 
0788 }
0789 
0790 #endif