diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh --- a/include/Rivet/AnalysisHandler.hh +++ b/include/Rivet/AnalysisHandler.hh @@ -1,358 +1,356 @@ // -*- C++ -*- #ifndef RIVET_RivetHandler_HH #define RIVET_RivetHandler_HH #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Particle.hh" #include "Rivet/AnalysisLoader.hh" #include "Rivet/Tools/RivetYODA.hh" namespace Rivet { // Forward declaration and smart pointer for Analysis class Analysis; typedef std::shared_ptr AnaHandle; // Needed to make smart pointers compare equivalent in the STL set struct CmpAnaHandle { bool operator() (const AnaHandle& a, const AnaHandle& b) const { return a.get() < b.get(); } }; /// A class which handles a number of analysis objects to be applied to /// generated events. An {@link Analysis}' AnalysisHandler is also responsible /// for handling the final writing-out of histograms. class AnalysisHandler { public: /// @name Constructors and destructors. */ //@{ /// Preferred constructor, with optional run name. AnalysisHandler(const string& runname=""); /// @brief Destructor /// The destructor is not virtual, as this class should not be inherited from. ~AnalysisHandler(); //@} private: /// Get a logger object. Log& getLog() const; public: - /// @name Run properties + /// @name Run properties and weights //@{ /// Get the name of this run. string runName() const; - /// Get the number of events seen. Should only really be used by external /// steering code or analyses in the finalize phase. size_t numEvents() const; /// @brief Access the sum of the event weights seen /// /// This is the weighted equivalent of the number of events. It should only /// be used by external steering code or analyses in the finalize phase. double sumW() const { return _eventCounter->sumW(); } /// Access to the sum of squared-weights double sumW2() const { return _eventCounter->sumW2(); } /// Names of event weight categories - const vector& weightNames() const { - return _weightNames; - } + const vector& weightNames() const { return _weightNames; } + /// Are any of the weights non-numeric? + size_t numWeights() const { return _weightNames.size(); } + + /// Are any of the weights non-numeric? + bool haveNamedWeights() const; + + /// Set the weight names from a GenEvent + void setWeightNames(const GenEvent& ge); /// Get the index of the nominal weight-stream - size_t defaultWeightIndex() const { - return _defaultWeightIdx; - } + size_t defaultWeightIndex() const { return _defaultWeightIdx; } + //@} + + + /// @name Cross-sections + //@{ + + /// Get the cross-section known to the handler. + Scatter1DPtr crossSection() const { return _xs; } /// Set the cross-section for the process being generated. AnalysisHandler& setCrossSection(double xs, double xserr); - /// Get the cross-section known to the handler. - Scatter1DPtr crossSection() const { return _xs; } - /// Get the nominal cross-section double nominalCrossSection() const { _xs.get()->setActiveWeightIdx(_defaultWeightIdx); const YODA::Scatter1D::Points& ps = _xs->points(); if (ps.size() != 1) { string errMsg = "cross section missing when requesting nominal cross section"; throw Error(errMsg); } double xs = ps[0].x(); _xs.get()->unsetActiveWeight(); return xs; } + //@} + + + /// @name Beams + //@{ /// Set the beam particles for this run AnalysisHandler& setRunBeams(const ParticlePair& beams) { _beams = beams; MSG_DEBUG("Setting run beams = " << beams << " @ " << sqrtS()/GeV << " GeV"); return *this; } /// Get the beam particles for this run, usually determined from the first event. const ParticlePair& beams() const { return _beams; } /// Get beam IDs for this run, usually determined from the first event. /// @deprecated Use standalone beamIds(ah.beams()), to clean AH interface PdgIdPair beamIds() const; /// Get energy for this run, usually determined from the first event. /// @deprecated Use standalone sqrtS(ah.beams()), to clean AH interface double sqrtS() const; /// Setter for _ignoreBeams void setIgnoreBeams(bool ignore=true); //@} /// @name Handle analyses //@{ /// Get a list of the currently registered analyses' names. std::vector analysisNames() const; /// Get the collection of currently registered analyses. const std::set& analyses() const { return _analyses; } /// Get a registered analysis by name. const AnaHandle analysis(const std::string& analysisname) const; /// Add an analysis to the run list by object AnalysisHandler& addAnalysis(Analysis* analysis); /// @brief Add an analysis to the run list using its name. /// /// The actual Analysis to be used will be obtained via /// AnalysisLoader::getAnalysis(string). If no matching analysis is found, /// no analysis is added (i.e. the null pointer is checked and discarded. AnalysisHandler& addAnalysis(const std::string& analysisname); /// @brief Add an analysis with a map of analysis options. AnalysisHandler& addAnalysis(const std::string& analysisname, std::map pars); /// @brief Add analyses to the run list using their names. /// /// The actual {@link Analysis}' to be used will be obtained via /// AnalysisHandler::addAnalysis(string), which in turn uses /// AnalysisLoader::getAnalysis(string). If no matching analysis is found /// for a given name, no analysis is added, but also no error is thrown. AnalysisHandler& addAnalyses(const std::vector& analysisnames); /// Remove an analysis from the run list using its name. AnalysisHandler& removeAnalysis(const std::string& analysisname); /// Remove analyses from the run list using their names. AnalysisHandler& removeAnalyses(const std::vector& analysisnames); //@} /// @name Main init/execute/finalise //@{ /// Initialize a run, with the run beams taken from the example event. void init(const GenEvent& event); /// @brief Analyze the given \a event by reference. /// /// This function will call the AnalysisBase::analyze() function of all /// included analysis objects. void analyze(const GenEvent& event); /// @brief Analyze the given \a event by pointer. /// /// This function will call the AnalysisBase::analyze() function of all /// included analysis objects, after checking the event pointer validity. void analyze(const GenEvent* event); /// Finalize a run. This function calls the AnalysisBase::finalize() /// functions of all included analysis objects. void finalize(); //@} /// @name Histogram / data object access //@{ /// Add a vector of analysis objects to the current state. void addData(const std::vector& aos); /// Read analysis plots into the histo collection (via addData) from the named file. void readData(const std::string& filename); /// Get all multi-weight Rivet analysis object wrappers vector getRivetAOs() const; /// Get single-weight YODA analysis objects vector getYodaAOs(bool includeorphans, bool includetmps, bool usefinalized) const; /// Get all analyses' plots as a vector of analysis objects. std::vector getData(bool includeorphans = false, bool includetmps = false, bool usefinalized = true) const; /// Write all analyses' plots (via getData) to the named file. void writeData(const std::string& filename) const; /// Tell Rivet to dump intermediate result to a file named @a /// dumpfile every @a period'th event. If @period is not positive, /// no dumping will be done. void dump(string dumpfile, int period) { _dumpPeriod = period; _dumpFile = dumpfile; } /// Take the vector of yoda files and merge them together using /// the cross section and weight information provided in each /// file. Each file in @a aofiles is assumed to have been produced /// by Rivet. By default the files are assumed to contain /// different processes (or the same processs but mutually /// exclusive cuts), but if @a equiv if ture, the files are /// assumed to contain output of completely equivalent (but /// statistically independent) Rivet runs. The corresponding /// analyses will be loaded and their analysis objects will be /// filled with the merged result. finalize() will be run on each /// relevant analysis. The resulting YODA file can then be rwitten /// out by writeData(). If delopts is non-empty, it is assumed to /// contain names different options to be merged into the same /// analysis objects. void mergeYodas(const vector & aofiles, const vector & delopts = vector(), bool equiv = false); /// Helper function to strip specific options from data object paths. void stripOptions(YODA::AnalysisObjectPtr ao, const vector & delopts) const; //@} - /// @name Event weight access - //@{ - - /// Are any of the weights non-numeric? - size_t numWeights() { return _weightNames.size(); } - - /// Are any of the weights non-numeric? - bool haveNamedWeights(); - - /// Set the weight names from a GenEvent - void setWeightNames(const GenEvent& ge); - - //@} - - /// Indicate which Rivet stage we're in. /// At the moment, only INIT is used to enable booking. enum class Stage { OTHER, INIT }; /// Which stage are we in? Stage stage() const { return _stage; } private: /// Current handler stage Stage _stage = Stage::OTHER; /// The collection of Analysis objects to be used. set _analyses; /// A vector of pre-loaded object which do not have a valid /// Analysis plugged in. vector _orphanedPreloads; /// A vector containing copies of analysis objects after /// finalize() has been run. vector _finalizedAOs; /// @name Run properties //@{ /// Weight names std::vector _weightNames; std::vector > _subEventWeights; size_t _numWeightTypes; // always == WeightVector.size() /// Run name std::string _runname; /// Event counter mutable CounterPtr _eventCounter; /// Cross-section known to AH Scatter1DPtr _xs; /// Beams used by this run. ParticlePair _beams; /// Flag to check if init has been called bool _initialised; /// Flag whether input event beams should be ignored in compatibility check bool _ignoreBeams; /// Current event number int _eventNumber; /// The index in the weight vector for the nominal weight stream size_t _defaultWeightIdx; /// Determines how often Rivet runs finalize() and writes the /// result to a YODA file. int _dumpPeriod; /// The name of a YODA file to which Rivet periodically dumps /// results. string _dumpFile; /// Flag to indicate periodic dumping is in progress bool _dumping; //@} private: /// The assignment operator is private and must never be called. /// In fact, it should not even be implemented. AnalysisHandler& operator=(const AnalysisHandler&); /// The copy constructor is private and must never be called. In /// fact, it should not even be implemented. AnalysisHandler(const AnalysisHandler&); }; } #endif diff --git a/include/Rivet/Tools/Correlators.hh b/include/Rivet/Tools/Correlators.hh --- a/include/Rivet/Tools/Correlators.hh +++ b/include/Rivet/Tools/Correlators.hh @@ -1,1136 +1,1147 @@ // -*- C++ -*- #ifndef RIVET_Correlators_HH #define RIVET_Correlators_HH #include "Rivet/Projection.hh" #include "Rivet/Projections/ParticleFinder.hh" #include #include "YODA/Scatter2D.h" #include "Rivet/Analysis.hh" -/* File contains tools for calculating flow coefficients using +/* File contains tools for calculating flow coefficients using * correlators. * Classes: * Correlators: Calculates single event correlators of a given harmonic. - * Cumulants: An additional base class for flow analyses + * Cumulants: An additional base class for flow analyses * (use as: class MyAnalysis : public Analysis, Cumulants {};) * Includes a framework for calculating cumulants and flow coefficients * from single event correlators, including automatic handling of statistical * errors. Contains multiple internal sub-classes: * CorBinBase: Base class for correlators binned in event or particle observables. * CorSingleBin: A simple bin for correlators. * CorBin: Has the interface of a simple bin, but does automatic calculation * of statistical errors by a bootstrap method. - * ECorrelator: Data type for event averaged correlators. + * ECorrelator: Data type for event averaged correlators. * @author Vytautas Vislavicius, Christine O. Rasmussen, Christian Bierlich. */ namespace Rivet { + using std::complex; + + /* @brief Projection for calculating correlators for flow measurements. - * + * * A projection which calculates Q-vectors and P-vectors, and projects - * them out as correlators. Implementation follows the description of - * the ''Generic Framework'' + * them out as correlators. Implementation follows the description of + * the ''Generic Framework'' * Phys. Rev. C 83 (2011) 044913, arXiv: 1010.0233 * Phys. Rev. C 89 (2014) 064904, arXiv: 1312.3572 - * + * */ class Correlators : public Projection { - + public: // Constructor. Takes parameters @parm fsp, the Final State - // projection correlators should be constructed from, @parm nMaxIn, - // the maximal sum of harmonics, eg. for + // projection correlators should be constructed from, @parm nMaxIn, + // the maximal sum of harmonics, eg. for // c_2{2} = {2,-2} = 2 + 2 = 4 // c_2{4} = {2,2,-2,-2} = 2 + 2 + 2 + 2 = 8 // c_4{2} = {4,-4} = 4 + 4 = 8 - // c_4{4} = {4,4,-4,-4} = 4 + 4 + 4 + 4 = 16. + // c_4{4} = {4,4,-4,-4} = 4 + 4 + 4 + 4 = 16. // @parm pMaxIn is the maximal number of particles - // you want to correlate and @parm pTbinEdgesIn is the (lower) + // you want to correlate and @parm pTbinEdgesIn is the (lower) // edges of pT bins, the last one the upper edge of the final bin. - Correlators(const ParticleFinder& fsp, int nMaxIn = 2, + Correlators(const ParticleFinder& fsp, int nMaxIn = 2, int pMaxIn = 0, vector pTbinEdgesIn = {}); - + // Constructor which takes a Scatter2D to estimate bin edges. - Correlators(const ParticleFinder& fsp, int nMaxIn, + Correlators(const ParticleFinder& fsp, int nMaxIn, int pMaxIn, const Scatter2DPtr hIn); - + /// @brief Integrated correlator of @parm n harmonic, with the /// number of powers being the size of @parm n. /// Eg. @parm n should be: /// <<2>>_2 => n = {2, -2} /// <<4>>_2 => n = {2, 2, -2, -2} /// <<2>>_4 => n = {4, -4} /// <<4>>_4 => n = {4, 4, -4, 4} and so on. const pair intCorrelator(vector n) const; /// @brief pT differential correlator of @parm n harmonic, with the /// number of powers being the size of @parm n. - /// The method can include overflow/underflow bins in the + /// The method can include overflow/underflow bins in the /// beginning/end of the returned vector, by toggling /// @parm overflow = true. const vector> pTBinnedCorrelators(vector n, bool overflow = false) const; /// @brief Integrated correlator of @parm n1 harmonic, with the /// number of powers being the size of @parm n1. This method - /// imposes an eta gap, correlating with another phase space, - /// where another Correlators projection @parm other should be - /// defined. The harmonics of the other phase space is given + /// imposes an eta gap, correlating with another phase space, + /// where another Correlators projection @parm other should be + /// defined. The harmonics of the other phase space is given /// as @parm n2. /// To get eg. integrated <<4>>_2, @parm n1 should be: /// n1 = {2, 2} and n2 = {-2, -2} - const pair intCorrelatorGap(const Correlators& other, + const pair intCorrelatorGap(const Correlators& other, vector n1, vector n2) const; /// @brief pT differential correlators of @parm n1 harmonic, with the /// number of powers being the size of @parm n1. This method - /// imposes an eta gap, correlating with another phase space, - /// where another Correlators projection @parm other should be - /// defined. The harmonics of the other phase space is given + /// imposes an eta gap, correlating with another phase space, + /// where another Correlators projection @parm other should be + /// defined. The harmonics of the other phase space is given /// as @parm n2. /// To get eg. differential <<4'>_2, @parm n1 should be: /// n1 = {2, 2} and @parm n2: n2 = {-2, -2}. /// To get eg. differential <<2'>>_4, @parm n1 should be: /// n1 = {4} and @parm n2: n2 = {-4}. /// The method can include overflow/underflow /// bins in the beginning/end of the returned vector, by toggling /// @parm overflow = true. const vector> pTBinnedCorrelatorsGap( - const Correlators& other, vector n1, vector n2, - bool oveflow = false) const; - + const Correlators& other, vector n1, vector n2, + bool oveflow = false) const; + /// @brief Construct a harmonic vectors from @parm n harmonics /// and @parm m number of particles. /// TODO: In C++14 this can be done much nicer with TMP. static vector hVec(int n, int m) { if (m%2 != 0) { cout << "Harmonic Vector: Number of particles " "must be an even number." << endl; return {}; } vector ret = {}; for (int i = 0; i < m; ++i) { if (i < m/2) ret.push_back(n); else ret.push_back(-n); } return ret; } /// @brief Return the maximal values for n, p to be used in the /// constructor of Correlators(xxx, nMax, pMax, xxxx) static pair getMaxValues(vector< vector >& hList){ int nMax = 0, pMax = 0; for (vector h : hList) { int tmpN = 0, tmpP = 0; for ( int i = 0; i < int(h.size()); ++i) { tmpN += abs(h[i]); ++tmpP; } if (tmpN > nMax) nMax = tmpN; if (tmpP > pMax) pMax = tmpP; } return make_pair(nMax,pMax); } // Clone on the heap. DEFAULT_RIVET_PROJ_CLONE(Correlators); - + protected: - + // @brief Project function. Loops over array and calculates Q vectors - // and P vectors if needed - void project(const Event& e); - + // and P vectors if needed + void project(const Event& e); + // @brief Compare against other projection. Test if harmonics, pT-bins // and underlying final state are similar. - int compare(const Projection& p) const { + CmpState compare(const Projection& p) const { const Correlators* other = dynamic_cast(&p); - if (nMax != other->nMax) return UNDEFINED; - if (pMax != other->pMax) return UNDEFINED; - if (pTbinEdges != other->pTbinEdges) return UNDEFINED; + if (nMax != other->nMax) return CmpState::NEQ; + if (pMax != other->pMax) return CmpState::NEQ; + if (pTbinEdges != other->pTbinEdges) return CmpState::NEQ; return mkPCmp(p, "FS"); - }; + }; - // @brief Calculate correlators from one particle. - void fillCorrelators(const Particle& p, const double& weight); - + // @brief Calculate correlators from one particle. + void fillCorrelators(const Particle& p, const double& weight); + // @brief Return a Q-vector. const complex getQ(int n, int p) const { bool isNeg = (n < 0); if (isNeg) return conj( qVec[abs(n)][p] ); else return qVec[n][p]; }; - - // @brief Return a P-vector. + + // @brief Return a P-vector. const complex getP(int n, int p, double pT = 0.) const { bool isNeg = (n < 0); map::const_iterator pTitr = pVec.lower_bound(pT); if (pTitr == pVec.end()) return DBL_NAN; if (isNeg) return conj( pTitr->second[abs(n)][p] ); else return pTitr->second[n][p]; }; private: - // Find correlators by recursion. Order = M (# of particles), + // Find correlators by recursion. Order = M (# of particles), // n's are harmonics, p's are the powers of the weights - const complex recCorr(int order, vector n, - vector p, bool useP, double pT = 0.) const; + const complex recCorr(int order, vector n, + vector p, bool useP, double pT = 0.) const; // Two-particle correlator Eq. (19) p. 6 // Flag if p-vectors or q-vectors should be used to // calculate the correlator. - const complex twoPartCorr(int n1, int n2, int p1 = 1, + const complex twoPartCorr(int n1, int n2, int p1 = 1, int p2 = 1, double pT = 0., bool useP = false) const; // Set elements in vectors to zero. void setToZero(); - + // Shorthands for setting and comparing to zero. const complex _ZERO = {0., 0.}; const double _TINY = 1e-10; // Shorthand typedefs for vec>. typedef vector< vector> > Vec2D; // Define Q-vectors and p-vectors - Vec2D qVec; // Q[n][p] + Vec2D qVec; // Q[n][p] map pVec; // p[pT][n][p] - + // The max values of n and p to be calculated. int nMax, pMax; // pT bin-edges vector pTbinEdges; bool isPtDiff; - + /// End class Correlators }; /// @brief Tools for flow analyses. /// The following are helper classes to construct event averaged correlators - /// as well as cummulants and flow coefficents from the basic event - // correlators defined above. They are all encapsulated in a Cumulants - // class, which can be used as a(nother) base class for flow analyses, + /// as well as cummulants and flow coefficents from the basic event + // correlators defined above. They are all encapsulated in a Cumulants + // class, which can be used as a(nother) base class for flow analyses, // to ensure access. - + class CumulantAnalysis : public Analysis { private: - // Number of bins used for bootstrap calculation of statistical + // Number of bins used for bootstrap calculation of statistical // uncertainties. It is hard coded, and shout NOT be changed unless there - // are good reasons to do so. + // are good reasons to do so. static const int BOOT_BINS = 9; // Enum for choosing the method of error calculation. enum Error { VARIANCE, ENVELOPE }; - + // The desired error method. Can be changed in the analysis constructor - // by setting it appropriately. + // by setting it appropriately. Error errorMethod; /// @brief Base class for correlator bins. class CorBinBase { public: CorBinBase() {} virtual ~CorBinBase() {}; // Derived class should have fill and mean defined. virtual void fill(const pair& cor, const double& weight) = 0; virtual const double mean() const = 0; }; + /// @brief The CorSingleBin is the basic quantity filled in an ECorrelator. /// It is a simple counter with an even simpler structure than normal /// YODA type DBNs, but added functionality to test for out of /// bounds correlators. class CorSingleBin : public CorBinBase { - public: + /// @brief The default constructor. CorSingleBin() : _sumWX(0.), _sumW(0.), _sumW2(0.), _numEntries(0.) {} - + ~CorSingleBin() {} - /// @brief Fill a correlator bin with the return type from a + /// @brief Fill a correlator bin with the return type from a /// Correlator (a pair giving numerator and denominator of _event). void fill(const pair& cor, const double& weight) { // Test if denominator for the single event average is zero. if (cor.second < 1e-10) return; // The single event average is then cor.first / cor.second. // With weights this becomes just: _sumWX += cor.first * weight; _sumW += weight * cor.second; _sumW2 += weight * weight * cor.second * cor.second; _numEntries += 1.; } const double mean() const { if (_sumW < 1e-10) return 0; return _sumWX / _sumW; } // @brief Sum of weights. const double sumW() const { return _sumW; } const double sumW2() const { return _sumW2; } // @brief Sum of weight * X. const double sumWX() const { return _sumWX; } // @brief Number of entries. const double numEntries() const { return _numEntries; } - + void addContent(double ne, double sw, double sw2, double swx) { _numEntries += ne; _sumW += sw; _sumW2 += sw2; _sumWX += swx; } private: double _sumWX, _sumW, _sumW2, _numEntries; - + }; // End of CorSingleBin sub-class. /// @brief The CorBin is the basic bin quantity in ECorrelators. - /// It consists of several CorSingleBins, to facilitate + /// It consists of several CorSingleBins, to facilitate /// bootstrapping calculation of statistical uncertainties. class CorBin : public CorBinBase { public: // @brief The constructor. nBins signifies the period of the bootstrap // calculation, and should never be changed here, but in its definition // above -- and only if there are good reasons to do so. CorBin() : binIndex(0), nBins(BOOT_BINS) { - for(size_t i = 0; i < nBins; ++i) bins.push_back(CorSingleBin()); + for(size_t i = 0; i < nBins; ++i) bins.push_back(CorSingleBin()); } // Destructor must be implemented. ~CorBin() {} // @brief Fill the correct underlying bin and take a step. void fill(const pair& cor, const double& weight) { // Test if denominator for the single event average is zero. if (cor.second < 1e-10) return; // Fill the correct bin. bins[binIndex].fill(cor, weight); if (binIndex == nBins - 1) binIndex = 0; else ++binIndex; } // @brief Calculate the total sample mean with all // available statistics. const double mean() const { double sow = 0; double sowx = 0; for(auto b : bins) { if (b.sumW() < 1e-10) continue; sow += b.sumW(); sowx += b.sumWX(); } return sowx / sow; } // @brief Return a copy of the bins. const vector getBins() const { return bins; } // @brief Return the bins as pointers to the base class. template const vector getBinPtrs() { vector ret(bins.size()); transform(bins.begin(), bins.end(), ret.begin(), [](CorSingleBin& b) {return &b;}); return ret; } private: vector bins; size_t binIndex; - size_t nBins; + size_t nBins; }; // End of CorBin sub-class. - + public: - /// @brief The ECorrelator is a helper class to calculate all event + /// @brief The ECorrelator is a helper class to calculate all event /// averages of correlators, in order to construct cumulants. /// It can be binned in any variable. class ECorrelator { - + public: /// @brief Constructor. Takes as argument the desired harmonic and number /// of correlated particles as a generic framework style vector, eg, /// {2, -2} for <<2>>_2, no binning. - /// TODO: Implement functionality for this if needed. + /// TODO: Implement functionality for this if needed. //ECorrelator(vector h) : h1(h), h2({}), - // binX(0), binContent(0), reference() { + // binX(0), binContent(0), reference() { //}; /// @brief Constructor. Takes as argument the desired harmonic and number /// of correlated particles as a generic framework style vector, eg, /// {2, -2} for <<2>>_2 and binning. ECorrelator(vector h, vector binIn) : h1(h), h2({}), binX(binIn), binContent(binIn.size() - 1), reference() {}; - /// @brief Constructor for gapped correlator. Takes as argument the - /// desired harmonics for the two final states, and binning. - ECorrelator(vector h1In, vector h2In, vector binIn) : - h1(h1In), h2(h2In), binX(binIn), binContent(binIn.size() - 1), + /// @brief Constructor for gapped correlator. Takes as argument the + /// desired harmonics for the two final states, and binning. + ECorrelator(vector h1In, vector h2In, vector binIn) : + h1(h1In), h2(h2In), binX(binIn), binContent(binIn.size() - 1), reference() {}; - + /// @brief Fill the appropriate bin given an input (per event) /// observable, eg. centrality. - void fill(const double& obs, const Correlators& c, + void fill(const double& obs, const Correlators& c, const double weight = 1.0) { int index = getBinIndex(obs); if (index < 0) return; binContent[index].fill(c.intCorrelator(h1), weight); } /// @brief Fill the appropriate bin given an input (per event) /// observable, eg. centrality. Using a rapidity gap between /// two Correlators. void fill(const double& obs, const Correlators& c1, const Correlators& c2, const double weight = 1.0) { if (!h2.size()) { cout << "Trying to fill gapped correlator, but harmonics behind " "the gap (h2) are not given!" << endl; return; } int index = getBinIndex(obs); if (index < 0) return; binContent[index].fill(c1.intCorrelatorGap(c2, h1, h2), weight); } - /// @brief Fill the bins with the appropriate correlator, taking the - /// binning directly from the Correlators object, and filling also the + /// @brief Fill the bins with the appropriate correlator, taking the + /// binning directly from the Correlators object, and filling also the /// reference flow. void fill(const Correlators& c, const double& weight = 1.0) { vector< pair > diffCorr = c.pTBinnedCorrelators(h1); // We always skip overflow when calculating the all event average. if (diffCorr.size() != binX.size() - 1) cout << "Tried to fill event with wrong binning (ungapped)" << endl; for (size_t i = 0; i < diffCorr.size(); ++i) { int index = getBinIndex(binX[i]); if (index < 0) return; binContent[index].fill(diffCorr[i], weight); } reference.fill(c.intCorrelator(h1), weight); } /// @brief Fill bins with the appropriate correlator, taking the binning /// directly from the Correlators object, and also the reference flow. /// Using a rapidity gap between two Correlators. - void fill(const Correlators& c1, const Correlators& c2, + void fill(const Correlators& c1, const Correlators& c2, const double& weight = 1.0) { if (!h2.size()) { cout << "Trying to fill gapped correlator, but harmonics behind " "the gap (h2) are not given!" << endl; return; } vector< pair > diffCorr = c1.pTBinnedCorrelatorsGap(c2, h1, h2); // We always skip overflow when calculating the all event average. if (diffCorr.size() != binX.size() - 1) cout << "Tried to fill event with wrong binning (gapped)" << endl; for (size_t i = 0; i < diffCorr.size(); ++i) { int index = getBinIndex(binX[i]); if (index < 0) return; binContent[index].fill(diffCorr[i], weight); } reference.fill(c1.intCorrelatorGap(c2, h1, h2), weight); } - + /// @brief Get a copy of the bin contents. const vector getBins() const { return binContent; } // @brief Return the bins as pointers to the base class. const vector getBinPtrs() { vector ret(binContent.size()); transform(binContent.begin(), binContent.end(), ret.begin(), [](CorBin& b) {return &b;}); return ret; } /// @brief Get a copy of the bin x-values. const vector getBinX() const { return binX; } /// @brief Get a copy of the @ret h1 harmonic vector. const vector getH1() const { return h1; } /// @brief Get a copy of the @ret h2 harmonic vector. const vector getH2() const { return h2; } /// @brief Replace reference flow bin with another reference /// flow bin, eg. calculated in another phase space or with /// other pid. void setReference(CorBin refIn) { reference = refIn; } /// @brief Extract the reference flow from a differential event /// averaged correlator. const CorBin getReference() const { if (reference.mean() < 1e-10) cout << "Warning: ECorrelator, reference bin is zero." << endl; return reference; } /// @brief set the @parm prIn list of profile histograms associated /// with the internal bins. Is called automatically when booking, no /// need to call it yourself. void setProfs(list prIn) { profs = prIn; } /// @brief Fill bins with content from preloaded histograms. void fillFromProfs() { list::iterator hItr = profs.begin(); auto refs = reference.getBinPtrs(); for (size_t i = 0; i < profs.size(); ++i, ++hItr) { for (size_t j = 0; j < binX.size() - 1; ++j) { const YODA::ProfileBin1D& pBin = (*hItr)->binAt(binX[j]); auto tmp = binContent[j].getBinPtrs(); tmp[i]->addContent(pBin.numEntries(), pBin.sumW(), pBin.sumW2(), pBin.sumWY()); } // Get the reference flow from the underflow bin of the histogram. const YODA::Dbn2D& uBin = (*hItr)->underflow(); refs[i]->addContent(uBin.numEntries(), uBin.sumW(), uBin.sumW2(), uBin.sumWY()); } // End loop of bootstrapped correlators. - + } /// @brief begin() iterator for the list of associated profile histograms. list::iterator profBegin() { return profs.begin(); - } - + } + /// @brief end() iterator for the list of associated profile histograms. list::iterator profEnd() { return profs.end(); - } + } private: /// @brief Get correct bin index for a given @parm obs value. const int getBinIndex(const double& obs) const { // Find the correct index of binContent. // If we are in overflow, just skip. if (obs >= binX.back()) return -1; // If we are in underflow, ditto. if (obs < binX[0]) return -1; int index = 0; for (int i = 0, N = binX.size() - 1; i < N; ++i, ++index) if (obs >= binX[i] && obs < binX[i + 1]) break; return index; } // The harmonics vectors. vector h1; vector h2; // The bins. vector binX; vector binContent; // The reference flow. CorBin reference; // The profile histograms associated with the CorBins for streaming. list profs; }; // End of ECorrelator sub-class. - // @brief Get the correct max N and max P for the set of + // @brief Get the correct max N and max P for the set of // booked correlators. const pair getMaxValues() const { vector< vector> harmVecs; for ( auto eItr = eCorrPtrs.begin(); eItr != eCorrPtrs.end(); ++eItr) { vector h1 = (*eItr)->getH1(); vector h2 = (*eItr)->getH2(); if (h1.size() > 0) harmVecs.push_back(h1); if (h2.size() > 0) harmVecs.push_back(h2); } if (harmVecs.size() == 0) { cout << "Warning: You tried to extract max values from harmonic " "vectors, but have not booked any." << endl; return pair(); } return Correlators::getMaxValues(harmVecs); } // Typedeffing shared pointer to ECorrelator. typedef shared_ptr ECorrPtr; - - // @brief Book an ECorrelator in the same way as a histogram. + + /// @brief Book an ECorrelator in the same way as a histogram. + /// @todo Rename to book(ECorrPtr, ...) ECorrPtr bookECorrelator(const string name, const vector& h, const Scatter2DPtr hIn) { vector binIn; for (auto b : hIn->points()) binIn.push_back(b.xMin()); binIn.push_back(hIn->points().back().xMax()); ECorrPtr ecPtr = ECorrPtr(new ECorrelator(h, binIn)); list eCorrProfs; for (int i = 0; i < BOOT_BINS; ++i) { - Profile1DPtr tmp = bookProfile1D(name+"-"+to_string(i),*hIn); + Profile1DPtr tmp; + book(tmp, name+"-"+to_string(i),*hIn); tmp->setPath(this->name()+"/FINAL/" + name+"-"+to_string(i)); //tmp->setPath(tmp->path()+"FINAL"); eCorrProfs.push_back(tmp); } ecPtr->setProfs(eCorrProfs); eCorrPtrs.push_back(ecPtr); return ecPtr; } - // @brief Book a gapped ECorrelator with two harmonic vectors. - ECorrPtr bookECorrelator(const string name, const vector& h1, + /// @brief Book a gapped ECorrelator with two harmonic vectors. + /// @todo Rename to book(ECorrPtr, ...) + ECorrPtr bookECorrelator(const string name, const vector& h1, const vector& h2, const Scatter2DPtr hIn ) { vector binIn; for (auto b : hIn->points()) binIn.push_back(b.xMin()); binIn.push_back(hIn->points().back().xMax()); ECorrPtr ecPtr = ECorrPtr(new ECorrelator(h1, h2, binIn)); list eCorrProfs; for (int i = 0; i < BOOT_BINS; ++i) { - Profile1DPtr tmp = bookProfile1D(name+"-"+to_string(i),*hIn); + Profile1DPtr tmp; + book(tmp, name+"-"+to_string(i),*hIn); tmp->setPath(this->name()+"/FINAL/" + name+"-"+to_string(i)); //tmp->setPath(tmp->path()+"FINAL"); eCorrProfs.push_back(tmp); } ecPtr->setProfs(eCorrProfs); eCorrPtrs.push_back(ecPtr); return ecPtr; } - - // @brief Short hand for gapped correlators which splits the harmonic - // vector into negative and positive components. - ECorrPtr bookECorrelatorGap (const string name, const vector& h, + + /// @brief Short hand for gapped correlators which splits the harmonic + /// vector into negative and positive components. + /// @todo Rename to book(ECorrPtr, ...) + ECorrPtr bookECorrelatorGap (const string name, const vector& h, const Scatter2DPtr hIn) { const vector h1(h.begin(), h.begin() + h.size() / 2); const vector h2(h.begin() + h.size() / 2, h.end()); return bookECorrelator(name, h1, h2, hIn); } - // @brief Templated version of correlator booking which takes - // @parm N desired harmonic and @parm M number of particles. + /// @brief Templated version of correlator booking which takes + /// @parm N desired harmonic and @parm M number of particles. + /// @todo Rename to book(ECorrPtr, ...) template ECorrPtr bookECorrelator(const string name, const Scatter2DPtr hIn) { return bookECorrelator(name, Correlators::hVec(N, M), hIn); } - // @brief Templated version of gapped correlator booking which takes - // @parm N desired harmonic and @parm M number of particles. + /// @brief Templated version of gapped correlator booking which takes + /// @parm N desired harmonic and @parm M number of particles. + /// @todo Rename to book(ECorrPtr, ...) template ECorrPtr bookECorrelatorGap(const string name, const Scatter2DPtr hIn) { return bookECorrelatorGap(name, Correlators::hVec(N, M), hIn); } - - // @brief The stream method MUST be called in finalize() if one wants to stream + + // @brief The stream method MUST be called in finalize() if one wants to stream // correlators to the yoda file, in order to do reentrant finalize // (ie. multi-histogram merging) for the analysis. void stream() { for (auto ecItr = eCorrPtrs.begin(); ecItr != eCorrPtrs.end(); ++ecItr){ (*ecItr)->fillFromProfs(); corrPlot(list((*ecItr)->profBegin(), (*ecItr)->profEnd()), *ecItr); } } private: - - /// Bookkeeping of the event averaged correlators. + + /// Bookkeeping of the event averaged correlators. list eCorrPtrs; - + public: - + // @brief Constructor. Use CumulantAnalysis as base class for the // analysis to have access to functionality. CumulantAnalysis (string n) : Analysis(n), errorMethod(VARIANCE) {}; // @brief Helper method for turning correlators into Scatter2Ds. // Takes @parm h a pointer to the resulting Scatter2D, @parm binx // the x-bins and a function @parm func defining the transformation. // Makes no attempt at statistical errors. // See usage in the methods below. - // Can also be used directly in the analysis if a user wants to - // perform an unforseen transformation from correlators to + // Can also be used directly in the analysis if a user wants to + // perform an unforseen transformation from correlators to // Scatter2D. template static void fillScatter(Scatter2DPtr h, vector& binx, T func) { vector points; for (int i = 0, N = binx.size() - 1; i < N; ++i) { double xMid = (binx[i] + binx[i + 1]) / 2.0; double xeMin = fabs(xMid - binx[i]); double xeMax = fabs(xMid - binx[i + 1]); double yVal = func(i); if (std::isnan(yVal)) yVal = 0.; double yErr = 0; points.push_back(YODA::Point2D(xMid, yVal, xeMin, xeMax, yErr, yErr)); } h->reset(); h->points().clear(); for (int i = 0, N = points.size(); i < N; ++i) h->addPoint(points[i]); } // @brief Helper method for turning correlators into Scatter2Ds // with error estimates. // Takes @parm h a pointer to the resulting Scatter2D, @parm binx - // the x-bins, a function @parm func defining the transformation + // the x-bins, a function @parm func defining the transformation // and a vector of errors @parm err. // See usage in the methods below. - // Can also be used directly in the analysis if a user wants to - // perform an unforseen transformation from correlators to + // Can also be used directly in the analysis if a user wants to + // perform an unforseen transformation from correlators to // Scatter2D. template const void fillScatter(Scatter2DPtr h, vector& binx, F func, vector >& yErr) const { vector points; for (int i = 0, N = binx.size() - 1; i < N; ++i) { double xMid = (binx[i] + binx[i + 1]) / 2.0; double xeMin = fabs(xMid - binx[i]); double xeMax = fabs(xMid - binx[i + 1]); double yVal = func(i); if (std::isnan(yVal)) points.push_back(YODA::Point2D(xMid, 0., xeMin, xeMax,0., 0.)); else points.push_back(YODA::Point2D(xMid, yVal, xeMin, xeMax, yErr[i].first, yErr[i].second)); } h->reset(); h->points().clear(); - for (int i = 0, N = points.size(); i < N; ++i) + for (int i = 0, N = points.size(); i < N; ++i) h->addPoint(points[i]); } /// @brief Take the @parm n th power of all points in @parm hIn, - /// and put the result in @parm hOut. Optionally put a + /// and put the result in @parm hOut. Optionally put a /// @parm k constant below the root. static void nthPow(Scatter2DPtr hOut, const Scatter2DPtr hIn, const double& n, const double& k = 1.0) { if (n == 0 || n == 1) { cout << "Error: Do not take the 0th or 1st power of a Scatter2D," " use scale instead." << endl; return; } if (hIn->points().size() != hOut->points().size()) { cout << "nthRoot: Scatterplots: " << hIn->name() << " and " << hOut->name() << " not initialized with same length." << endl; return; } vector points; // The error pre-factor is k^(1/n) / n by Taylors formula. double eFac = pow(k,1./n)/n; for (auto b : hIn->points()) { double yVal = pow(k * b.y(),n); if (std::isnan(yVal)) points.push_back(YODA::Point2D(b.x(), 0., b.xErrMinus(), b.xErrPlus(), 0, 0 )); else { double yemin = abs(eFac * pow(yVal,1./(n - 1.))) * b.yErrMinus(); if (std::isnan(yemin)) yemin = b.yErrMinus(); double yemax = abs(eFac * pow(yVal,1./(n - 1.))) * b.yErrPlus(); if (std::isnan(yemax)) yemax = b.yErrPlus(); - points.push_back(YODA::Point2D(b.x(), yVal, b.xErrMinus(), + points.push_back(YODA::Point2D(b.x(), yVal, b.xErrMinus(), b.xErrPlus(), yemin, yemax )); - } + } } hOut->reset(); hOut->points().clear(); - for (int i = 0, N = points.size(); i < N; ++i) + for (int i = 0, N = points.size(); i < N; ++i) hOut->addPoint(points[i]); } - + /// @brief Take the @parm n th power of all points in @parm h, - /// and put the result back in the same Scatter2D. Optionally put a + /// and put the result back in the same Scatter2D. Optionally put a /// @parm k constant below the root. - static void nthPow(Scatter2DPtr h, const double& n, + static void nthPow(Scatter2DPtr h, const double& n, const double& k = 1.0) { if (n == 0 || n == 1) { cout << "Error: Do not take the 0th or 1st power of a Scatter2D," " use scale instead." << endl; return; } vector points; vector pIn = h->points(); // The error pre-factor is k^(1/n) / n by Taylors formula. double eFac = pow(k,1./n)/n; for (auto b : pIn) { double yVal = pow(k * b.y(),n); if (std::isnan(yVal)) - points.push_back(YODA::Point2D(b.x(), 0., b.xErrMinus(), + points.push_back(YODA::Point2D(b.x(), 0., b.xErrMinus(), b.xErrPlus(), 0, 0 )); else { double yemin = abs(eFac * pow(yVal,1./(n - 1.))) * b.yErrMinus(); if (std::isnan(yemin)) yemin = b.yErrMinus(); double yemax = abs(eFac * pow(yVal,1./(n - 1.))) * b.yErrPlus(); if (std::isnan(yemax)) yemax = b.yErrPlus(); - points.push_back(YODA::Point2D(b.x(), yVal, b.xErrMinus(), + points.push_back(YODA::Point2D(b.x(), yVal, b.xErrMinus(), b.xErrPlus(), yemin, yemax )); - } + } } h->reset(); h->points().clear(); for (int i = 0, N = points.size(); i < N; ++i) h->addPoint(points[i]); } - + // @brief Calculate the bootstrapped sample variance for the observable - // given by correlators and the transformation @parm func. + // given by correlators and the transformation @parm func. template static pair sampleVariance(T func) { // First we calculate the mean (two pass calculation). double avg = 0.; for (int i = 0; i < BOOT_BINS; ++i) avg += func(i); avg /= double(BOOT_BINS); // Then we find the variance. double var = 0.; for (int i = 0; i < BOOT_BINS; ++i) var += pow(func(i) - avg, 2.); var /= (double(BOOT_BINS) - 1); return pair(var, var); } // @brief Calculate the bootstrapped sample envelope for the observable - // given by correlators and the transformation @parm func. + // given by correlators and the transformation @parm func. template static pair sampleEnvelope(T func) { // First we calculate the mean. double avg = 0.; for (int i = 0; i < BOOT_BINS; ++i) avg += func(i); avg /= double(BOOT_BINS); double yMax = avg; double yMin = avg; // Then we find the envelope using the mean as initial value. for (int i = 0; i < BOOT_BINS; ++i) { double yVal = func(i); if (yMin > yVal) yMin = yVal; else if (yMax < yVal) yMax = yVal; } return pair(fabs(avg - yMin), fabs(yMax - avg)); } // @brief Selection method for which sample error to use, given // in the constructor. template const pair sampleError(T func) const { if (errorMethod == VARIANCE) return sampleVariance(func); else if (errorMethod == ENVELOPE) return sampleEnvelope(func); else cout << "Error: Error method not found!" << endl; return pair(0.,0.); - } + } // @brief Two particle integrated cn. const void cnTwoInt(Scatter2DPtr h, ECorrPtr e2) const { vector bins = e2->getBins(); vector binx = e2->getBinX(); // Assert bin size. if (binx.size() - 1 != bins.size()){ cout << "cnTwoInt: Bin size (x,y) differs!" << endl; return; } vector binPtrs; // The mean value of the cumulant. auto cn = [&] (int i) { return binPtrs[i]->mean(); }; // Error calculation. vector > yErr; for (int j = 0, N = bins.size(); j < N; ++j) { binPtrs = bins[j].getBinPtrs(); yErr.push_back(sampleError(cn)); } binPtrs = e2->getBinPtrs(); fillScatter(h, binx, cn, yErr); } - + // @brief Two particle integrated vn. const void vnTwoInt(Scatter2DPtr h, ECorrPtr e2) const { cnTwoInt(h, e2); nthPow(h, 0.5); } // @brief Put an event averaged correlator into a Scatter2D. // Reduces to cnTwoInt, but better with a proper name. const void corrPlot(Scatter2DPtr h, ECorrPtr e) const { cnTwoInt(h, e); } // @brief Put an event averaged correlator into Profile1Ds, one // for each bootstrapping bin. const void corrPlot(list profs, ECorrPtr e) const { vector corBins = e->getBins(); vector binx = e->getBinX(); auto ref = e->getReference(); auto refBins = ref.getBinPtrs(); // Assert bin size. if (binx.size() - 1 != corBins.size()){ cout << "corrPlot: Bin size (x,y) differs!" << endl; return; } list::iterator hItr = profs.begin(); // Loop over the boostrapped correlators. for (size_t i = 0; i < profs.size(); ++i, ++hItr) { vector profBins; // Numbers for the summary distribution double ne = 0., sow = 0., sow2 = 0.; for (size_t j = 0, N = binx.size() - 1; j < N; ++j) { - vector binPtrs = + vector binPtrs = corBins[j].getBinPtrs(); // Construct bin of the profiled quantities. We have no information - // (and no desire to add it) of sumWX of the profile, so really + // (and no desire to add it) of sumWX of the profile, so really // we should use a Dbn1D - but that does not work for Profile1D's. profBins.push_back( YODA::ProfileBin1D((*hItr)->bin(j).xEdges(), - YODA::Dbn2D( binPtrs[i]->numEntries(), binPtrs[i]->sumW(), + YODA::Dbn2D( binPtrs[i]->numEntries(), binPtrs[i]->sumW(), binPtrs[i]->sumW2(), 0., 0., binPtrs[i]->sumWX(), 0, 0))); ne += binPtrs[i]->numEntries(); sow += binPtrs[i]->sumW(); sow2 += binPtrs[i]->sumW2(); } // Reset the bins of the profiles. (*hItr)->reset(); (*hItr)->bins().clear(); // Add our new bins. // The total distribution (*hItr)->setTotalDbn(YODA::Dbn2D(ne,sow,sow2,0.,0.,0.,0.,0.)); // The bins. for (int j = 0, N = profBins.size(); j < N; ++j) (*hItr)->addBin(profBins[j]); // The reference flow in the underflow bin. (*hItr)->setUnderflow(YODA::Dbn2D(refBins[i]->numEntries(), refBins[i]->sumW(), refBins[i]->sumW2(), 0., 0., refBins[i]->sumWX(), 0., 0.)); } // End loop of bootstrapped correlators. } // @brief Four particle integrated cn. const void cnFourInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4) const { auto e2bins = e2->getBins(); auto e4bins = e4->getBins(); auto binx = e2->getBinX(); if (binx.size() - 1 != e2bins.size()){ cout << "cnFourInt: Bin size (x,y) differs!" << endl; return; } if (binx != e4->getBinX()) { cout << "Error in cnFourInt: Correlator x-binning differs!" << endl; return; } vector e2binPtrs; vector e4binPtrs; auto cn = [&] (int i) { double e22 = e2binPtrs[i]->mean() * e2binPtrs[i]->mean(); return e4binPtrs[i]->mean() - 2. * e22; }; // Error calculation. vector > yErr; for (int j = 0, N = e2bins.size(); j < N; ++j) { e2binPtrs = e2bins[j].getBinPtrs(); e4binPtrs = e4bins[j].getBinPtrs(); yErr.push_back(sampleError(cn)); } // Put the bin ptrs back in place. e2binPtrs = e2->getBinPtrs(); e4binPtrs = e4->getBinPtrs(); - fillScatter(h, binx, cn, yErr); + fillScatter(h, binx, cn, yErr); } // @brief Four particle integrated vn const void vnFourInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4) const { cnFourInt(h, e2, e4); nthPow(h, 0.25, -1.0); } - + // @brief Six particle integrated cn. - const void cnSixInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4, + const void cnSixInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4, ECorrPtr e6) const { auto e2bins = e2->getBins(); auto e4bins = e4->getBins(); auto e6bins = e6->getBins(); auto binx = e2->getBinX(); if (binx.size() - 1 != e2bins.size()){ cout << "cnSixInt: Bin size (x,y) differs!" << endl; return; } if (binx != e4->getBinX() || binx != e6->getBinX()) { cout << "Error in cnSixInt: Correlator x-binning differs!" << endl; return; } vector e2binPtrs; vector e4binPtrs; vector e6binPtrs; auto cn = [&] (int i) { - double e23 = pow(e2binPtrs[i]->mean(), 3.0); - return e6binPtrs[i]->mean() - 9.*e2binPtrs[i]->mean()*e4binPtrs[i]->mean() + + double e23 = pow(e2binPtrs[i]->mean(), 3.0); + return e6binPtrs[i]->mean() - 9.*e2binPtrs[i]->mean()*e4binPtrs[i]->mean() + 12.*e23; }; // Error calculation. vector > yErr; for (int j = 0, N = e2bins.size(); j < N; ++j) { e2binPtrs = e2bins[j].getBinPtrs(); e4binPtrs = e4bins[j].getBinPtrs(); e6binPtrs = e6bins[j].getBinPtrs(); yErr.push_back(sampleError(cn)); } // Put the bin ptrs back in place. e2binPtrs = e2->getBinPtrs(); e4binPtrs = e4->getBinPtrs(); e6binPtrs = e6->getBinPtrs(); fillScatter(h, binx, cn, yErr); } - + // @brief Six particle integrated vn const void vnSixInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4, ECorrPtr e6) const { cnSixInt(h, e2, e4, e6); nthPow(h, 1./6., 0.25); } - + // @brief Eight particle integrated cn. const void cnEightInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4, ECorrPtr e6, ECorrPtr e8) const { auto e2bins = e2->getBins(); auto e4bins = e4->getBins(); auto e6bins = e6->getBins(); auto e8bins = e8->getBins(); auto binx = e2->getBinX(); if (binx.size() - 1 != e2bins.size()){ cout << "cnEightInt: Bin size (x,y) differs!" << endl; return; } - if (binx != e4->getBinX() || binx != e6->getBinX() || + if (binx != e4->getBinX() || binx != e6->getBinX() || binx != e8->getBinX()) { cout << "Error in cnEightInt: Correlator x-binning differs!" << endl; return; } vector e2binPtrs; vector e4binPtrs; vector e6binPtrs; vector e8binPtrs; auto cn = [&] (int i ) { double e22 = e2binPtrs[i]->mean() * e2binPtrs[i]->mean(); double e24 = e22 * e22; double e42 = e4binPtrs[i]->mean() * e4binPtrs[i]->mean(); - return e8binPtrs[i]->mean() - 16. * e6binPtrs[i]->mean() * + return e8binPtrs[i]->mean() - 16. * e6binPtrs[i]->mean() * e2binPtrs[i]->mean() - 18. * e42 + 144. * e4binPtrs[i]->mean()*e22 - - 144. * e24; + - 144. * e24; }; // Error calculation. vector > yErr; for (int j = 0, N = e2bins.size(); j < N; ++j) { e2binPtrs = e2bins[j].getBinPtrs(); e4binPtrs = e4bins[j].getBinPtrs(); e6binPtrs = e6bins[j].getBinPtrs(); e8binPtrs = e8bins[j].getBinPtrs(); yErr.push_back(sampleError(cn)); } // Put the bin ptrs back in place. e2binPtrs = e2->getBinPtrs(); e4binPtrs = e4->getBinPtrs(); e6binPtrs = e6->getBinPtrs(); e8binPtrs = e8->getBinPtrs(); fillScatter(h, binx, cn, yErr); } // @brief Eight particle integrated vn const void vnEightInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4, ECorrPtr e6, ECorrPtr e8) const { cnEightInt(h, e2, e4, e6, e8); nthPow(h, 1./8., -1./33.); } // @brief Two particle differential vn. const void vnTwoDiff(Scatter2DPtr h, ECorrPtr e2Dif) const { auto e2bins = e2Dif->getBins(); auto ref = e2Dif->getReference(); auto binx = e2Dif->getBinX(); if (binx.size() -1 != e2bins.size()) { cout << "vnTwoDif: Bin size (x,y) differs!" << endl; return; } vector e2binPtrs; vector refPtrs; auto vn = [&] (int i) { // Test reference flow. if (ref.mean() <= 0) return 0.; return e2binPtrs[i]->mean() / sqrt(ref.mean()); }; // We need here a seperate error function, as we don't // iterate over the reference flow. auto vnerr = [&] (int i) { // Test reference flow. if (refPtrs[i]->mean() <=0) return 0.; return e2binPtrs[i]->mean() / sqrt(refPtrs[i]->mean()); }; // Error calculation. vector > yErr; refPtrs = ref.getBinPtrs(); for (int j = 0, N = e2bins.size(); j < N; ++j) { e2binPtrs = e2bins[j].getBinPtrs(); yErr.push_back(sampleError(vnerr)); } // Put the e2binPtrs back in place. e2binPtrs = e2Dif->getBinPtrs(); fillScatter(h, binx, vn); } // @brief Four particle differential vn. - const void vnFourDiff(Scatter2DPtr h, ECorrPtr e2Dif, + const void vnFourDiff(Scatter2DPtr h, ECorrPtr e2Dif, ECorrPtr e4Dif) const { auto e2bins = e2Dif->getBins(); auto e4bins = e4Dif->getBins(); auto ref2 = e2Dif->getReference(); auto ref4 = e4Dif->getReference(); auto binx = e2Dif->getBinX(); if (binx.size() - 1 != e2bins.size()){ cout << "vnFourDif: Bin size (x,y) differs!" << endl; return; } if (binx != e4Dif->getBinX()) { cout << "Error in vnFourDif: Correlator x-binning differs!" << endl; return; } vector e2binPtrs; vector e4binPtrs; vector ref2Ptrs; vector ref4Ptrs; double denom = 2 * ref2.mean() * ref2.mean() - ref4.mean(); auto vn = [&] (int i) { // Test denominator. if (denom <= 0 ) return 0.; return ((2 * ref2.mean() * e2bins[i].mean() - e4bins[i].mean()) / pow(denom, 0.75)); }; // We need here a seperate error function, as we don't // iterate over the reference flow. auto vnerr = [&] (int i) { double denom2 = 2 * ref2Ptrs[i]->mean() * ref2Ptrs[i]->mean() - ref4Ptrs[i]->mean(); // Test denominator. if (denom2 <= 0) return 0.; - return ((2 * ref2Ptrs[i]->mean() * e2binPtrs[i]->mean() - + return ((2 * ref2Ptrs[i]->mean() * e2binPtrs[i]->mean() - e4binPtrs[i]->mean()) / pow(denom2, 0.75)); }; // Error calculation. vector > yErr; ref2Ptrs = ref2.getBinPtrs(); ref4Ptrs = ref4.getBinPtrs(); for (int j = 0, N = e2bins.size(); j < N; ++j) { e2binPtrs = e2bins[j].getBinPtrs(); e4binPtrs = e4bins[j].getBinPtrs(); yErr.push_back(sampleError(vnerr)); } // Put the binPtrs back in place. e2binPtrs = e2Dif->getBinPtrs(); e4binPtrs = e4Dif->getBinPtrs(); fillScatter(h, binx, vn, yErr); } }; // End class CumulantAnalysis. } // End Rivet namespace. #endif diff --git a/include/Rivet/Tools/RivetYODA.hh b/include/Rivet/Tools/RivetYODA.hh --- a/include/Rivet/Tools/RivetYODA.hh +++ b/include/Rivet/Tools/RivetYODA.hh @@ -1,578 +1,579 @@ #ifndef RIVET_RIVETYODA_HH #define RIVET_RIVETYODA_HH #include "Rivet/Config/RivetCommon.hh" #include "YODA/AnalysisObject.h" #include "YODA/Counter.h" #include "YODA/Histo1D.h" #include "YODA/Histo2D.h" #include "YODA/Profile1D.h" #include "YODA/Profile2D.h" #include "YODA/Scatter1D.h" #include "YODA/Scatter2D.h" #include "YODA/Scatter3D.h" #include #include namespace YODA { + typedef std::shared_ptr AnalysisObjectPtr; typedef std::shared_ptr CounterPtr; typedef std::shared_ptr Histo1DPtr; typedef std::shared_ptr Histo2DPtr; typedef std::shared_ptr Profile1DPtr; typedef std::shared_ptr Profile2DPtr; typedef std::shared_ptr Scatter1DPtr; typedef std::shared_ptr Scatter2DPtr; typedef std::shared_ptr Scatter3DPtr; } namespace Rivet { class AnalysisObjectWrapper { public: virtual ~AnalysisObjectWrapper() {} virtual YODA::AnalysisObject* operator->() = 0; virtual YODA::AnalysisObject* operator->() const = 0; virtual const YODA::AnalysisObject & operator*() const = 0; /// @todo Rename to setActive(idx) virtual void setActiveWeightIdx(unsigned int iWeight) = 0; bool operator ==(const AnalysisObjectWrapper& p) { return (this == &p); } bool operator !=(const AnalysisObjectWrapper& p) { return (this != &p); } protected: /// @todo do we need this? // virtual void reset() = 0; }; /// @todo /// implement scatter1dptr and scatter2dptr here /// these need to be multi-weighted eventually. /* class Scatter1DPtr : public AnalysisObjectPtr { public: Scatter1DPtr() : _persistent() { } Scatter1DPtr(size_t len_of_weightvec, const YODA::Scatter1D& p) { for (size_t m = 0; m < len_of_weightvec; ++m) _persistent.push_back(make_shared(p)); } bool operator!() const { return !_persistent; } explicit operator bool() const { return bool(_persistent); } YODA::Scatter1D* operator->() { return _persistent.get(); } YODA::Scatter1D* operator->() const { return _persistent.get(); } YODA::Scatter1D & operator*() { return *_persistent; } const YODA::Scatter1D & operator*() const { return *_persistent; } protected: vector _persistent; }; class Scatter2DPtr : public AnalysisObjectPtr { public: Scatter2DPtr(size_t len_of_weightvec, const YODA::Scatter2D& p) { for (size_t m = 0; m < len_of_weightvec; ++m) _persistent.push_back(make_shared(p)); } Scatter2DPtr() : _persistent() { } bool operator!() { return !_persistent; } explicit operator bool() { return bool(_persistent); } YODA::Scatter2D* operator->() { return _persistent.get(); } YODA::Scatter2D* operator->() const { return _persistent.get(); } YODA::Scatter2D & operator*() { return *_persistent; } const YODA::Scatter2D & operator*() const { return *_persistent; } protected: vector _persistent; }; class Scatter3DPtr : public AnalysisObjectPtr { public: Scatter3DPtr(size_t len_of_weightvec, const YODA::Scatter3D& p) { for (size_t m = 0; m < len_of_weightvec; ++m) _persistent.push_back(make_shared(p)); } Scatter3DPtr() : _persistent() { } bool operator!() { return !_persistent; } explicit operator bool() { return bool(_persistent); } YODA::Scatter3D* operator->() { return _persistent.get(); } YODA::Scatter3D* operator->() const { return _persistent.get(); } YODA::Scatter3D & operator*() { return *_persistent; } const YODA::Scatter3D & operator*() const { return *_persistent; } protected: vector _persistent; }; */ class MultiweightAOWrapper : public AnalysisObjectWrapper { public: using Inner = YODA::AnalysisObject; virtual void newSubEvent() = 0; virtual void pushToPersistent(const vector >& weight) = 0; virtual YODA::AnalysisObjectPtr activeYODAPtr() const = 0; }; using Weight = double; template using Fill = pair; template using Fills = multiset>; // TODO TODO TODO // need to override the old fill method too! // otherwise we can't intercept existing fill calls in analysis code // TODO TODO TODO /// Wrappers for analysis objects to store all fills unaggregated, until collapsed template class TupleWrapper; template<> class TupleWrapper : public YODA::Counter { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Counter & h) : YODA::Counter(h) {} // todo: do we need to deal with users using fractions directly? void fill( double weight=1.0, double fraction=1.0 ) { fills_.insert( {YODA::Counter::FillType(),weight} ); } void reset() { fills_.clear(); } const Fills & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Histo1D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Histo1D & h) : YODA::Histo1D(h) {} // todo: do we need to deal with users using fractions directly? void fill( double x, double weight=1.0, double fraction=1.0 ) { if ( std::isnan(x) ) throw YODA::RangeError("X is NaN"); fills_.insert( { x , weight } ); } void reset() { fills_.clear(); } const Fills & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Profile1D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Profile1D & h) : YODA::Profile1D(h) {} // todo: do we need to deal with users using fractions directly? void fill( double x, double y, double weight=1.0, double fraction=1.0 ) { if ( std::isnan(x) ) throw YODA::RangeError("X is NaN"); if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN"); fills_.insert( { YODA::Profile1D::FillType{x,y}, weight } ); } void reset() { fills_.clear(); } const Fills & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Histo2D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Histo2D & h) : YODA::Histo2D(h) {} // todo: do we need to deal with users using fractions directly? void fill( double x, double y, double weight=1.0, double fraction=1.0 ) { if ( std::isnan(x) ) throw YODA::RangeError("X is NaN"); if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN"); fills_.insert( { YODA::Histo2D::FillType{x,y}, weight } ); } void reset() { fills_.clear(); } const Fills & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Profile2D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Profile2D & h) : YODA::Profile2D(h) {} // todo: do we need to deal with users using fractions directly? void fill( double x, double y, double z, double weight=1.0, double fraction=1.0 ) { if ( std::isnan(x) ) throw YODA::RangeError("X is NaN"); if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN"); if ( std::isnan(z) ) throw YODA::RangeError("Z is NaN"); fills_.insert( { YODA::Profile2D::FillType{x,y,z}, weight } ); } void reset() { fills_.clear(); } const Fills & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Scatter1D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Scatter1D & h) : YODA::Scatter1D(h) {} }; template<> class TupleWrapper : public YODA::Scatter2D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Scatter2D & h) : YODA::Scatter2D(h) {} }; template<> class TupleWrapper : public YODA::Scatter3D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Scatter3D & h) : YODA::Scatter3D(h) {} }; template class Wrapper : public MultiweightAOWrapper { friend class Analysis; public: using Inner = T; /* @todo * some things are not really well-defined here * for instance: fill() in the finalize() method and integral() in * the analyze() method. */ Wrapper() = default; Wrapper(const vector& weightnames, const T & p); ~Wrapper(); typename T::Ptr active() const; /* @todo this probably need to loop over all? */ bool operator!() const { return !_active; } // Don't use active() here, assert will catch explicit operator bool() const { return static_cast(_active); } // Don't use active() here, assert will catch T * operator->() { return active().get(); } T * operator->() const { return active().get(); } T & operator*() { return *active(); } const T & operator*() const { return *active(); } /* @todo * these need to be re-thought out. void reset() { active()->reset(); } */ /* @todo * these probably need to loop over all? * do we even want to provide equality? */ /* @todo * how about no. friend bool operator==(Wrapper a, Wrapper b){ if (a._persistent.size() != b._persistent.size()) return false; for (size_t i = 0; i < a._persistent.size(); i++) { if (a._persistent.at(i) != b._persistent.at(i)) { return false; } } return true; } friend bool operator!=(Wrapper a, Wrapper b){ return !(a == b); } friend bool operator<(Wrapper a, Wrapper b){ if (a._persistent.size() >= b._persistent.size()) return false; for (size_t i = 0; i < a._persistent.size(); i++) { if (*(a._persistent.at(i)) >= *(b._persistent.at(i))) { return false; } } return true; } */ private: void setActiveWeightIdx(unsigned int iWeight) { _active = _persistent.at(iWeight); } /* this is for dev only---we shouldn't need this in real runs. */ void unsetActiveWeight() { _active.reset(); } void newSubEvent(); virtual YODA::AnalysisObjectPtr activeYODAPtr() const { return _active; } const vector & persistent() const { return _persistent; } /* to be implemented for each type */ void pushToPersistent(const vector >& weight); /* M of these, one for each weight */ vector _persistent; /* N of these, one for each event in evgroup */ vector::Ptr> _evgroup; typename T::Ptr _active; // do we need implicit cast? // operator typename T::Ptr () { // return _active; // } friend class AnalysisHandler; }; /// We need our own shared_ptr class, so we can dispatch -> and * /// all the way down to the inner YODA analysis objects /// /// TODO: provide remaining functionality that shared_ptr has (not needed right now) /// template class rivet_shared_ptr { public: typedef T value_type; rivet_shared_ptr() = default; rivet_shared_ptr(decltype(nullptr)) : _p(nullptr) {} /// Convenience constructor, pass through to the Wrapper constructor rivet_shared_ptr(const vector& weightNames, const typename T::Inner & p) : _p( make_shared(weightNames, p) ) {} template rivet_shared_ptr(const shared_ptr & p) : _p(p) {} template rivet_shared_ptr(const rivet_shared_ptr & p) : _p(p.get()) {} // Goes right through to the active YODA object's members T & operator->() { return *_p; } const T & operator->() const { return *_p; } // The active YODA object typename T::Inner & operator*() { return **_p; } const typename T::Inner & operator*() const { return **_p; } bool operator!() const { return !_p || !(*_p); } explicit operator bool() const { return _p && bool(*_p); } template bool operator==(const rivet_shared_ptr & other) const { return _p == other._p; } template bool operator!=(const rivet_shared_ptr & other) const { return _p != other._p; } template bool operator<(const rivet_shared_ptr & other) const { return _p < other._p; } template bool operator>(const rivet_shared_ptr & other) const { return _p > other._p; } template bool operator<=(const rivet_shared_ptr & other) const { return _p <= other._p; } template bool operator>=(const rivet_shared_ptr & other) const { return _p >= other._p; } shared_ptr get() const { return _p; } private: shared_ptr _p; }; // every object listed here needs a virtual fill method in YODA, // otherwise the Tuple fakery won't work. using MultiweightAOPtr = rivet_shared_ptr; using Histo1DPtr = rivet_shared_ptr>; using Histo2DPtr = rivet_shared_ptr>; using Profile1DPtr = rivet_shared_ptr>; using Profile2DPtr = rivet_shared_ptr>; using CounterPtr = rivet_shared_ptr>; using Scatter1DPtr = rivet_shared_ptr>; using Scatter2DPtr = rivet_shared_ptr>; using Scatter3DPtr = rivet_shared_ptr>; using YODA::Counter; using YODA::Histo1D; using YODA::HistoBin1D; using YODA::Histo2D; using YODA::HistoBin2D; using YODA::Profile1D; using YODA::ProfileBin1D; using YODA::Profile2D; using YODA::ProfileBin2D; using YODA::Scatter1D; using YODA::Point1D; using YODA::Scatter2D; using YODA::Point2D; using YODA::Scatter3D; using YODA::Point3D; /// Function to get a map of all the refdata in a paper with the /// given @a papername. map getRefData(const string& papername); /// @todo Also provide a Scatter3D getRefData() version? /// Get the file system path to the reference file for this paper. string getDatafilePath(const string& papername); /// Traits class to access the type of the AnalysisObject in the /// reference files. template struct ReferenceTraits {}; template<> struct ReferenceTraits { typedef Counter RefT; }; template<> struct ReferenceTraits { typedef Scatter1D RefT; }; template<> struct ReferenceTraits { typedef Scatter2D RefT; }; template<> struct ReferenceTraits { typedef Scatter2D RefT; }; template<> struct ReferenceTraits { typedef Scatter2D RefT; }; template<> struct ReferenceTraits { typedef Scatter3D RefT; }; template<> struct ReferenceTraits { typedef Scatter3D RefT; }; template<> struct ReferenceTraits { typedef Scatter3D RefT; }; /// If @a dst and @a src both are of same subclass T, copy the /// contents of @a src into @a dst and return true. Otherwise return /// false. template inline bool aocopy(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst) { shared_ptr tsrc = dynamic_pointer_cast(src); if ( !tsrc ) return false; shared_ptr tdst = dynamic_pointer_cast(dst); if ( !tdst ) return false; *tdst = *tsrc; return true; } /// If @a dst and @a src both are of same subclass T, add the /// contents of @a src into @a dst and return true. Otherwise return /// false. template inline bool aoadd(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale) { shared_ptr tsrc = dynamic_pointer_cast(src); if ( !tsrc ) return false; shared_ptr tdst = dynamic_pointer_cast(dst); if ( !tdst ) return false; tsrc->scaleW(scale); *tdst += *tsrc; return true; } /// If @a dst is the same subclass as @a src, copy the contents of @a /// src into @a dst and return true. Otherwise return false. bool copyao(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst); /// If @a dst is the same subclass as @a src, scale the contents of /// @a src with @a scale and add it to @a dst and return true. Otherwise /// return false. bool addaos(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale); /// Check if two analysis objects have the same binning or, if not /// binned, are in other ways compatible. template inline bool bookingCompatible(TPtr a, TPtr b) { return a->sameBinning(*b); } inline bool bookingCompatible(CounterPtr, CounterPtr) { return true; } inline bool bookingCompatible(Scatter1DPtr a, Scatter1DPtr b) { return a->numPoints() == b->numPoints(); } inline bool bookingCompatible(Scatter2DPtr a, Scatter2DPtr b) { return a->numPoints() == b->numPoints(); } inline bool bookingCompatible(Scatter3DPtr a, Scatter3DPtr b) { return a->numPoints() == b->numPoints(); } } #endif diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc --- a/src/Core/AnalysisHandler.cc +++ b/src/Core/AnalysisHandler.cc @@ -1,755 +1,757 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/AnalysisHandler.hh" #include "Rivet/Analysis.hh" #include "Rivet/Tools/ParticleName.hh" #include "Rivet/Tools/BeamConstraint.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/Beam.hh" #include "YODA/IO.h" #include #include using std::cout; using std::cerr; namespace { inline std::vector split(const std::string& input, const std::string& regex) { // passing -1 as the submatch index parameter performs splitting std::regex re(regex); std::sregex_token_iterator first{input.begin(), input.end(), re, -1}, last; return {first, last}; } } namespace Rivet { AnalysisHandler::AnalysisHandler(const string& runname) : _runname(runname), _eventCounter(vector(), Counter()), _xs(vector(), Scatter1D()), _initialised(false), _ignoreBeams(false), _defaultWeightIdx(0) {} AnalysisHandler::~AnalysisHandler() { } Log& AnalysisHandler::getLog() const { return Log::getLog("Rivet.Analysis.Handler"); } /// http://stackoverflow.com/questions/4654636/how-to-determine-if-a-string-is-a-number-with-c - bool is_number(const std::string& s) { - std::string::const_iterator it = s.begin(); - while (it != s.end() && std::isdigit(*it)) ++it; - return !s.empty() && it == s.end(); + namespace { + bool is_number(const std::string& s) { + std::string::const_iterator it = s.begin(); + while (it != s.end() && std::isdigit(*it)) ++it; + return !s.empty() && it == s.end(); + } } /// Check if any of the weightnames is not a number - bool AnalysisHandler::haveNamedWeights() { + bool AnalysisHandler::haveNamedWeights() const { bool dec=false; for (unsigned int i=0;i<_weightNames.size();++i) { string s = _weightNames[i]; if (!is_number(s)) { dec=true; break; } } return dec; } void AnalysisHandler::init(const GenEvent& ge) { if (_initialised) throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!"); /// @todo Should the Rivet analysis objects know about weight names? setRunBeams(Rivet::beams(ge)); MSG_DEBUG("Initialising the analysis handler"); _eventNumber = ge.event_number(); setWeightNames(ge); if (haveNamedWeights()) MSG_INFO("Using named weights"); else MSG_INFO("NOT using named weights. Using first weight as nominal weight"); _eventCounter = CounterPtr(weightNames(), Counter("_EVTCOUNT")); // Set the cross section based on what is reported by this event. if (ge.cross_section()) { MSG_TRACE("Getting cross section."); double xs = ge.cross_section()->cross_section(); double xserr = ge.cross_section()->cross_section_error(); setCrossSection(xs, xserr); } // Check that analyses are beam-compatible, and remove those that aren't const size_t num_anas_requested = analysisNames().size(); vector anamestodelete; for (const AnaHandle a : _analyses) { if (!_ignoreBeams && !a->isCompatible(beams())) { //MSG_DEBUG(a->name() << " requires beams " << a->requiredBeams() << " @ " << a->requiredEnergies() << " GeV"); anamestodelete.push_back(a->name()); } } for (const string& aname : anamestodelete) { MSG_WARNING("Analysis '" << aname << "' is incompatible with the provided beams: removing"); removeAnalysis(aname); } if (num_anas_requested > 0 && analysisNames().empty()) { cerr << "All analyses were incompatible with the first event's beams\n" << "Exiting, since this probably wasn't intentional!" << endl; exit(1); } // Warn if any analysis' status is not unblemished for (const AnaHandle a : analyses()) { if (toUpper(a->status()) == "PRELIMINARY") { MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!"); } else if (toUpper(a->status()) == "OBSOLETE") { MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!"); } else if (toUpper(a->status()).find("UNVALIDATED") != string::npos) { MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!"); } } // Initialize the remaining analyses _stage = Stage::INIT; for (AnaHandle a : _analyses) { MSG_DEBUG("Initialising analysis: " << a->name()); try { // Allow projection registration in the init phase onwards a->_allowProjReg = true; a->init(); //MSG_DEBUG("Checking consistency of analysis: " << a->name()); //a->checkConsistency(); } catch (const Error& err) { cerr << "Error in " << a->name() << "::init method: " << err.what() << endl; exit(1); } MSG_DEBUG("Done initialising analysis: " << a->name()); } _stage = Stage::OTHER; _initialised = true; MSG_DEBUG("Analysis handler initialised"); } void AnalysisHandler::setWeightNames(const GenEvent& ge) { /// reroute the print output to a std::stringstream and process /// The iteration is done over a map in hepmc2 so this is safe std::ostringstream stream; ge.weights().print(stream); // Super lame, I know string str = stream.str(); std::regex re("(([^()]+))"); // Regex for stuff enclosed by parentheses () size_t idx = 0; for (std::sregex_iterator i = std::sregex_iterator(str.begin(), str.end(), re); i != std::sregex_iterator(); ++i ) { std::smatch m = *i; vector temp = ::split(m.str(), "[,]"); if (temp.size() ==2) { MSG_DEBUG("Name of weight #" << _weightNames.size() << ": " << temp[0]); // store the default weight based on weight names if (temp[0] == "Weight" || temp[0] == "0" || temp[0] == "Default") { MSG_DEBUG(_weightNames.size() << " is being used as the nominal."); _weightNames.push_back(""); _defaultWeightIdx = idx; } else _weightNames.push_back(temp[0]); idx++; } } } void AnalysisHandler::analyze(const GenEvent& ge) { // Call init with event as template if not already initialised if (!_initialised) init(ge); assert(_initialised); // Ensure that beam details match those from the first event (if we're checking beams) if ( !_ignoreBeams ) { const PdgIdPair beams = Rivet::beamIds(ge); const double sqrts = Rivet::sqrtS(ge); if (!compatible(beams, _beams) || !fuzzyEquals(sqrts, sqrtS())) { cerr << "Event beams mismatch: " << PID::toBeamsString(beams) << " @ " << sqrts/GeV << " GeV" << " vs. first beams " << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl; exit(1); } } // Create the Rivet event wrapper /// @todo Filter/normalize the event here Event event(ge); // Set the cross section based on what is reported by this event. MSG_TRACE("Getting cross-section."); if (ge.cross_section()) { MSG_TRACE("Getting cross-section from GenEvent."); double xs = ge.cross_section()->cross_section(); double xserr = ge.cross_section()->cross_section_error(); setCrossSection(xs, xserr); } // Won't happen for first event because _eventNumber is set in init() if (_eventNumber != ge.event_number()) { /// @todo Can we get away with not passing a matrix? MSG_TRACE("AnalysisHandler::analyze(): Pushing _eventCounter to persistent."); _eventCounter.get()->pushToPersistent(_subEventWeights); // if this is indeed a new event, push the temporary // histograms and reset for (const AnaHandle& a : _analyses) { for (auto ao : a->analysisObjects()) { MSG_TRACE("AnalysisHandler::analyze(): Pushing " << a->name() << "'s " << ao->name() << " to persistent."); ao.get()->pushToPersistent(_subEventWeights); } MSG_TRACE("AnalysisHandler::analyze(): finished pushing " << a->name() << "'s objects to persistent."); } _eventNumber = ge.event_number(); MSG_DEBUG("nominal event # " << _eventCounter.get()->_persistent[0]->numEntries()); MSG_DEBUG("nominal sum of weights: " << _eventCounter.get()->_persistent[0]->sumW()); MSG_DEBUG("Event has " << _subEventWeights.size() << " sub events."); _subEventWeights.clear(); } MSG_TRACE("starting new sub event"); _eventCounter.get()->newSubEvent(); for (const AnaHandle& a : _analyses) { for (auto ao : a->analysisObjects()) { ao.get()->newSubEvent(); } } _subEventWeights.push_back(event.weights()); MSG_DEBUG("Analyzing subevent #" << _subEventWeights.size() - 1 << "."); _eventCounter->fill(); // Run the analyses for (AnaHandle a : _analyses) { MSG_TRACE("About to run analysis " << a->name()); try { a->analyze(event); } catch (const Error& err) { cerr << "Error in " << a->name() << "::analyze method: " << err.what() << endl; exit(1); } MSG_TRACE("Finished running analysis " << a->name()); } if ( _dumpPeriod > 0 && numEvents()%_dumpPeriod == 0 ) { MSG_INFO("Dumping intermediate results to " << _dumpFile << "."); _dumping = true; finalize(); _dumping = false; writeData(_dumpFile); } } void AnalysisHandler::analyze(const GenEvent* ge) { if (ge == nullptr) { MSG_ERROR("AnalysisHandler received null pointer to GenEvent"); //throw Error("AnalysisHandler received null pointer to GenEvent"); } analyze(*ge); } void AnalysisHandler::finalize() { if (!_initialised) return; MSG_INFO("Finalising analyses"); // First push all analyses' objects to persistent MSG_TRACE("AnalysisHandler::finalize(): Pushing analysis objects to persistent."); _eventCounter.get()->pushToPersistent(_subEventWeights); for (const AnaHandle& a : _analyses) { for (auto ao : a->analysisObjects()) ao.get()->pushToPersistent(_subEventWeights); } // First we make copies of all analysis objects. map backupAOs; for (auto ao : getYodaAOs(false, true, false) ) backupAOs[ao->path()] = YODA::AnalysisObjectPtr(ao->newclone()); // Finalize all the histograms for (const AnaHandle& a : _analyses) { // a->setCrossSection(_xs); for (size_t iW = 0; iW < numWeights(); iW++) { _eventCounter.get()->setActiveWeightIdx(iW); _xs.get()->setActiveWeightIdx(iW); for (auto ao : a->analysisObjects()) ao.get()->setActiveWeightIdx(iW); MSG_TRACE("Running " << a->name() << "::finalize() for weight " << iW << "."); try { if ( !_dumping || a->info().reentrant() ) a->finalize(); else if ( _dumping == 1 && iW == 0 ) MSG_INFO("Skipping periodic dump of " << a->name() << " as it is not declared reentrant."); } catch (const Error& err) { cerr << "Error in " << a->name() << "::finalize method: " << err.what() << endl; exit(1); } } } // Now we copy all analysis objects to the list of finalized // ones, and restore the value to their original ones. _finalizedAOs.clear(); for ( auto ao : getYodaAOs(false, false, false) ) _finalizedAOs.push_back(YODA::AnalysisObjectPtr(ao->newclone())); for ( auto ao : getYodaAOs(false, true, false) ) { // TODO: This should be possible to do in a nicer way, with a flag etc. if (ao->path().find("/FINAL") != std::string::npos) continue; auto aoit = backupAOs.find(ao->path()); if ( aoit == backupAOs.end() ) { AnaHandle ana = analysis(split(ao->path(), "/")[0]); if ( ana ) ana->removeAnalysisObject(ao->path()); } else copyao(aoit->second, ao); } // Print out number of events processed const int nevts = numEvents(); MSG_INFO("Processed " << nevts << " event" << (nevts != 1 ? "s" : "")); // // Delete analyses // MSG_DEBUG("Deleting analyses"); // _analyses.clear(); // Print out MCnet boilerplate cout << '\n'; cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << '\n'; cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:1003.0694 (http://arxiv.org/abs/1003.0694)" << endl; } AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname, std::map pars) { // Make an option handle. std::string parHandle = ""; for (map::iterator par = pars.begin(); par != pars.end(); ++par) { parHandle +=":"; parHandle += par->first + "=" + par->second; } return addAnalysis(analysisname + parHandle); } AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname) { // Check for a duplicate analysis /// @todo Might we want to be able to run an analysis twice, with different params? /// Requires avoiding histo tree clashes, i.e. storing the histos on the analysis objects. string ananame = analysisname; vector anaopt = split(analysisname, ":"); if ( anaopt.size() > 1 ) ananame = anaopt[0]; AnaHandle analysis( AnalysisLoader::getAnalysis(ananame) ); if (analysis.get() != 0) { // < Check for null analysis. MSG_DEBUG("Adding analysis '" << analysisname << "'"); map opts; for ( int i = 1, N = anaopt.size(); i < N; ++i ) { vector opt = split(anaopt[i], "="); if ( opt.size() != 2 ) { MSG_WARNING("Error in option specification. Skipping analysis " << analysisname); return *this; } if ( !analysis->info().validOption(opt[0], opt[1]) ) { MSG_WARNING("Cannot set option '" << opt[0] << "' to '" << opt[1] << "'. Skipping analysis " << analysisname); return *this; } opts[opt[0]] = opt[1]; } for ( auto opt: opts) { analysis->_options[opt.first] = opt.second; analysis->_optstring += ":" + opt.first + "=" + opt.second; } for (const AnaHandle& a : _analyses) { if (a->name() == analysis->name() ) { MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate"); return *this; } } analysis->_analysishandler = this; _analyses.insert(analysis); } else { MSG_WARNING("Analysis '" << analysisname << "' not found."); } // MSG_WARNING(_analyses.size()); // for (const AnaHandle& a : _analyses) MSG_WARNING(a->name()); return *this; } AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) { std::shared_ptr toremove; for (const AnaHandle a : _analyses) { if (a->name() == analysisname) { toremove = a; break; } } if (toremove.get() != 0) { MSG_DEBUG("Removing analysis '" << analysisname << "'"); _analyses.erase(toremove); } return *this; } void AnalysisHandler::addData(const std::vector& aos) { for (const YODA::AnalysisObjectPtr ao : aos) { string path = ao->path(); if ( path.substr(0, 5) != "/RAW/" ) { _orphanedPreloads.push_back(ao); continue; } path = path.substr(4); ao->setPath(path); if (path.size() > 1) { // path > "/" try { const string ananame = ::split(path, "/")[0]; AnaHandle a = analysis(ananame); /// @todo FIXXXXX //MultiweightAOPtr mao = ????; /// @todo generate right Multiweight object from ao //a->addAnalysisObject(mao); /// @todo Need to statistically merge... } catch (const Error& e) { MSG_TRACE("Adding analysis object " << path << " to the list of orphans."); _orphanedPreloads.push_back(ao); } } } } void AnalysisHandler::stripOptions(YODA::AnalysisObjectPtr ao, const vector & delopts) const { string path = ao->path(); string ananame = split(path, "/")[0]; vector anaopts = split(ananame, ":"); for ( int i = 1, N = anaopts.size(); i < N; ++i ) for ( auto opt : delopts ) if ( opt == "*" || anaopts[i].find(opt + "=") == 0 ) path.replace(path.find(":" + anaopts[i]), (":" + anaopts[i]).length(), ""); ao->setPath(path); } - void AnalysisHandler:: - mergeYodas(const vector & aofiles, const vector & delopts, bool equiv) { - vector< vector > aosv; - vector xsecs; - vector xsecerrs; - vector sows; - set ananames; - _eventCounter->reset(); + // void AnalysisHandler::mergeYodas(const vector & aofiles, + // const vector & delopts, bool equiv) { + // vector< vector > aosv; + // vector xsecs; + // vector xsecerrs; + // vector sows; + // set ananames; + // _eventCounter->reset(); - // First scan all files and extract analysis objects and add the - // corresponding analyses.. - for (const string& file : aofiles ) { - Scatter1DPtr xsec; - CounterPtr sow; + // // First scan all files and extract analysis objects and add the + // // corresponding analyses.. + // for (const string& file : aofiles ) { + // Scatter1DPtr xsec; + // CounterPtr sow; - // For each file make sure that cross section and sum-of-weights - // objects are present and stor all RAW ones in a vector; - vector aos; - try { - /// @todo Use new YODA SFINAE to fill the smart ptr vector directly - vector aos_raw; - YODA::read(file, aos_raw); - for (YODA::AnalysisObject* aor : aos_raw) { - YODA::AnalysisObjectPtr ao(aor); - if ( ao->path().substr(0, 5) != "/RAW/" ) continue; - ao->setPath(ao->path().substr(4)); - if ( ao->path() == "/_XSEC" ) - xsec = dynamic_pointer_cast(ao); - else if ( ao->path() == "/_EVTCOUNT" ) - sow = dynamic_pointer_cast(ao); - else { - stripOptions(ao, delopts); - string ananame = split(ao->path(), "/")[0]; - if ( ananames.insert(ananame).second ) addAnalysis(ananame); - aos.push_back(ao); - } - } - if ( !xsec || !sow ) { - MSG_ERROR( "Error in AnalysisHandler::mergeYodas: The file " << file - << " did not contain weights and cross section info."); - exit(1); - } - xsecs.push_back(xsec->point(0).x()); - sows.push_back(sow); - xsecerrs.push_back(sqr(xsec->point(0).xErrAvg())); - _eventCounter->operator+=(*sow); //< HAHAHAHAHA! - sows.push_back(sow); - aosv.push_back(aos); - } catch (...) { //< YODA::ReadError& - throw UserError("Unexpected error in reading file: " + file); - } - } + // // For each file make sure that cross section and sum-of-weights + // // objects are present and stor all RAW ones in a vector; + // vector aos; + // try { + // /// @todo Use new YODA SFINAE to fill the smart ptr vector directly + // vector aos_raw; + // YODA::read(file, aos_raw); + // for (YODA::AnalysisObject* aor : aos_raw) { + // YODA::AnalysisObjectPtr ao(aor); + // if ( ao->path().substr(0, 5) != "/RAW/" ) continue; + // ao->setPath(ao->path().substr(4)); + // if ( ao->path() == "/_XSEC" ) + // xsec = dynamic_pointer_cast(ao); + // else if ( ao->path() == "/_EVTCOUNT" ) + // sow = dynamic_pointer_cast(ao); + // else { + // stripOptions(ao, delopts); + // string ananame = split(ao->path(), "/")[0]; + // if ( ananames.insert(ananame).second ) addAnalysis(ananame); + // aos.push_back(ao); + // } + // } + // if ( !xsec || !sow ) { + // MSG_ERROR( "Error in AnalysisHandler::mergeYodas: The file " << file + // << " did not contain weights and cross section info."); + // exit(1); + // } + // xsecs.push_back(xsec->point(0).x()); + // sows.push_back(sow); + // xsecerrs.push_back(sqr(xsec->point(0).xErrAvg())); + // _eventCounter->operator+=(*sow); //< HAHAHAHAHA! + // sows.push_back(sow); + // aosv.push_back(aos); + // } catch (...) { //< YODA::ReadError& + // throw UserError("Unexpected error in reading file: " + file); + // } + // } - // Now calculate the scale to be applied for all bins in a file - // and get the common cross section and sum of weights. - _xs = _xserr = 0.0; - for ( int i = 0, N = sows.size(); i < N; ++i ) { - double effnent = sows[i]->effNumEntries(); - _xs += (equiv? effnent: 1.0)*xsecs[i]; - _xserr += (equiv? sqr(effnent): 1.0)*xsecerrs[i]; - } + // // Now calculate the scale to be applied for all bins in a file + // // and get the common cross section and sum of weights. + // _xs = _xserr = 0.0; + // for ( int i = 0, N = sows.size(); i < N; ++i ) { + // double effnent = sows[i]->effNumEntries(); + // _xs += (equiv? effnent: 1.0)*xsecs[i]; + // _xserr += (equiv? sqr(effnent): 1.0)*xsecerrs[i]; + // } - vector scales(sows.size(), 1.0); - if ( equiv ) { - _xs /= _eventCounter.effNumEntries(); - _xserr = sqrt(_xserr)/_eventCounter.effNumEntries(); - } else { - _xserr = sqrt(_xserr); - for ( int i = 0, N = sows.size(); i < N; ++i ) - scales[i] = (_eventCounter.sumW()/sows[i]->sumW())*(xsecs[i]/_xs); - } + // vector scales(sows.size(), 1.0); + // if ( equiv ) { + // _xs /= _eventCounter.effNumEntries(); + // _xserr = sqrt(_xserr)/_eventCounter.effNumEntries(); + // } else { + // _xserr = sqrt(_xserr); + // for ( int i = 0, N = sows.size(); i < N; ++i ) + // scales[i] = (_eventCounter.sumW()/sows[i]->sumW())*(xsecs[i]/_xs); + // } - // Initialize the analyses allowing them to book analysis objects. - for (AnaHandle a : _analyses) { - MSG_DEBUG("Initialising analysis: " << a->name()); - if ( !a->info().reentrant() ) - MSG_WARNING("Analysis " << a->name() << " has not been validated to have " - << "a reentrant finalize method. The result is unpredictable."); - try { - // Allow projection registration in the init phase onwards - a->_allowProjReg = true; - cerr << "sqrtS " << sqrtS() << endl; - a->init(); - //MSG_DEBUG("Checking consistency of analysis: " << a->name()); - //a->checkConsistency(); - } catch (const Error& err) { - cerr << "Error in " << a->name() << "::init method: " << err.what() << endl; - exit(1); - } - MSG_DEBUG("Done initialising analysis: " << a->name()); - } - _initialised = true; - // Get a list of all anaysis objects to handle. - map current; - for ( auto ao : getData(false, true) ) current[ao->path()] = ao; - // Go through all objects to be merged and add them to current - // after appropriate scaling. - for ( int i = 0, N = aosv.size(); i < N; ++i) - for ( auto ao : aosv[i] ) { - if ( ao->path() == "/_XSEC" || ao->path() == "_EVTCOUNT" ) continue; - auto aoit = current.find(ao->path()); - if ( aoit == current.end() ) { - MSG_WARNING("" << ao->path() << " was not properly booked."); - continue; - } - if ( !addaos(aoit->second, ao, scales[i]) ) - MSG_WARNING("Cannot merge objects with path " << ao->path() - <<" of type " << ao->annotation("Type") ); - } - // Now we can simply finalize() the analysis, leaving the - // controlling program to write it out some yoda-file. - finalize(); - } + // // Initialize the analyses allowing them to book analysis objects. + // for (AnaHandle a : _analyses) { + // MSG_DEBUG("Initialising analysis: " << a->name()); + // if ( !a->info().reentrant() ) + // MSG_WARNING("Analysis " << a->name() << " has not been validated to have " + // << "a reentrant finalize method. The result is unpredictable."); + // try { + // // Allow projection registration in the init phase onwards + // a->_allowProjReg = true; + // cerr << "sqrtS " << sqrtS() << endl; + // a->init(); + // //MSG_DEBUG("Checking consistency of analysis: " << a->name()); + // //a->checkConsistency(); + // } catch (const Error& err) { + // cerr << "Error in " << a->name() << "::init method: " << err.what() << endl; + // exit(1); + // } + // MSG_DEBUG("Done initialising analysis: " << a->name()); + // } + // _initialised = true; + // // Get a list of all anaysis objects to handle. + // map current; + // for ( auto ao : getData(false, true) ) current[ao->path()] = ao; + // // Go through all objects to be merged and add them to current + // // after appropriate scaling. + // for ( int i = 0, N = aosv.size(); i < N; ++i) + // for ( auto ao : aosv[i] ) { + // if ( ao->path() == "/_XSEC" || ao->path() == "_EVTCOUNT" ) continue; + // auto aoit = current.find(ao->path()); + // if ( aoit == current.end() ) { + // MSG_WARNING("" << ao->path() << " was not properly booked."); + // continue; + // } + // if ( !addaos(aoit->second, ao, scales[i]) ) + // MSG_WARNING("Cannot merge objects with path " << ao->path() + // <<" of type " << ao->annotation("Type") ); + // } + // // Now we can simply finalize() the analysis, leaving the + // // controlling program to write it out some yoda-file. + // finalize(); + // } void AnalysisHandler::readData(const string& filename) { vector aos; try { /// @todo Use new YODA SFINAE to fill the smart ptr vector directly vector aos_raw; - YODA::ReaderYODA::read(filename, aos_raw); + YODA::read(filename, aos_raw); for (YODA::AnalysisObject* aor : aos_raw) aos.push_back(YODA::AnalysisObjectPtr(aor)); //} catch (const YODA::ReadError & e) { } catch (...) { //< YODA::ReadError& throw UserError("Unexpected error in reading file: " + filename); } if (!aos.empty()) addData(aos); } vector AnalysisHandler::getRivetAOs() const { vector rtn; for (AnaHandle a : _analyses) { for (const auto & ao : a->analysisObjects()) { rtn.push_back(ao); } } rtn.push_back(_eventCounter); rtn.push_back(_xs); return rtn; } vector AnalysisHandler::getYodaAOs(bool includeorphans, bool includetmps, bool usefinalized) const { vector rtn; if (usefinalized) rtn = _finalizedAOs; else { for (auto rao : getRivetAOs()) { // need to set the index // before we can search the PATH rao.get()->setActiveWeightIdx(_defaultWeightIdx); // Exclude paths from final write-out if they contain a "TMP" layer (i.e. matching "/TMP/") if (!includetmps && rao->path().find("/TMP/") != string::npos) continue; for (size_t iW = 0; iW < numWeights(); iW++) { rao.get()->setActiveWeightIdx(iW); rtn.push_back(rao.get()->activeYODAPtr()); } } } // Sort histograms alphanumerically by path before write-out sort(rtn.begin(), rtn.end(), [](YODA::AnalysisObjectPtr a, YODA::AnalysisObjectPtr b) { return a->path() < b->path(); } ); return rtn; } vector AnalysisHandler::getData(bool includeorphans, bool includetmps, bool usefinalized) const { return getYodaAOs(includeorphans, includetmps, usefinalized); } void AnalysisHandler::writeData(const string& filename) const { vector out = _finalizedAOs; set finalana; for ( auto ao : out) finalana.insert(ao->path()); out.reserve(2*out.size()); vector aos = getData(false, true, false); if ( _dumping ) { for ( auto ao : aos ) { if ( finalana.find(ao->path()) == finalana.end() ) - out.push_back(AnalysisObjectPtr(ao->newclone())); + out.push_back(YODA::AnalysisObjectPtr(ao->newclone())); } } for ( auto ao : aos ) { ao = YODA::AnalysisObjectPtr(ao->newclone()); ao->setPath("/RAW" + ao->path()); out.push_back(ao); } try { - YODA::WriterYODA::write(filename, aos.begin(), aos.end()); + YODA::write(filename, aos.begin(), aos.end()); } catch (...) { //< YODA::WriteError& throw UserError("Unexpected error in writing file: " + filename); } } string AnalysisHandler::runName() const { return _runname; } size_t AnalysisHandler::numEvents() const { return _eventCounter->numEntries(); } std::vector AnalysisHandler::analysisNames() const { std::vector rtn; for (AnaHandle a : _analyses) { rtn.push_back(a->name()); } return rtn; } const AnaHandle AnalysisHandler::analysis(const std::string& analysisname) const { for (const AnaHandle a : analyses()) if (a->name() == analysisname) return a; throw Error("No analysis named '" + analysisname + "' registered in AnalysisHandler"); } AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector& analysisnames) { for (const string& aname : analysisnames) { //MSG_DEBUG("Adding analysis '" << aname << "'"); addAnalysis(aname); } return *this; } AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector& analysisnames) { for (const string& aname : analysisnames) { removeAnalysis(aname); } return *this; } AnalysisHandler& AnalysisHandler::setCrossSection(double xs, double xserr) { _xs = Scatter1DPtr(weightNames(), Scatter1D("_XSEC")); _eventCounter.get()->setActiveWeightIdx(_defaultWeightIdx); - double nomwgt = sumOfWeights(); + double nomwgt = sumW(); // The cross section of each weight variation is the nominal cross section - // times the sumOfWeights(variation) / sumOfWeights(nominal). + // times the sumW(variation) / sumW(nominal). // This way the cross section will work correctly for (size_t iW = 0; iW < numWeights(); iW++) { _eventCounter.get()->setActiveWeightIdx(iW); - double s = sumOfWeights() / nomwgt; + double s = sumW() / nomwgt; _xs.get()->setActiveWeightIdx(iW); _xs->addPoint(xs*s, xserr*s); } _eventCounter.get()->unsetActiveWeight(); _xs.get()->unsetActiveWeight(); return *this; } AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) { analysis->_analysishandler = this; _analyses.insert(AnaHandle(analysis)); return *this; } PdgIdPair AnalysisHandler::beamIds() const { return Rivet::beamIds(beams()); } double AnalysisHandler::sqrtS() const { return Rivet::sqrtS(beams()); } void AnalysisHandler::setIgnoreBeams(bool ignore) { _ignoreBeams=ignore; } } diff --git a/src/Tools/Correlators.cc b/src/Tools/Correlators.cc --- a/src/Tools/Correlators.cc +++ b/src/Tools/Correlators.cc @@ -1,233 +1,234 @@ // -*- C++ -*- #include "Rivet/Tools/Correlators.hh" namespace Rivet { - + + // Constructor - Correlators::Correlators(const ParticleFinder& fsp, int nMaxIn, - int pMaxIn, vector pTbinEdgesIn) : + Correlators::Correlators(const ParticleFinder& fsp, int nMaxIn, + int pMaxIn, vector pTbinEdgesIn) : nMax(nMaxIn + 1), pMax(pMaxIn + 1), pTbinEdges(pTbinEdgesIn) { setName("Correlators"); declareProjection(fsp, "FS"); isPtDiff = !pTbinEdges.empty(); if (isPtDiff) { vector::iterator underflow = pTbinEdges.begin(); pTbinEdges.insert(underflow,pTbinEdges[0]-1); - } + } setToZero(); } // Alternative constructor. - Correlators::Correlators(const ParticleFinder& fsp, int nMaxIn, + Correlators::Correlators(const ParticleFinder& fsp, int nMaxIn, int pMaxIn, const Scatter2DPtr hIn) : nMax(nMaxIn + 1), pMax(pMaxIn + 1) { for (auto b : hIn->points()) pTbinEdges.push_back(b.xMin()); pTbinEdges.push_back(hIn->points().back().xMax()); setName("Correlators"); declareProjection(fsp, "FS"); isPtDiff = !pTbinEdges.empty(); if (isPtDiff) { vector::iterator underflow = pTbinEdges.begin(); pTbinEdges.insert(underflow,pTbinEdges[0]-1); - } + } setToZero(); } - + // Set all elements in vectors to zero void Correlators::setToZero(){ vector< complex > pTmp(pMax, _ZERO); Vec2D qTmp(nMax, pTmp); qVec = qTmp; if (isPtDiff) { pVec.erase(pVec.begin(), pVec.end()); for (double pT : pTbinEdges) pVec.insert(pair(pT, qVec)); } } // Functions for output: const pair Correlators::intCorrelator(vector n) const { // Create vector of zeros for normalisation and vector of initial // powers int m = n.size(); vector powers(m, 1); vector zeros(m, 0); complex num = recCorr(m, n, powers, false); complex den = recCorr(m, zeros, powers, false); pair ret; ret.second = (den.real() < _TINY) ? 0. : den.real(); ret.first = num.real(); return ret; - } - + } + const vector> Correlators::pTBinnedCorrelators(vector n, bool overflow) const { // Create vector of zeros for normalisation and vector of initial // powers - if (!isPtDiff) + if (!isPtDiff) cout << "You must book the correlator with a binning if you want to" " extract binned correlators! Failing." << endl; int m = n.size(); vector powers(m, 1); vector zeros(m, 0); vector> ret; for (double pT : pTbinEdges){ complex num = recCorr(m, n, powers, true, pT); complex den = recCorr(m, zeros, powers, true, pT); pair tmp; tmp.second = (den.real() < _TINY) ? 0. : den.real(); tmp.first = num.real(); ret.push_back(tmp); } - if (!overflow) + if (!overflow) return vector > (ret.begin() + 1, ret.end() - 1); return ret; - } - + } + // M-particle correlation with eta-gap const pair Correlators::intCorrelatorGap(const Correlators& other, vector n1, vector n2) const { // Create vectors of zeros for normalisation and vectors of initial // powers int m1 = n1.size(); int m2 = n2.size(); // Return if too few particles in event. vector zero1(m1, 0); vector zero2(m2, 0); vector p1(m1, 1); vector p2(m2, 1); complex num1 = recCorr(m1, n1, p1, false); complex den1 = recCorr(m1, zero1, p1, false); complex num2 = other.recCorr(m2, n2, p2, false); complex den2 = other.recCorr(m2, zero2, p2, false); complex num = num1 * num2; complex den = den1 * den2; pair ret; ret.second = (den1.real() < _TINY || den2.real() < _TINY) ? 0. : den.real(); ret.first = num.real(); return ret; - } - + } + // M-particle correlation with eta-gap const vector> Correlators::pTBinnedCorrelatorsGap( const Correlators& other, vector n1, vector n2, bool overflow) const { - if (!isPtDiff) + if (!isPtDiff) cout << "You must book the correlator with a binning if you want to" " extract binned correlators! Failing." << endl; // Create vectors of zeros for normalisation and vectors of initial // powers int m1 = n1.size(); int m2 = n2.size(); vector zero1(m1, 0); vector zero2(m2, 0); vector p1(m1, 1); vector p2(m2, 1); vector> ret; for (double pT : pTbinEdges) { complex num1 = recCorr(m1, n1, p1, true, pT); complex den1 = recCorr(m1, zero1, p1, true, pT); complex num2 = other.recCorr(m2, n2, p2, false); complex den2 = other.recCorr(m2, zero2, p2, false); complex num = num1 * num2; complex den = den1 * den2; pair tmp; - tmp.second = (den1.real() < _TINY || den2.real() < _TINY) + tmp.second = (den1.real() < _TINY || den2.real() < _TINY) ? 0. : den.real(); tmp.first = num.real(); ret.push_back(tmp); } //If we don't want to include underflow/overflow, remove them here - if (!overflow) + if (!overflow) return vector > (ret.begin() + 1, ret.end() - 1); return ret; - } + } // Project function. Loops over array and calculates Q vectors void Correlators::project(const Event& e) { setToZero(); // @TODO: Weight could be implemented to account for non-uniform // acceptance if detector simulation is needed. If not, the weight // can be unity. Note that this weight is not the MC event weight, which // should be used when filling histograms (as usual). const double w = 1.0;; const Particles& parts = applyProjection(e, "FS").particles(); // Check that we have at least two particles in the event if (parts.size() > 2) { for(const Particle& p : parts) fillCorrelators(p, w); } - } - + } + // Calculate correlators from one particle void Correlators::fillCorrelators(const Particle& p, const double& weight = 1.) { - for (int iN = 0; iN < nMax; ++iN) + for (int iN = 0; iN < nMax; ++iN) for (int iP = 0; iP < pMax; ++iP) { double real = cos(iN * p.phi()); double imag = sin(iN * p.phi()); complex expi(real, imag); complex tmp = pow(weight, iP) * expi; qVec[iN][iP] += tmp; if (isPtDiff) { map::iterator pTitr = pVec.lower_bound(p.pT()); // Move to the correct bin. - if (pTitr != pVec.begin()) pTitr--; + if (pTitr != pVec.begin()) pTitr--; pTitr->second[iN][iP] += tmp; } - } - } - + } + } + // Two-particle correlator Eq. (19) p. 6 in Generic Fr. paper. - const complex Correlators::twoPartCorr(int n1, int n2, int p1, + const complex Correlators::twoPartCorr(int n1, int n2, int p1, int p2, double pT, bool useP) const { - complex tmp1 = (!useP) ? getQ(n1, p1) : + complex tmp1 = (!useP) ? getQ(n1, p1) : getP(n1, p1, pT); complex tmp2 = getQ(n2, p2); - complex tmp3 = (!useP) ? getQ(n1+n2, p1+p2) : - getP(n1+n2, p1+p2, pT); + complex tmp3 = (!useP) ? getQ(n1+n2, p1+p2) : + getP(n1+n2, p1+p2, pT); complex sum = tmp1 * tmp2 - tmp3; return sum; } - // Find correlators by recursion. Order = M (# of particles), + // Find correlators by recursion. Order = M (# of particles), // n's are harmonics, p's are the powers of the weights - const complex Correlators::recCorr(int order, vector n, + const complex Correlators::recCorr(int order, vector n, vector p, bool useP, double pT) const { // Sanity checks int nUsed = 0; - for (int i = 0, N = n.size(); i < N; ++i) nUsed += n[i]; + for (int i = 0, N = n.size(); i < N; ++i) nUsed += n[i]; if (nMax < nUsed) cout <<"Requested n = " << nUsed << ", nMax = " << nMax << endl; if (int(p.size()) > pMax) cout << "Requested p = " << p.size() << ", pMax = " << pMax << endl; - // If order is 1, then return Q/p vector (important when dealing + // If order is 1, then return Q/p vector (important when dealing // with gaps and one side has only one particle - if (order < 2) + if (order < 2) return (!useP) ? getQ(n[0], p[0]) : getP(n[0], p[0], pT); // Return 2-p correlator explicitly. if ( order < 3 ) - return twoPartCorr(n[0], n[1], p[0], p[1], pT, useP); + return twoPartCorr(n[0], n[1], p[0], p[1], pT, useP); - // Else find nth order harmonics by recursion + // Else find nth order harmonics by recursion // at order M - 1 int orderNow = order - 1; int nNow = n[orderNow]; int pNow = p[orderNow]; complex recNow = getQ(n[orderNow], p[orderNow]); - recNow *= recCorr(orderNow, n, p, useP, pT); + recNow *= recCorr(orderNow, n, p, useP, pT); for (int i = 0; i < orderNow; ++i){ vector tmpN, tmpP; for (int j = 0; j < orderNow; ++j){ tmpN.push_back(n[j]); tmpP.push_back(p[j]); } tmpN[i] += nNow; tmpP[i] += pNow; recNow -= recCorr(orderNow, tmpN, tmpP, useP, pT); } return recNow; } } diff --git a/src/Tools/Random.cc b/src/Tools/Random.cc --- a/src/Tools/Random.cc +++ b/src/Tools/Random.cc @@ -1,64 +1,66 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include #if defined(_OPENMP) #include "omp.h" #endif namespace Rivet { + using namespace std; + // Return a thread-safe random number generator mt19937& rng() { #if defined(_OPENMP) static map gens; const int nthread = omp_get_thread_num(); if (gens.find(nthread) == gens.end()) { // Make seeds for each thread, either via the standard seed generator or based on a fixed seed from the environment vector seeds(nthread+1); const uint32_t envseed = getEnvParam("RIVET_RANDOM_SEED", 0); //cout << "RIVET_RANDOM_SEED = " << envseed << endl; if (envseed > 0) { std::iota(seeds.begin(), seeds.end(), envseed); } else { seed_seq seq{1,2,3,4,5}; seq.generate(seeds.begin(), seeds.end()); } gens[nthread] = mt19937(seeds[nthread]); //cout << "Thread " << nthread+1 << ", seed=" << seeds[nthread] << " (" << gens.size() << " RNGs)" << endl; } mt19937& g = gens[nthread]; #else static mt19937 g(12345); #endif return g; } // Return a uniformly sampled random number between 0 and 1 double rand01() { const double x = generate_canonical(rng()); ///< @todo What's the "correct" number of bits of randomness? //cout << "RAND01 -> " << x << endl; return x; } // Return a Gaussian/normal sampled random number with the given mean and width double randnorm(double loc, double scale) { normal_distribution<> d(loc, scale); const double x = d(rng()); //cout << "RANDNORM -> " << x << endl; return x; } // Return a log-normal sampled random number double randlognorm(double loc, double scale) { lognormal_distribution<> d(loc, scale); const double x = d(rng()); //cout << "RANDLOGNORM -> " << x << endl; return x; } } diff --git a/src/Tools/RivetYODA.cc b/src/Tools/RivetYODA.cc --- a/src/Tools/RivetYODA.cc +++ b/src/Tools/RivetYODA.cc @@ -1,406 +1,407 @@ #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Tools/RivetYODA.hh" #include "Rivet/Tools/RivetPaths.hh" #include "YODA/ReaderYODA.h" #include "YODA/ReaderAIDA.h" // use execinfo for backtrace if available #include "DummyConfig.hh" #ifdef HAVE_EXECINFO_H #include #endif using namespace std; namespace Rivet { + template Wrapper::~Wrapper() {} template Wrapper::Wrapper(const vector& weightNames, const T & p) { for (const string& weightname : weightNames) { _persistent.push_back(make_shared(p)); auto obj = _persistent.back(); if (weightname != "") obj->setPath(obj->path() + "[" + weightname + "]"); } } template typename T::Ptr Wrapper::active() const { if ( !_active ) { #ifdef HAVE_BACKTRACE void * buffer[4]; backtrace(buffer, 4); backtrace_symbols_fd(buffer, 4 , 1); #endif assert(false && "No active pointer set. Was this object booked in init()?"); } return _active; } template void Wrapper::newSubEvent() { typename TupleWrapper::Ptr tmp = make_shared>(_persistent[0]->clone()); tmp->reset(); _evgroup.push_back( tmp ); _active = _evgroup.back(); assert(_active); } string getDatafilePath(const string& papername) { /// Try to find YODA otherwise fall back to try AIDA const string path1 = findAnalysisRefFile(papername + ".yoda"); if (!path1.empty()) return path1; const string path2 = findAnalysisRefFile(papername + ".aida"); if (!path2.empty()) return path2; throw Rivet::Error("Couldn't find ref data file '" + papername + ".yoda" + " in data path, '" + getRivetDataPath() + "', or '.'"); } map getRefData(const string& papername) { const string datafile = getDatafilePath(papername); // Make an appropriate data file reader and read the data objects /// @todo Remove AIDA support some day... YODA::Reader& reader = (datafile.find(".yoda") != string::npos) ? \ YODA::ReaderYODA::create() : YODA::ReaderAIDA::create(); vector aovec; reader.read(datafile, aovec); // Return value, to be populated map rtn; for ( YODA::AnalysisObject* ao : aovec ) { YODA::AnalysisObjectPtr refdata(ao); if (!refdata) continue; const string plotpath = refdata->path(); // Split path at "/" and only return the last field, i.e. the histogram ID const size_t slashpos = plotpath.rfind("/"); const string plotname = (slashpos+1 < plotpath.size()) ? plotpath.substr(slashpos+1) : ""; rtn[plotname] = refdata; } return rtn; } } namespace { using Rivet::Fill; using Rivet::Fills; using Rivet::TupleWrapper; template double get_window_size(const typename T::Ptr & histo, typename T::BinType x) { // the bin index we fall in const auto binidx = histo->binIndexAt(x); // gaps, overflow, underflow don't contribute if ( binidx == -1 ) return 0; const auto & b = histo->bin(binidx); // if we don't have a valid neighbouring bin, // we use infinite width typename T::Bin b1(-1.0/0.0, 1.0/0.0); // points in the top half compare to the upper neighbour if ( x > b.xMid() ) { size_t nextidx = binidx + 1; if ( nextidx < histo->bins().size() ) b1 = histo->bin(nextidx); } else { // compare to the lower neighbour int nextidx = binidx - 1; if ( nextidx >= 0 ) b1 = histo->bin(nextidx); } // the factor 2 is arbitrary, could poss. be smaller return min( b.width(), b1.width() ) / 2.0; } template typename T::BinType fillT2binT(typename T::FillType a) { return a; } template <> YODA::Profile1D::BinType fillT2binT(YODA::Profile1D::FillType a) { return get<0>(a); } template <> YODA::Profile2D::BinType fillT2binT(YODA::Profile2D::FillType a) { return YODA::Profile2D::BinType{ get<0>(a), get<1>(a) }; } template void commit(vector & persistent, const vector< vector> > & tuple, const vector> & weights ) { // TODO check if all the xs are in the same bin anyway! // Then no windowing needed assert(persistent.size() == weights[0].size()); for ( const auto & x : tuple ) { double maxwindow = 0.0; for ( const auto & xi : x ) { // TODO check for NOFILL here // persistent[0] has the same binning as all the other persistent objects double window = get_window_size(persistent[0], fillT2binT(xi.first)); if ( window > maxwindow ) maxwindow = window; } const double wsize = maxwindow; // all windows have same size set edgeset; // bin edges need to be in here! for ( const auto & xi : x ) { edgeset.insert(fillT2binT(xi.first) - wsize); edgeset.insert(fillT2binT(xi.first) + wsize); } vector< std::tuple,double> > hfill; double sumf = 0.0; auto edgit = edgeset.begin(); double ehi = *edgit; while ( ++edgit != edgeset.end() ) { double elo = ehi; ehi = *edgit; valarray sumw(0.0, persistent.size()); // need m copies of this bool gap = true; // Check for gaps between the sub-windows. for ( size_t i = 0; i < x.size(); ++i ) { // check equals comparisons here! if ( fillT2binT(x[i].first) + wsize >= ehi && fillT2binT(x[i].first) - wsize <= elo ) { sumw += x[i].second * weights[i]; gap = false; } } if ( gap ) continue; hfill.push_back( make_tuple( (ehi + elo)/2.0, sumw, (ehi - elo) ) ); sumf += ehi - elo; } for ( auto f : hfill ) for ( size_t m = 0; m < persistent.size(); ++m ) persistent[m]->fill( get<0>(f), get<1>(f)[m], get<2>(f)/sumf ); // Note the scaling to one single fill } } template<> void commit(vector & persistent, const vector< vector> > & tuple, const vector> & weights) {} template<> void commit(vector & persistent, const vector< vector> > & tuple, const vector> & weights) {} template double distance(T a, T b) { return abs(a - b); } template <> double distance >(tuple a, tuple b) { return Rivet::sqr(get<0>(a) - get<0>(b)) + Rivet::sqr(get<1>(a) - get<1>(b)); } bool copyao(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst) { for (const std::string& a : src->annotations()) dst->setAnnotation(a, src->annotation(a)); if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; if ( aocopy(src,dst) ) return true; return false; } bool addaos(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale) { if ( aoadd(dst,src,scale) ) return true; if ( aoadd(dst,src,scale) ) return true; if ( aoadd(dst,src,scale) ) return true; if ( aoadd(dst,src,scale) ) return true; if ( aoadd(dst,src,scale) ) return true; return false; } } /// fills is a vector of sub-event with an ordered set of x-values of /// the fills in each sub-event. NOFILL should be an "impossible" /// value for this histogram. Returns a vector of sub-events with /// an ordered vector of fills (including NOFILLs) for each sub-event. template vector< vector > > match_fills(const vector::Ptr> & evgroup, const Fill & NOFILL) { vector< vector > > matched; // First just copy subevents into vectors and find the longest vector. unsigned int maxfill = 0; // length of biggest vector int imax = 0; // index position of biggest vector for ( const auto & it : evgroup ) { const auto & subev = it->fills(); if ( subev.size() > maxfill ) { maxfill = subev.size(); imax = matched.size(); } matched.push_back(vector >(subev.begin(), subev.end())); } // Now, go through all subevents with missing fills. const vector> & full = matched[imax]; // the longest one for ( auto & subev : matched ) { if ( subev.size() == maxfill ) continue; // Add NOFILLs to the end; while ( subev.size() < maxfill ) subev.push_back(NOFILL); // Iterate from the back and shift all fill values backwards by // swapping with NOFILLs so that they better match the full // subevent. for ( int i = maxfill - 1; i >= 0; --i ) { if ( subev[i] == NOFILL ) continue; size_t j = i; while ( j + 1 < maxfill && subev[j + 1] == NOFILL && distance(fillT2binT(subev[j].first), fillT2binT(full[j].first)) > distance(fillT2binT(subev[j].first), fillT2binT(full[j + 1].first)) ) { swap(subev[j], subev[j + 1]); ++j; } } } // transpose vector>> result(maxfill,vector>(matched.size())); for (size_t i = 0; i < matched.size(); ++i) for (size_t j = 0; j < maxfill; ++j) result.at(j).at(i) = matched.at(i).at(j); return result; } namespace Rivet { template void Wrapper::pushToPersistent(const vector >& weight) { assert( _evgroup.size() == weight.size() ); // have we had subevents at all? const bool have_subevents = _evgroup.size() > 1; if ( ! have_subevents ) { // simple replay of all tuple entries // each recorded fill is inserted into all persistent weightname histos for ( size_t m = 0; m < _persistent.size(); ++m ) for ( const auto & f : _evgroup[0]->fills() ) _persistent[m]->fill( f.first, f.second * weight[0][m] ); } else { // outer index is subevent, inner index is jets in the event vector>> linedUpXs = match_fills(_evgroup, {typename T::FillType(), 0.0}); commit( _persistent, linedUpXs, weight ); } _evgroup.clear(); _active.reset(); } template <> void Wrapper::pushToPersistent(const vector >& weight) { for ( size_t m = 0; m < _persistent.size(); ++m ) { for ( size_t n = 0; n < _evgroup.size(); ++n ) { for ( const auto & f : _evgroup[n]->fills() ) { _persistent[m]->fill( f.second * weight[n][m] ); } } } _evgroup.clear(); _active.reset(); } template <> void Wrapper::pushToPersistent(const vector >& weight) { _evgroup.clear(); _active.reset(); } template <> void Wrapper::pushToPersistent(const vector >& weight) { _evgroup.clear(); _active.reset(); } template <> void Wrapper::pushToPersistent(const vector >& weight) { _evgroup.clear(); _active.reset(); } // explicitly instantiate all wrappers template class Wrapper; template class Wrapper; template class Wrapper; template class Wrapper; template class Wrapper; template class Wrapper; template class Wrapper; template class Wrapper; }