diff --git a/analyses/pluginMC/TEST.cc b/analyses/pluginMC/TEST.cc --- a/analyses/pluginMC/TEST.cc +++ b/analyses/pluginMC/TEST.cc @@ -1,80 +1,74 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/PrimaryParticles.hh" #include "Rivet/Tools/Correlators.hh" namespace Rivet { class TEST : public CumulantAnalysis { public: /// @name Constructors etc. //@{ /// Constructor TEST() : CumulantAnalysis("TEST") { } //@} - public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { ChargedFinalState cfs(-1.0, 1.0); declare(cfs, "CFS"); ChargedFinalState pp(Cuts::abseta < 2.0); declare(pp, "PP"); h_c22 = bookScatter2D("c22",120,0,120); h_c23 = bookScatter2D("c23",120,0,120); ec22 = bookECorrelator<2,2>("ec22",h_c22); ec23 = bookECorrelator<3,2>("ec32",h_c22); pair max = getMaxValues(); // Declare correlator projections. declare(Correlators(pp, max.first, max.second),"CRS"); } - /// Perform the per-event analysis void analyze(const Event& event) { ec22->fill(apply(event,"CFS").particles().size(), apply(event,"CRS"), event.weight()); ec23->fill(apply(event,"CFS").particles().size(), apply(event,"CRS"), event.weight()); } - - /// Normalise histograms etc., after the run void finalize() { CumulantAnalysis::finalize(); cnTwoInt(h_c22,ec22); } //@} + private: - private: /// @name Histograms //@{ - Scatter2DPtr h_c22; ECorrPtr ec22; Scatter2DPtr h_c23; ECorrPtr ec23; //@} - }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(TEST); } 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,1087 +1,1126 @@ // -*- 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 * correlators. * Classes: * Correlators: Calculates single event correlators of a given harmonic. * 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. * @author Vytautas Vislavicius, Christine O. Rasmussen, Christian Bierlich. */ namespace Rivet { /* @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'' * 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 // 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. // @parm pMaxIn is the maximal number of particles // 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, int pMaxIn = 0, vector pTbinEdgesIn = {}); // Constructor which takes a Scatter2D to estimate bin edges. 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 /// 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 /// 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, 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 /// 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; /// @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); // @brief Compare against other projection. Test if harmonics, pT-bins // and underlying final state are similar. int 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; return mkPCmp(p, "FS"); }; // @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. 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), // 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; // 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, 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] 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, // to ensure access. class CumulantAnalysis : public Analysis { private: // 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. 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. Error errorMethod; /// @brief Base class for correlator bins. class CorBinBase { public: - virtual ~CorBinBase() = 0; + 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 /// Correlator (a pair giving numerator and denominator of _event). void fill(const pair& cor, const double& weight) { _numEntries += 1.; // 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; } 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 /// 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()); } + // 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; }; // End of CorBin sub-class. public: /// @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. //ECorrelator(vector h) : h1(h), h2({}), // 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), reference() {}; /// @brief Fill the appropriate bin given an input (per event) /// observable, eg. centrality. 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 /// 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, 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; - cout << "TEST: " << (*prIn.begin())->effNumEntries() << endl; + } + + /// @brief Fill bins with content from preloaded histograms. + void fillFromProfs() { + list::iterator hItr = profs.begin(); + 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()); + } + } // 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 // 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. 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) - eCorrProfs.push_back(bookProfile1D(name+"-e"+to_string(i),*hIn)); + for (int i = 0; i < BOOT_BINS; ++i) { + Profile1DPtr tmp = bookProfile1D(name+"-"+to_string(i),*hIn); + tmp->setPath(this->name()+"/CORR/" + name+"-"+to_string(i)); + //tmp->setPath(tmp->path()+"CORR"); + 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, 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) - eCorrProfs.push_back(bookProfile1D(name+"-e"+to_string(i),*hIn)); + for (int i = 0; i < BOOT_BINS; ++i) { + Profile1DPtr tmp = bookProfile1D(name+"-"+to_string(i),*hIn); + tmp->setPath(this->name()+"/CORR/" + name+"-"+to_string(i)); + //tmp->setPath(tmp->path()+"CORR"); + 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, 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. 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. template ECorrPtr bookECorrelatorGap(const string name, const Scatter2DPtr hIn) { return bookECorrelatorGap(name, Correlators::hVec(N, M), hIn); } // @brief Finalize MUST be explicitly called for this base class with a - // CumulantsAnalysis::finalize() call as the first thing in each analysis, - // in order to stream the contents of ECorrelators to a yoda file for - // reentrant finalize. + // CumulantsAnalysis::finalize() call as the first thing in each analysis' + // finalize step, in order to stream the contents of ECorrelators to a yoda + // file for reentrant finalize. void finalize() { - for (auto ecItr = eCorrPtrs.begin(); ecItr != eCorrPtrs.end(); ++ecItr) + 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. list eCorrPtrs; public: // @brief Constructor. Use CumulantAnalysis as base class for the // analysis to have access to functionality. CumulantAnalysis (string n) : Analysis(n), errorMethod(ENVELOPE) {}; // @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 // 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 // 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 // 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) 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 /// @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(), b.xErrPlus(), yemin, yemax )); } } hOut->reset(); hOut->points().clear(); 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 /// @parm k constant below the root. 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(), 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(), 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. 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. 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(); // 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 = 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 // 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(), 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. + (*hItr)->setTotalDbn(YODA::Dbn2D(ne,sow,sow2,0.,0.,0.,0.,0.)); + for (int j = 0, N = profBins.size(); j < N; ++j) { (*hItr)->addBin(profBins[j]); } } // 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); } // @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, 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() + 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() || 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() * e2binPtrs[i]->mean() - 18. * e42 + 144. * e4binPtrs[i]->mean()*e22 - 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, 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() - 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/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc --- a/src/Core/AnalysisHandler.cc +++ b/src/Core/AnalysisHandler.cc @@ -1,571 +1,574 @@ // -*- 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" namespace Rivet { AnalysisHandler::AnalysisHandler(const string& runname) : _runname(runname), _eventcounter("/_EVTCOUNT"), _xs(NAN), _xserr(NAN), _initialised(false), _ignoreBeams(false), _dumpPeriod(0), _dumping(false) { } AnalysisHandler::~AnalysisHandler() { } Log& AnalysisHandler::getLog() const { return Log::getLog("Rivet.Analysis.Handler"); } void AnalysisHandler::init(const GenEvent& ge) { if (_initialised) throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!"); setRunBeams(Rivet::beams(ge)); MSG_DEBUG("Initialising the analysis handler"); _eventcounter.reset(); // Check that analyses are beam-compatible, and remove those that aren't const size_t num_anas_requested = analysisNames().size(); vector anamestodelete; for (const AnaHandle a : _analyses) { if (!_ignoreBeams && !a->isCompatible(beams())) { //MSG_DEBUG(a->name() << " requires beams " << a->requiredBeams() << " @ " << a->requiredEnergies() << " GeV"); anamestodelete.push_back(a->name()); } } for (const string& aname : anamestodelete) { MSG_WARNING("Analysis '" << aname << "' is incompatible with the provided beams: removing"); removeAnalysis(aname); } if (num_anas_requested > 0 && analysisNames().empty()) { cerr << "All analyses were incompatible with the first event's beams\n" << "Exiting, since this probably wasn't intentional!" << endl; exit(1); } // Warn if any analysis' status is not unblemished for (const AnaHandle a : analyses()) { if (toUpper(a->status()) == "PRELIMINARY") { MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!"); } else if (toUpper(a->status()) == "OBSOLETE") { MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!"); } else if (toUpper(a->status()).find("UNVALIDATED") != string::npos) { MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!"); } } // Initialize the remaining analyses for (AnaHandle a : _analyses) { MSG_DEBUG("Initialising analysis: " << a->name()); try { // Allow projection registration in the init phase onwards a->_allowProjReg = true; a->init(); //MSG_DEBUG("Checking consistency of analysis: " << a->name()); //a->checkConsistency(); } catch (const Error& err) { cerr << "Error in " << a->name() << "::init method: " << err.what() << endl; exit(1); } MSG_DEBUG("Done initialising analysis: " << a->name()); } _initialised = true; MSG_DEBUG("Analysis handler initialised"); } void AnalysisHandler::analyze(const GenEvent& ge) { // Call init with event as template if not already initialised if (!_initialised) init(ge); assert(_initialised); // Ensure that beam details match those from the first event (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); // Weights /// @todo Drop this / just report first weight when we support multiweight events _eventcounter.fill(event.weight()); MSG_DEBUG("Event #" << _eventcounter.numEntries() << " weight = " << event.weight()); // Cross-section #ifdef HEPMC_HAS_CROSS_SECTION if (ge.cross_section()) { _xs = ge.cross_section()->cross_section(); _xserr = ge.cross_section()->cross_section_error(); } #endif // Run the analyses for (AnaHandle a : _analyses) { MSG_TRACE("About to run analysis " << a->name()); try { a->analyze(event); } catch (const Error& err) { cerr << "Error in " << a->name() << "::analyze method: " << err.what() << endl; exit(1); } MSG_TRACE("Finished running analysis " << a->name()); } 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; // First we make copies of all analysis objects. map backupAOs; for (auto ao : getData(false, true) ) backupAOs[ao->path()] = AnalysisObjectPtr(ao->newclone()); // Now we run the (re-entrant) finalize() functions for all analyses. MSG_INFO("Finalising analyses"); for (AnaHandle a : _analyses) { a->setCrossSection(_xs); try { if ( !_dumping || a->info().reentrant() ) a->finalize(); else if ( _dumping ) 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 : getData() ) _finalizedAOs.push_back(AnalysisObjectPtr(ao->newclone())); for ( auto ao : getData(false, true) ) { + // TODO: This should be possible to do in a nicer way, with a flag etc. + if (ao->path().find("/CORR") != 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 = _eventcounter.numEntries(); MSG_INFO("Processed " << nevts << " event" << (nevts != 1 ? "s" : "")); // // Delete analyses // MSG_DEBUG("Deleting analyses"); // _analyses.clear(); // Print out MCnet boilerplate cout << endl; cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << endl; cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:1003.0694 (http://arxiv.org/abs/1003.0694)" << endl; } 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 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); a->addAnalysisObject(ao); /// @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(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(); // First scan all files and extract analysis objects and add the // corresponding anayses.. for ( auto 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 (AnalysisObject* aor : aos_raw) { AnalysisObjectPtr ao = AnalysisObjectPtr(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()); - xsecerrs.push_back(sqr(xsec->point(0).xErrAvg())); + sows.push_back(sow); + xsecerrs.push_back(sqr(xsec->point(0).xErrAvg())); _eventcounter += *sow; 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]; } 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()); + 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::read(filename, aos_raw); for (AnalysisObject* aor : aos_raw) aos.push_back(AnalysisObjectPtr(aor)); } catch (...) { //< YODA::ReadError& throw UserError("Unexpected error in reading file: " + filename); } if (!aos.empty()) addData(aos); } vector AnalysisHandler:: getData(bool includeorphans, bool includetmps) const { vector rtn; // Event counter rtn.push_back( make_shared(_eventcounter) ); // Cross-section + err as scatter YODA::Scatter1D::Points pts; pts.insert(YODA::Point1D(_xs, _xserr)); rtn.push_back( make_shared(pts, "/_XSEC") ); // Analysis histograms for (const AnaHandle a : analyses()) { vector aos = a->analysisObjects(); // MSG_WARNING(a->name() << " " << aos.size()); for (const AnalysisObjectPtr ao : aos) { // Exclude paths from final write-out if they contain a "TMP" layer (i.e. matching "/TMP/") /// @todo This needs to be much more nuanced for re-entrant histogramming if ( !includetmps && ao->path().find("/TMP/" ) != string::npos) continue; rtn.push_back(ao); } } // Sort histograms alphanumerically by path before write-out sort(rtn.begin(), rtn.end(), [](AnalysisObjectPtr a, AnalysisObjectPtr b) {return a->path() < b->path();}); if ( includeorphans ) rtn.insert(rtn.end(), _orphanedPreloads.begin(), _orphanedPreloads.end()); return rtn; } void AnalysisHandler::writeData(const string& filename) const { vector out = _finalizedAOs; out.reserve(2*out.size()); vector aos = getData(false, true); for ( auto ao : aos ) { ao = AnalysisObjectPtr(ao->newclone()); ao->setPath("/RAW" + ao->path()); out.push_back(ao); } - + try { YODA::write(filename, out.begin(), out.end()); } catch (...) { //< YODA::WriteError& throw UserError("Unexpected error in writing file: " + filename); } } std::vector AnalysisHandler::analysisNames() const { std::vector rtn; for (AnaHandle a : _analyses) { rtn.push_back(a->name()); } return rtn; } const AnaHandle AnalysisHandler::analysis(const std::string& analysisname) const { for (const AnaHandle a : analyses()) if (a->name() == analysisname) return a; throw Error("No analysis named '" + analysisname + "' registered in AnalysisHandler"); } AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector& analysisnames) { for (const string& aname : analysisnames) { //MSG_DEBUG("Adding analysis '" << aname << "'"); addAnalysis(aname); } return *this; } AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector& analysisnames) { for (const string& aname : analysisnames) { removeAnalysis(aname); } return *this; } bool AnalysisHandler::needCrossSection() const { bool rtn = false; for (const AnaHandle a : _analyses) { if (!rtn) rtn = a->needsCrossSection(); if (rtn) break; } return rtn; } AnalysisHandler& AnalysisHandler::setCrossSection(double xs) { _xs = xs; return *this; } bool AnalysisHandler::hasCrossSection() const { return (!std::isnan(crossSection())); } AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) { analysis->_analysishandler = this; _analyses.insert(AnaHandle(analysis)); return *this; } PdgIdPair AnalysisHandler::beamIds() const { return Rivet::beamIds(beams()); } double AnalysisHandler::sqrtS() const { return Rivet::sqrtS(beams()); } void AnalysisHandler::setIgnoreBeams(bool ignore) { _ignoreBeams=ignore; } }