rivet is hosted by Hepforge, IPPP Durham
Rivet 4.0.0
Correlators.hh
1// -*- C++ -*-
2#ifndef RIVET_Correlators_HH
3#define RIVET_Correlators_HH
4
5// Tools for calculating flow coefficients using correlators.
6//
7// Classes:
8// Correlators: Calculates single event correlators of a given harmonic.
9// Cumulants: An additional base class for flow analyses
10// Use as: class MyAnalysis : public Analysis, Cumulants {};
11// Includes a framework for calculating cumulants and flow coefficients
12// from single event correlators, including automatic handling of
13// statistical errors. Contains multiple internal sub-classes:
14// CorBinBase: Base class for correlators binned in event or particle observables.
15// CorSingleBin: A simple bin for correlators.
16// CorBin: Has the interface of a simple bin, but does automatic calculation
17// of statistical errors by a bootstrap method.
18// ECorrelator: Data type for event averaged correlators.
19//
21
22#include "Rivet/Analysis.hh"
23#include "Rivet/Projection.hh"
24#include "Rivet/Projections/ParticleFinder.hh"
25#include "YODA/Scatter.h"
26
27namespace Rivet {
28
29
37 class Correlators : public Projection {
38 public:
39
53 Correlators(const ParticleFinder& fsp, int nMaxIn = 2,
54 int pMaxIn = 0, vector<double> pTbinEdgesIn = {});
55
56 // Constructor which takes an Estimate1D to estimate bin edges.
57 Correlators(const ParticleFinder& fsp, int nMaxIn,
58 int pMaxIn, const YODA::Estimate1D hIn);
59
61 using Projection::operator =;
62
63
71 const pair<double,double> intCorrelator(vector<int> n) const;
72
77 const vector<pair<double,double>> pTBinnedCorrelators(vector<int> n,
78 bool overflow = false) const;
79
88 const pair<double,double> intCorrelatorGap(const Correlators& other,
89 vector<int> n1, vector<int> n2) const;
90
102 const vector<pair<double,double>>
103 pTBinnedCorrelatorsGap(const Correlators& other, vector<int> n1, vector<int> n2, bool overflow=false) const;
104
109 static vector<int> hVec(int n, int m) {
110 if (m % 2 != 0) {
111 cout << "Harmonic Vector: Number of particles must be an even number." << endl;
112 return {};
113 }
114 vector<int> ret = {};
115 for (int i = 0; i < m; ++i) {
116 if (i < m/2) ret.push_back(n);
117 else ret.push_back(-n);
118 }
119 return ret;
120 }
121
124 static pair<int,int> getMaxValues(vector< vector<int> >& hList){
125 int nMax = 0, pMax = 0;
126 for (vector<int> h : hList) {
127 int tmpN = 0, tmpP = 0;
128 for ( int i = 0; i < int(h.size()); ++i) {
129 tmpN += abs(h[i]);
130 ++tmpP;
131 }
132 if (tmpN > nMax) nMax = tmpN;
133 if (tmpP > pMax) pMax = tmpP;
134 }
135 return make_pair(nMax,pMax);
136 }
137
138 // Clone on the heap.
139 RIVET_DEFAULT_PROJ_CLONE(Correlators);
140
141
142 protected:
143
145 void project(const Event& e);
146
148 CmpState compare(const Projection& p) const {
149 const Correlators* other = dynamic_cast<const Correlators*>(&p);
150 if (nMax != other->nMax) return CmpState::NEQ;
151 if (pMax != other->pMax) return CmpState::NEQ;
152 if (pTbinEdges != other->pTbinEdges) return CmpState::NEQ;
153 return mkPCmp(p, "FS");
154 }
155
157 void fillCorrelators(const Particle& p, const double& weight);
158
160 const complex<double> getQ(int n, int p) const {
161 bool isNeg = (n < 0);
162 if (isNeg) return conj( qVec[abs(n)][p] );
163 else return qVec[n][p];
164 }
165
167 const complex<double> getP(int n, int p, double pT = 0.) const {
168 bool isNeg = (n < 0);
169 map<double, Vec2D>::const_iterator pTitr = pVec.lower_bound(pT);
170 if (pTitr == pVec.end()) return DBL_NAN;
171 if (isNeg) return conj( pTitr->second[abs(n)][p] );
172 else return pTitr->second[n][p];
173 }
174
175
176 private:
177
181 const complex<double> recCorr(int order, vector<int> n,
182 vector<int> p, bool useP, double pT = 0.) const;
183
188 const complex<double> twoPartCorr(int n1, int n2, int p1 = 1,
189 int p2 = 1, double pT = 0., bool useP = false) const;
190
192 void setToZero();
193
195 const complex<double> _ZERO = {0., 0.};
196 const double _TINY = 1e-10;
197
199 typedef vector< vector<complex<double>> > Vec2D;
200
202 Vec2D qVec; // Q[n][p]
203 map<double, Vec2D> pVec; // p[pT][n][p]
204
206 int nMax, pMax;
207
209 vector<double> pTbinEdges;
210
211 bool isPtDiff;
212
213 };
214
215
216
224 class CumulantAnalysis : public Analysis {
225 private:
226
227 // Number of bins used for bootstrap calculation of statistical
228 // uncertainties. It is hard coded, and shout NOT be changed unless there
229 // are good reasons to do so.
230 static const int BOOT_BINS = 9;
231
232 // Enum for choosing the method of error calculation.
233 enum Error {
234 VARIANCE,
235 ENVELOPE
236 };
237
238 // The desired error method. Can be changed in the analysis constructor
239 // by setting it appropriately.
240 Error errorMethod;
241
242
244 class CorBinBase {
245 public:
246 CorBinBase() {}
247 virtual ~CorBinBase() {};
248 // Derived class should have fill and mean defined.
249 virtual void fill(const pair<double, double>& cor, const double& weight = 1.0) = 0;
250 virtual double mean() const = 0;
251 };
252
253
259 class CorSingleBin : public CorBinBase {
260 public:
261
263 CorSingleBin()
264 : _sumWX(0.), _sumW(0.), _sumW2(0.), _numEntries(0.)
265 { }
266
267 // Destructor does nothing but must be implemented (?)
268 ~CorSingleBin() {}
269
273 void fill(const pair<double, double>& cor, const double& weight = 1.0) {
274 // Test if denominator for the single event average is zero.
275 if (cor.second < 1e-10) return;
276 // The single event average <M> is then cor.first / cor.second.
277 // With weights this becomes just:
278 _sumWX += cor.first * weight;
279 _sumW += weight * cor.second;
280 _sumW2 += weight * weight * cor.second * cor.second;
281 _numEntries += 1.;
282 }
283
285 double mean() const {
286 if (_sumW < 1e-10) return 0;
287 return _sumWX / _sumW;
288 }
289
291 double sumW() const {
292 return _sumW;
293 }
294
296 double sumW2() const {
297 return _sumW2;
298 }
299
301 double sumWX() const {
302 return _sumWX;
303 }
304
306 double numEntries() const {
307 return _numEntries;
308 }
309
311 void addContent(double ne, double sw, double sw2, double swx) {
312 _numEntries += ne;
313 _sumW += sw;
314 _sumW2 += sw2;
315 _sumWX += swx;
316 }
317
318
319 private:
320
321 double _sumWX, _sumW, _sumW2, _numEntries;
322
323 };
324
325
330 class CorBin : public CorBinBase {
331 public:
332
338 CorBin() : binIndex(0), nBins(BOOT_BINS) {
339 for(size_t i = 0; i < nBins; ++i) bins.push_back(CorSingleBin());
340 }
341
342 // Destructor does nothing but must be implemented (?)
343 ~CorBin() {}
344
346 void fill(const pair<double, double>& cor, const double& weight = 1.0) {
347 // Test if denominator for the single event average is zero.
348 if (cor.second < 1e-10) return;
349 // Fill the correct bin.
350 bins[binIndex].fill(cor, weight);
351 if (binIndex == nBins - 1) binIndex = 0;
352 else ++binIndex;
353 }
354
356 double mean() const {
357 double sow = 0;
358 double sowx = 0;
359 for (auto b : bins) {
360 if (b.sumW() < 1e-10) continue;
361 sow += b.sumW();
362 sowx += b.sumWX();
363 }
364 return sowx / sow;
365 }
366
368 const vector<CorSingleBin>& getBins() const {
369 return bins;
370 }
371
373 template<class T=CorBinBase>
374 vector<T*> getBinPtrs() {
375 vector<T*> ret(bins.size());
376 transform(bins.begin(), bins.end(), ret.begin(), [](CorSingleBin& b) {return &b;});
377 return ret;
378 }
379
381 template<class T=CorBinBase>
382 vector<const T*> getBinPtrs() const {
383 vector<const T*> ret(bins.size());
384 transform(bins.begin(), bins.end(), ret.begin(), [](const CorSingleBin& b) {return &b;});
385 return ret;
386 }
387
388 private:
389
390 vector<CorSingleBin> bins;
391 size_t binIndex;
392 size_t nBins;
393
394 };
395
396
397 public:
398
404 public:
405
411 //ECorrelator(vector<int> h) : h1(h), h2({}),
412 // binX(0), binContent(0), reference() {
413 //};
414
420 ECorrelator(const vector<int>& h, const vector<double>& binIn)
421 : h1(h), h2({}), binX(binIn), binContent(binIn.size() - 1), reference()
422 { }
423
428 ECorrelator(const vector<int>& h1In, const vector<int>& h2In, const vector<double>& binIn)
429 : h1(h1In), h2(h2In), binX(binIn), binContent(binIn.size() - 1), reference()
430 { }
431
433 void fill(const double& obs, const Correlators& c, double weight=1.0) {
434 int index = getBinIndex(obs);
435 if (index < 0) return;
436 binContent[index].fill(c.intCorrelator(h1), weight);
437 }
438
442 void fill(const double& obs, const Correlators& c1, const Correlators& c2, double weight=1.0) {
443 if (!h2.size()) {
444 cout << "Trying to fill gapped correlator, but harmonics behind the gap (h2) are not given!" << endl;
445 return;
446 }
447 int index = getBinIndex(obs);
448 if (index < 0) return;
449 binContent[index].fill(c1.intCorrelatorGap(c2, h1, h2), weight);
450 }
451
456 void fill(const Correlators& c, const double weight=1.0) {
457 vector< pair<double, double> > diffCorr = c.pTBinnedCorrelators(h1);
458 // We always skip overflow when calculating the all-event average.
459 if (diffCorr.size() != binX.size() - 1)
460 cout << "Tried to fill event with wrong binning (ungapped)" << endl;
461 for (size_t i = 0; i < diffCorr.size(); ++i) {
462 int index = getBinIndex(binX[i]);
463 if (index < 0) return;
464 binContent[index].fill(diffCorr[i], weight);
465 }
466 reference.fill(c.intCorrelator(h1), weight);
467 }
468
473 void fill(const Correlators& c1, const Correlators& c2, const double weight = 1.0) {
474 if (!h2.size()) {
475 cout << "Trying to fill gapped correlator, but harmonics behind "
476 "the gap (h2) are not given!" << endl;
477 return;
478 }
479 vector< pair<double, double> > diffCorr = c1.pTBinnedCorrelatorsGap(c2, h1, h2);
480 // We always skip overflow when calculating the all event average.
481 if (diffCorr.size() != binX.size() - 1)
482 cout << "Tried to fill event with wrong binning (gapped)" << endl;
483 for (size_t i = 0; i < diffCorr.size(); ++i) {
484 int index = getBinIndex(binX[i]);
485 if (index < 0) return;
486 binContent[index].fill(diffCorr[i], weight);
487 }
488 reference.fill(c1.intCorrelatorGap(c2, h1, h2), weight);
489 }
490
492 const vector<CorBin>& getBins() const {
493 return binContent;
494 }
495
497 vector<const CorBinBase*> getBinPtrs() const {
498 vector<const CorBinBase*> ret(binContent.size());
499 transform(binContent.begin(), binContent.end(), ret.begin(), [](const CorBin& b) {return &b;});
500 return ret;
501 }
502
504 const vector<double>& getBinX() const {
505 return binX;
506 }
507
509 const vector<int>& getH1() const {
510 return h1;
511 }
512
514 const vector<int>& getH2() const {
515 return h2;
516 }
517
519 void setReference(CorBin refIn) {
520 reference = refIn;
521 }
522
524 CorBin getReference() const {
525 if (reference.mean() < 1e-10)
526 cout << "Warning: ECorrelator, reference bin is zero." << endl;
527 return reference;
528 }
529
532 void setProfs(vector<string> prIn) {
533 profs = prIn;
534 }
535
537 bool fillFromProfile(YODA::AnalysisObjectPtr yao, string name) {
538 auto refs = reference.getBinPtrs<CorSingleBin>();
539 for (size_t i = 0; i < profs.size(); ++i) {
540 if (yao->path() == "/RAW/"+name+"/TMP/"+profs[i]) {
541 YODA::Profile1DPtr pPtr = dynamic_pointer_cast<YODA::Profile1D>(yao);
542 for (size_t j = 0; j < binX.size() - 1; ++j) {
543 const YODA::Dbn2D& pBin = pPtr->binAt(binX[j]);
544 auto tmp = binContent[j].getBinPtrs<CorSingleBin>();
545 tmp[i]->addContent(pBin.numEntries(), pBin.sumW(), pBin.sumW2(),
546 pBin.sumWY());
547 }
548 // Get the reference flow from the underflow bin of the histogram.
549 const YODA::Dbn2D& uBin = pPtr->bin(0);
550 refs[i]->addContent(uBin.numEntries(), uBin.sumW(), uBin.sumW2(),
551 uBin.sumWY());
552 return true;
553 }
554 } // End loop of bootstrapped correlators.
555 return false;
556 }
557
558 private:
559
560 // Get correct bin index for a given @param obs value
561 int getBinIndex(const double& obs) const {
562 // Find the correct index of binContent.
563 // If we are in overflow, just skip.
564 if (obs >= binX.back()) return -1;
565 // If we are in underflow, ditto.
566 if (obs < binX[0]) return -1;
567 int index = 0;
568 for (int i = 0, N = binX.size() - 1; i < N; ++i, ++index)
569 if (obs >= binX[i] && obs < binX[i + 1]) break;
570 return index;
571 }
572
573 // The harmonics vectors.
574 vector<int> h1;
575 vector<int> h2;
576
577 // The bins.
578 vector<double> binX;
579 vector<CorBin> binContent;
580 // The reference flow.
581 CorBin reference;
582 public:
583 // The profile histograms associated with the CorBins for streaming.
584 vector <string> profs;
585
586 };
587
588
590 const pair<int, int> getMaxValues() const {
591 vector< vector<int>> harmVecs;
592 for ( auto eItr = eCorrPtrs.begin(); eItr != eCorrPtrs.end(); ++eItr) {
593 const vector<int>& h1 = (*eItr)->getH1();
594 const vector<int>& h2 = (*eItr)->getH2();
595 if (h1.size() > 0) harmVecs.push_back(h1);
596 if (h2.size() > 0) harmVecs.push_back(h2);
597 }
598 if (harmVecs.size() == 0) {
599 cout << "Warning: You tried to extract max values from harmonic "
600 "vectors, but have not booked any." << endl;
601 return pair<int,int>();
602 }
603 return Correlators::getMaxValues(harmVecs);
604 }
605
607 typedef shared_ptr<ECorrelator> ECorrPtr;
608
611 ECorrPtr bookECorrelator(const string name, const vector<int>& h, const YODA::Estimate1D& hIn) {
612 vector<double> binIn;
613 const YODA::Scatter2D s = hIn.mkScatter();
614 for (const auto& p : s.points()) binIn.push_back(p.xMin());
615 binIn.push_back(s.points().back().xMax());
616 return bookECorrelator(name, h, binIn);
617 }
618
621 ECorrPtr bookECorrelator(const string name, const vector<int>& h, const vector<double>& binIn) {
622 ECorrPtr ecPtr = ECorrPtr(new ECorrelator(h, binIn));
623 vector<string> eCorrProfs;
624 for (int i = 0; i < BOOT_BINS; ++i) {
625 Profile1DPtr tmp;
626 book(tmp,"TMP/"+name+"-"+to_string(i),binIn);
627 eCorrProfs.push_back(name+"-"+to_string(i));
628 }
629 ecPtr->setProfs(eCorrProfs);
630 eCorrPtrs.push_back(ecPtr);
631 return ecPtr;
632 }
633
636 ECorrPtr bookECorrelator(const string name, const vector<int>& h1,
637 const vector<int>& h2, const vector<double>& binIn) {
638 ECorrPtr ecPtr = ECorrPtr(new ECorrelator(h1, h2, binIn));
639 vector<string> eCorrProfs;
640 Profile1DPtr tmp;
641 for (int i = 0; i < BOOT_BINS; ++i) {
642 book(tmp, "TMP/"+name+"-"+to_string(i), binIn);
643 eCorrProfs.push_back(name+"-"+to_string(i));
644 }
645 ecPtr->setProfs(eCorrProfs);
646 eCorrPtrs.push_back(ecPtr);
647 return ecPtr;
648 }
649
652 ECorrPtr bookECorrelator(const string& name, const vector<int>& h1,
653 const vector<int>& h2, const YODA::Estimate1D& hIn) {
654 vector<double> binIn;
655 const YODA::Scatter2D s = hIn.mkScatter();
656 for (const auto& p : s.points()) binIn.push_back(p.xMin());
657 binIn.push_back(s.points().back().xMax());
658 return bookECorrelator(name, h1, h2, binIn);
659 }
660
665 ECorrPtr bookECorrelatorGap(const string& name, const vector<int>& h,
666 const YODA::Estimate1D& hIn) {
667 const vector<int> h1(h.begin(), h.begin() + h.size() / 2);
668 const vector<int> h2(h.begin() + h.size() / 2, h.end());
669 return bookECorrelator(name, h1, h2, hIn);
670 }
671
676 template<unsigned int N, unsigned int M>
677 ECorrPtr bookECorrelator(const string& name, const vector<double>& binIn) {
678 return bookECorrelator(name, Correlators::hVec(N, M), binIn);
679 }
680
685 template<unsigned int N, unsigned int M>
686 ECorrPtr bookECorrelator(const string& name, const YODA::Estimate1D& hIn) {
687 return bookECorrelator(name, Correlators::hVec(N, M), hIn);
688 }
689
694 template<unsigned int N, unsigned int M>
695 ECorrPtr bookECorrelatorGap(const string& name, const YODA::Estimate1D& hIn) {
696 const vector<int> h = Correlators::hVec(N,M);
697 const vector<int> h1(h.begin(), h.begin() + h.size() / 2);
698 const vector<int> h2(h.begin() + h.size() / 2, h.end());
699 return bookECorrelator(name, h1, h2, hIn);
700
701 }
702
703 protected:
704
705 // Bookkeeping of the event averaged correlators.
706 list<ECorrPtr> eCorrPtrs;
707
708
709 public:
710
715 CumulantAnalysis(const string& n)
716 : Analysis(n), errorMethod(VARIANCE)
717 { }
718
726 template<typename T>
727 static void fillScatter(Scatter2DPtr h, const vector<double>& binx, const T& func) {
728 vector<YODA::Point2D> points;
729 // Test if we have proper bins from a booked histogram.
730 bool hasBins = (h->points().size() > 0);
731 for (int i = 0, N = binx.size() - 1; i < N; ++i) {
732 double xMid = (binx[i] + binx[i + 1]) / 2.0;
733 double xeMin = fabs(xMid - binx[i]);
734 double xeMax = fabs(xMid - binx[i + 1]);
735 if (hasBins) {
736 xMid = h->points()[i].x();
737 xeMin = h->points()[i].xErrMinus();
738 xeMax = h->points()[i].xErrPlus();
739 }
740 double yVal = func(i);
741 if (std::isnan(yVal)) yVal = 0.;
742 double yErr = 0;
743 points.push_back(YODA::Point2D(xMid, yVal, xeMin, xeMax, yErr, yErr));
744 }
745 h->reset();
746 h->points().clear();
747 for (int i = 0, N = points.size(); i < N; ++i) h->addPoint(points[i]);
748 }
749
757 template<typename F>
758 void fillScatter(Scatter2DPtr h, const vector<double>& binx, const F func,
759 vector<pair<double, double> >& yErr) const {
760 vector<YODA::Point2D> points;
761 // Test if we have proper bins from a booked histogram.
762 const bool hasBins = (h->points().size() > 0);
763 for (int i = 0, N = binx.size() - 1; i < N; ++i) {
764 double xMid = (binx[i] + binx[i + 1]) / 2.0;
765 double xeMin = fabs(xMid - binx[i]);
766 double xeMax = fabs(xMid - binx[i + 1]);
767 if (hasBins) {
768 xMid = h->points()[i].x();
769 xeMin = h->points()[i].xErrMinus();
770 xeMax = h->points()[i].xErrPlus();
771 }
772 const double yVal = func(i);
773 if (std::isnan(yVal))
774 points.push_back(YODA::Point2D(xMid, 0., xeMin, xeMax,0., 0.));
775 else
776 points.push_back(YODA::Point2D(xMid, yVal, xeMin, xeMax,
777 yErr[i].first, yErr[i].second));
778 }
779 h->reset();
780
781 for (int i = 0, N = points.size(); i < N; ++i) {
782 h->addPoint(points[i]);
783 }
784 }
785
786
790 static void nthPow(Scatter2DPtr hOut, const Scatter2DPtr hIn, const double& n, const double& k = 1.0) {
791 if (n == 0 || n == 1) {
792 cout << "Error: Do not take the 0th or 1st power of a Scatter2D,"
793 " use scale instead." << endl;
794 return;
795 }
796 if (hIn->points().size() != hOut->points().size()) {
797 cout << "nthRoot: Scatterplots: " << hIn->name() << " and " <<
798 hOut->name() << " not initialized with same length." << endl;
799 return;
800 }
801 vector<YODA::Point2D> points;
802 // The error pre-factor is k^(1/n) / n by Taylors formula.
803 double eFac = pow(k,1./n)/n;
804 for (auto b : hIn->points()) {
805 double yVal = pow(k * b.y(),n);
806 if (std::isnan(yVal))
807 points.push_back(YODA::Point2D(b.x(), 0., b.xErrMinus(),
808 b.xErrPlus(), 0, 0 ));
809 else {
810 double yemin = abs(eFac * pow(yVal,1./(n - 1.))) * b.yErrMinus();
811 if (std::isnan(yemin)) yemin = b.yErrMinus();
812 double yemax = abs(eFac * pow(yVal,1./(n - 1.))) * b.yErrPlus();
813 if (std::isnan(yemax)) yemax = b.yErrPlus();
814 points.push_back(YODA::Point2D(b.x(), yVal, b.xErrMinus(),
815 b.xErrPlus(), yemin, yemax ));
816 }
817 }
818 hOut->reset();
819 for (int i = 0, N = points.size(); i < N; ++i)
820 hOut->addPoint(points[i]);
821 }
822
823
827 static void nthPow(Scatter2DPtr h, const double n, const double k = 1.0) {
828 if (n == 0 || n == 1) {
829 cout << "Error: Do not take the 0th or 1st power of a Scatter2D,"
830 " use scale instead." << endl;
831 return;
832 }
833 vector<YODA::Point2D> points;
834 vector<YODA::Point2D> pIn = h->points();
835 // The error pre-factor is k^(1/n) / n by Taylors formula.
836 double eFac = pow(k,1./n)/n;
837 for (const auto& b : pIn) {
838 const double yVal = pow(k * b.y(), n);
839 if (std::isnan(yVal)) {
840 points.push_back(YODA::Point2D(b.x(), 0., b.xErrMinus(),
841 b.xErrPlus(), 0, 0 ));
842 }
843 else {
844 double yemin = abs(eFac * pow(yVal,1./(n - 1.))) * b.yErrMinus();
845 if (std::isnan(yemin)) yemin = b.yErrMinus();
846 double yemax = abs(eFac * pow(yVal,1./(n - 1.))) * b.yErrPlus();
847 if (std::isnan(yemax)) yemax = b.yErrPlus();
848 points.push_back(YODA::Point2D(b.x(), yVal, b.xErrMinus(),
849 b.xErrPlus(), yemin, yemax ));
850 }
851 }
852 h->reset();
853 for (int i = 0, N = points.size(); i < N; ++i) h->addPoint(points[i]);
854 }
855
856
861 template<typename T>
862 static pair<double, double> sampleVariance(T func) {
863 // First we calculate the mean (two pass calculation).
864 double avg = 0.;
865 for (int i = 0; i < BOOT_BINS; ++i) avg += func(i);
866 avg /= double(BOOT_BINS);
867 // Then we find the variance.
868 double var = 0.;
869 for (int i = 0; i < BOOT_BINS; ++i) var += pow(func(i) - avg, 2.);
870 var /= (double(BOOT_BINS) - 1);
871 return pair<double, double>(var, var);
872 }
873
874
879 template<typename T>
880 static pair<double, double> sampleEnvelope(T func) {
881 // First we calculate the mean.
882 double avg = 0.;
883 for (int i = 0; i < BOOT_BINS; ++i) avg += func(i);
884 avg /= double(BOOT_BINS);
885 double yMax = avg;
886 double yMin = avg;
887 // Then we find the envelope using the mean as initial value.
888 for (int i = 0; i < BOOT_BINS; ++i) {
889 double yVal = func(i);
890 if (yMin > yVal) yMin = yVal;
891 else if (yMax < yVal) yMax = yVal;
892 }
893 return pair<double, double>(fabs(avg - yMin), fabs(yMax - avg));
894 }
895
896
898 template<typename T>
899 const pair<double, double> sampleError(T func) const {
900 if (errorMethod == VARIANCE) return sampleVariance(func);
901 else if (errorMethod == ENVELOPE) return sampleEnvelope(func);
902 else
903 cout << "Error: Error method not found!" << endl;
904 return pair<double, double>(0.,0.);
905 }
906
907
909 void cnTwoInt(Scatter2DPtr h, ECorrPtr e2) const {
910 const vector<CorBin>& bins = e2->getBins();
911 const vector<double>& binx = e2->getBinX();
912 // Assert bin size.
913 if (binx.size() - 1 != bins.size()){
914 cout << "cnTwoInt: Bin size (x,y) differs!" << endl;
915 return;
916 }
917 vector<const CorBinBase*> binPtrs;
918 // The mean value of the cumulant.
919 auto cn = [&](const int i) { return binPtrs[i]->mean(); };
920 // Error calculation.
921 vector<pair<double,double>> yErr;
922 for (int j = 0, N = bins.size(); j < N; ++j) {
923 binPtrs = bins[j].getBinPtrs();
924 yErr.push_back(sampleError(cn));
925 }
926 binPtrs = e2->getBinPtrs();
927 fillScatter(h, binx, cn, yErr);
928 }
929
930
932 void vnTwoInt(Scatter2DPtr h, ECorrPtr e2) const {
933 cnTwoInt(h, e2);
934 nthPow(h, 0.5);
935 }
936
937
941 void corrPlot(Scatter2DPtr h, ECorrPtr e) const {
942 cnTwoInt(h, e);
943 }
944
945
946
947
948 // TODO Use full path for lookup, change to single AU in output, rename.
949 void rawHookIn(YODA::AnalysisObjectPtr yao) final {
950 // Fill the corresponding ECorrelator.
951 for (auto ec : eCorrPtrs) if(ec->fillFromProfile(yao, name())) break;;
952 }
953
958 void rawHookOut(const vector<MultiplexAOPtr>& raos, size_t iW) final {
959 // Loop over the correlators and extract the numbers.
960 for (auto ec : eCorrPtrs) {
961 const vector<CorBin>& corBins = ec->getBins();
962 const vector<double>& binx = ec->getBinX();
963 auto ref = ec->getReference();
964 auto refBins = ref.getBinPtrs<CorSingleBin>();
965 // Assert bin size.
966 if (binx.size() - 1 != corBins.size()){
967 cout << "corrPlot: Bin size (x,y) differs!" << endl;
968 return;
969 }
970 // Loop over the booked histograms using their names.
971 for (int i = 0, N = ec->profs.size(); i < N; ++i) {
972 for (auto rao : raos) {
973 if (rao->path() != "/"+name()+"/TMP/"+ec->profs[i]) continue;
974 // Get a pointer to the active profile.
975 rao.get()->setActiveWeightIdx(iW);
976 YODA::Profile1DPtr pPtr = dynamic_pointer_cast<YODA::Profile1D>(rao.get()->activeAO());
977 // New bins.
978 vector<YODA::Dbn2D> profBins;
979 // Add reference flow in the underflow bin.
980 pPtr->bin(0).set( YODA::Dbn2D(refBins[i]->numEntries(), refBins[i]->sumW(),
981 refBins[i]->sumW2(), 0., 0., refBins[i]->sumWX(), 0., 0.) );
982 for (size_t j = 0, N = binx.size() - 1; j < N; ++j) {
983 vector<const CorSingleBin*> binPtrs = corBins[j].getBinPtrs<CorSingleBin>();
984 // Construct bin of the profiled quantities and put the
985 // ECorrelator into the raw histogram. We have no information
986 // (and no desire to add it) of sumWX of the profile, so really
987 // we should use a Dbn1D - but that does not work for Profile1D's.
988 pPtr->bin(j).set( YODA::Dbn2D(binPtrs[i]->numEntries(), binPtrs[i]->sumW(),
989 binPtrs[i]->sumW2(), 0., 0., binPtrs[i]->sumWX(), 0, 0) );
990 }
991 }
992 }
993 }
994 }
995
997 void cnFourInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4) const {
998 const auto& e2bins = e2->getBins();
999 const auto& e4bins = e4->getBins();
1000 const vector<double>& binx = e2->getBinX();
1001 if (binx.size() - 1 != e2bins.size()){
1002 cout << "cnFourInt: Bin size (x,y) differs!" << endl;
1003 return;
1004 }
1005 if (binx != e4->getBinX()) {
1006 cout << "Error in cnFourInt: Correlator x-binning differs!" << endl;
1007 return;
1008 }
1009 vector<const CorBinBase*> e2binPtrs;
1010 vector<const CorBinBase*> e4binPtrs;
1011 auto cn = [&] (const int i) {
1012 double e22 = e2binPtrs[i]->mean() * e2binPtrs[i]->mean();
1013 return e4binPtrs[i]->mean() - 2. * e22;
1014 };
1015 // Error calculation.
1016 vector<pair<double, double> > yErr;
1017 for (int j = 0, N = e2bins.size(); j < N; ++j) {
1018 e2binPtrs = e2bins[j].getBinPtrs();
1019 e4binPtrs = e4bins[j].getBinPtrs();
1020 yErr.push_back(sampleError(cn));
1021 }
1022 // Put the bin ptrs back in place.
1023 e2binPtrs = e2->getBinPtrs();
1024 e4binPtrs = e4->getBinPtrs();
1025 fillScatter(h, binx, cn, yErr);
1026 }
1027
1028
1030 void vnFourInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4) const {
1031 cnFourInt(h, e2, e4);
1032 nthPow(h, 0.25, -1.0);
1033 }
1034
1035
1038 ECorrPtr e6) const {
1039 const auto& e2bins = e2->getBins();
1040 const auto& e4bins = e4->getBins();
1041 const auto& e6bins = e6->getBins();
1042 const auto& binx = e2->getBinX();
1043 if (binx.size() - 1 != e2bins.size()){
1044 cout << "cnSixInt: Bin size (x,y) differs!" << endl;
1045 return;
1046 }
1047 if (binx != e4->getBinX() || binx != e6->getBinX()) {
1048 cout << "Error in cnSixInt: Correlator x-binning differs!" << endl;
1049 return;
1050 }
1051 vector<const CorBinBase*> e2binPtrs;
1052 vector<const CorBinBase*> e4binPtrs;
1053 vector<const CorBinBase*> e6binPtrs;
1054 auto cn = [&] (const int i) {
1055 const double e2mean = e2binPtrs[i]->mean();
1056 const double e4mean = e4binPtrs[i]->mean();
1057 const double e6mean = e6binPtrs[i]->mean();
1058 const double e23 = pow(e2mean, 3.0);
1059 return e6mean - 9.*e2mean*e4mean + 12.*e23;
1060 };
1061 // Error calculation.
1062 vector<pair<double, double> > yErr;
1063 for (int j = 0, N = e2bins.size(); j < N; ++j) {
1064 e2binPtrs = e2bins[j].getBinPtrs();
1065 e4binPtrs = e4bins[j].getBinPtrs();
1066 e6binPtrs = e6bins[j].getBinPtrs();
1067 yErr.push_back(sampleError(cn));
1068 }
1069 // Put the bin ptrs back in place.
1070 e2binPtrs = e2->getBinPtrs();
1071 e4binPtrs = e4->getBinPtrs();
1072 e6binPtrs = e6->getBinPtrs();
1073 fillScatter(h, binx, cn, yErr);
1074 }
1075
1076
1079 ECorrPtr e6) const {
1080 cnSixInt(h, e2, e4, e6);
1081 nthPow(h, 1./6., 0.25);
1082 }
1083
1084
1087 ECorrPtr e6, ECorrPtr e8) const {
1088 const auto& e2bins = e2->getBins();
1089 const auto& e4bins = e4->getBins();
1090 const auto& e6bins = e6->getBins();
1091 const auto& e8bins = e8->getBins();
1092 const vector<double>& binx = e2->getBinX();
1093 if (binx.size() - 1 != e2bins.size()){
1094 cout << "cnEightInt: Bin size (x,y) differs!" << endl;
1095 return;
1096 }
1097 if (binx != e4->getBinX() || binx != e6->getBinX() || binx != e8->getBinX()) {
1098 cout << "Error in cnEightInt: Correlator x-binning differs!" << endl;
1099 return;
1100 }
1101 vector<const CorBinBase*> e2binPtrs;
1102 vector<const CorBinBase*> e4binPtrs;
1103 vector<const CorBinBase*> e6binPtrs;
1104 vector<const CorBinBase*> e8binPtrs;
1105 auto cn = [&] (const int i) {
1106 const double e2mean = e2binPtrs[i]->mean();
1107 const double e4mean = e4binPtrs[i]->mean();
1108 const double e6mean = e6binPtrs[i]->mean();
1109 const double e8mean = e8binPtrs[i]->mean();
1110 const double e22 = sqr(e2mean);
1111 const double e24 = sqr(e22);
1112 const double e42 = sqr(e4mean);
1113 return e8mean - 16. * e6mean * e2mean - 18. * e42 + 144. * e4mean*e22 - 144. * e24;
1114 };
1115 // Error calculation.
1116 vector<pair<double, double> > yErr;
1117 for (int j = 0, N = e2bins.size(); j < N; ++j) {
1118 e2binPtrs = e2bins[j].getBinPtrs();
1119 e4binPtrs = e4bins[j].getBinPtrs();
1120 e6binPtrs = e6bins[j].getBinPtrs();
1121 e8binPtrs = e8bins[j].getBinPtrs();
1122 yErr.push_back(sampleError(cn));
1123 }
1124 // Put the bin ptrs back in place.
1125 e2binPtrs = e2->getBinPtrs();
1126 e4binPtrs = e4->getBinPtrs();
1127 e6binPtrs = e6->getBinPtrs();
1128 e8binPtrs = e8->getBinPtrs();
1129 fillScatter(h, binx, cn, yErr);
1130 }
1131
1132
1135 cnEightInt(h, e2, e4, e6, e8);
1136 nthPow(h, 1./8., -1./33.);
1137 }
1138
1139
1141 void vnTwoDiff(Scatter2DPtr h, ECorrPtr e2Dif) const {
1142 const auto& e2bins = e2Dif->getBins();
1143 const auto& ref = e2Dif->getReference();
1144 const auto& binx = e2Dif->getBinX();
1145 if (binx.size() -1 != e2bins.size()) {
1146 cout << "vnTwoDif: Bin size (x,y) differs!" << endl;
1147 return;
1148 }
1149 vector<const CorBinBase*> e2binPtrs;
1150 vector<const CorBinBase*> refPtrs;
1151 auto vn = [&] (const int i) {
1152 // Test reference flow.
1153 if (ref.mean() <= 0) return 0.;
1154 return e2binPtrs[i]->mean() / sqrt(ref.mean());
1155 };
1156 // We need here a separate error function, as we don't iterate over the reference flow.
1157 auto vnerr = [&] (const int i) {
1158 // Test reference flow.
1159 if (refPtrs[i]->mean() <=0) return 0.;
1160 return e2binPtrs[i]->mean() / sqrt(refPtrs[i]->mean());
1161 };
1162 // Error calculation.
1163 vector<pair<double,double>> yErr;
1164 refPtrs = ref.getBinPtrs();
1165 for (int j = 0, N = e2bins.size(); j < N; ++j) {
1166 e2binPtrs = e2bins[j].getBinPtrs();
1167 yErr.push_back(sampleError(vnerr));
1168 }
1169 // Put the e2binPtrs back in place.
1170 e2binPtrs = e2Dif->getBinPtrs();
1171 fillScatter(h, binx, vn);
1172 }
1173
1174
1176 void vnFourDiff(Scatter2DPtr h, ECorrPtr e2Dif, ECorrPtr e4Dif) const {
1177 const auto& e2bins = e2Dif->getBins();
1178 const auto& e4bins = e4Dif->getBins();
1179 const auto& ref2 = e2Dif->getReference();
1180 const auto& ref4 = e4Dif->getReference();
1181 const auto& binx = e2Dif->getBinX();
1182 if (binx.size() - 1 != e2bins.size()){
1183 cout << "vnFourDif: Bin size (x,y) differs!" << endl;
1184 return;
1185 }
1186 if (binx != e4Dif->getBinX()) {
1187 cout << "Error in vnFourDif: Correlator x-binning differs!" << endl;
1188 return;
1189 }
1190 vector<const CorBinBase*> e2binPtrs;
1191 vector<const CorBinBase*> e4binPtrs;
1192 vector<const CorBinBase*> ref2Ptrs;
1193 vector<const CorBinBase*> ref4Ptrs;
1194 double denom = 2 * ref2.mean() * ref2.mean() - ref4.mean();
1195 auto vn = [&] (const int i) {
1196 // Test denominator.
1197 if (denom <= 0 ) return 0.;
1198 return ((2 * ref2.mean() * e2bins[i].mean() - e4bins[i].mean()) / pow(denom, 0.75));
1199 };
1200 // We need here a separate error function, as we don't iterate over the reference flow.
1201 auto vnerr = [&] (const int i) {
1202 double denom2 = 2 * ref2Ptrs[i]->mean() * ref2Ptrs[i]->mean() -
1203 ref4Ptrs[i]->mean();
1204 // Test denominator.
1205 if (denom2 <= 0) return 0.;
1206 return ((2 * ref2Ptrs[i]->mean() * e2binPtrs[i]->mean() - e4binPtrs[i]->mean()) / pow(denom2, 0.75));
1207 };
1208 // Error calculation.
1209 vector<pair<double,double> > yErr;
1210 ref2Ptrs = ref2.getBinPtrs();
1211 ref4Ptrs = ref4.getBinPtrs();
1212 for (int j = 0, N = e2bins.size(); j < N; ++j) {
1213 e2binPtrs = e2bins[j].getBinPtrs();
1214 e4binPtrs = e4bins[j].getBinPtrs();
1215 yErr.push_back(sampleError(vnerr));
1216 }
1217 // Put the binPtrs back in place.
1218 e2binPtrs = e2Dif->getBinPtrs();
1219 e4binPtrs = e4Dif->getBinPtrs();
1220 fillScatter(h, binx, vn, yErr);
1221 }
1222
1223 };
1224
1225
1226}
1227
1228#endif
This is the base class of all analysis classes in Rivet.
Definition Analysis.hh:67
double sumW2() const
Get the sum of squared event weights seen (via the analysis handler).
double sumW() const
Get the sum of event weights seen (via the analysis handler).
virtual std::string name() const
Get the name of the analysis.
Definition Analysis.hh:142
Projection for calculating correlators for flow measurements.
Definition Correlators.hh:37
const vector< pair< double, double > > pTBinnedCorrelatorsGap(const Correlators &other, vector< int > n1, vector< int > n2, bool overflow=false) const
pT differential correlators of n1 harmonic, for number n1.size()
const complex< double > getP(int n, int p, double pT=0.) const
Return a P-vector.
Definition Correlators.hh:167
CmpState compare(const Projection &p) const
Compare to other projection, testing harmonics, pT bins and underlying final state similarity.
Definition Correlators.hh:148
void project(const Event &e)
Loop over array and calculates Q and P vectors if needed.
void fillCorrelators(const Particle &p, const double &weight)
Calculate correlators from one particle.
const vector< pair< double, double > > pTBinnedCorrelators(vector< int > n, bool overflow=false) const
pT differential correlator of n harmonic, for number of powers n.size()
const pair< double, double > intCorrelator(vector< int > n) const
Integrated correlator of n harmonic, with the number of powers being the size of n....
const complex< double > getQ(int n, int p) const
Return a Q-vector.
Definition Correlators.hh:160
Correlators(const ParticleFinder &fsp, int nMaxIn=2, int pMaxIn=0, vector< double > pTbinEdgesIn={})
const pair< double, double > intCorrelatorGap(const Correlators &other, vector< int > n1, vector< int > n2) const
Integrated correlator of n1 harmonic, for number of powers n1.size()
static vector< int > hVec(int n, int m)
Construct a harmonic vectors from n harmonics and m number of particles.
Definition Correlators.hh:109
static pair< int, int > getMaxValues(vector< vector< int > > &hList)
Return the maximal values for n, p to be used in the constructor of Correlators(xxx,...
Definition Correlators.hh:124
A helper class to calculate all event averages of correlators.
Definition Correlators.hh:403
const vector< int > & getH1() const
Get a copy of the h1 harmonic vector.
Definition Correlators.hh:509
void setProfs(vector< string > prIn)
Set the prIn list of profile histograms associated with the internal bins.
Definition Correlators.hh:532
const vector< double > & getBinX() const
Get a copy of the bin x-values.
Definition Correlators.hh:504
void setReference(CorBin refIn)
Replace reference flow bin with another, e.g. calculated in another phase space or with other pid.
Definition Correlators.hh:519
const vector< CorBin > & getBins() const
Get the bin contents.
Definition Correlators.hh:492
void fill(const double &obs, const Correlators &c, double weight=1.0)
Fill the appropriate bin given an input (per event) observable, e.g. centrality.
Definition Correlators.hh:433
ECorrelator(const vector< int > &h1In, const vector< int > &h2In, const vector< double > &binIn)
Constructor for gapped correlator.
Definition Correlators.hh:428
const vector< int > & getH2() const
Get a copy of the h2 harmonic vector.
Definition Correlators.hh:514
bool fillFromProfile(YODA::AnalysisObjectPtr yao, string name)
Fill bins with content from preloaded histograms.
Definition Correlators.hh:537
vector< const CorBinBase * > getBinPtrs() const
Return the bins as pointers to the base class.
Definition Correlators.hh:497
void fill(const double &obs, const Correlators &c1, const Correlators &c2, double weight=1.0)
Fill the appropriate bin given an input (per event) observable, e.g. centrality, with a rapidity gap ...
Definition Correlators.hh:442
CorBin getReference() const
Extract the reference flow from a differential event averaged correlator.
Definition Correlators.hh:524
void fill(const Correlators &c1, const Correlators &c2, const double weight=1.0)
Fill bins with the appropriate correlator, and a rapidity gap between two Correlators.
Definition Correlators.hh:473
void fill(const Correlators &c, const double weight=1.0)
Fill the bins with the appropriate correlator.
Definition Correlators.hh:456
ECorrelator(const vector< int > &h, const vector< double > &binIn)
Constructor. Takes as argument the desired harmonic and number of correlated particles as a generic f...
Definition Correlators.hh:420
Tools for flow analyses.
Definition Correlators.hh:224
void vnTwoDiff(Scatter2DPtr h, ECorrPtr e2Dif) const
Two-particle differential vn.
Definition Correlators.hh:1141
void cnFourInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4) const
Four particle integrated cn.
Definition Correlators.hh:997
void vnSixInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4, ECorrPtr e6) const
Six-particle integrated vn.
Definition Correlators.hh:1078
ECorrPtr bookECorrelatorGap(const string &name, const vector< int > &h, const YODA::Estimate1D &hIn)
Definition Correlators.hh:665
void vnFourDiff(Scatter2DPtr h, ECorrPtr e2Dif, ECorrPtr e4Dif) const
Four-particle differential vn.
Definition Correlators.hh:1176
ECorrPtr bookECorrelatorGap(const string &name, const YODA::Estimate1D &hIn)
Templated version of gapped correlator booking which takes N desired harmonic and M number of particl...
Definition Correlators.hh:695
ECorrPtr bookECorrelator(const string name, const vector< int > &h, const vector< double > &binIn)
Book an ECorrelator in the same way as a histogram.
Definition Correlators.hh:621
const pair< double, double > sampleError(T func) const
Selection method for which sample error to use, given in the constructor.
Definition Correlators.hh:899
const pair< int, int > getMaxValues() const
Get the correct max N and max P for the set of booked correlators.
Definition Correlators.hh:590
ECorrPtr bookECorrelator(const string name, const vector< int > &h, const YODA::Estimate1D &hIn)
Book an ECorrelator in the same way as a histogram.
Definition Correlators.hh:611
static void fillScatter(Scatter2DPtr h, const vector< double > &binx, const T &func)
Helper method for turning correlators into Scatter2Ds.
Definition Correlators.hh:727
void cnSixInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4, ECorrPtr e6) const
Six-particle integrated cn.
Definition Correlators.hh:1037
ECorrPtr bookECorrelator(const string &name, const vector< int > &h1, const vector< int > &h2, const YODA::Estimate1D &hIn)
Book a gapped ECorrelator with two harmonic vectors.
Definition Correlators.hh:652
ECorrPtr bookECorrelator(const string &name, const YODA::Estimate1D &hIn)
Templated version of correlator booking which takes N desired harmonic and M number of particles.
Definition Correlators.hh:686
CumulantAnalysis(const string &n)
Constructor.
Definition Correlators.hh:715
void vnFourInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4) const
Four-particle integrated vn.
Definition Correlators.hh:1030
static void nthPow(Scatter2DPtr h, const double n, const double k=1.0)
Take the n th power of all points in h, and put the result back in the same Scatter2D.
Definition Correlators.hh:827
void corrPlot(Scatter2DPtr h, ECorrPtr e) const
Put an event-averaged correlator into a Scatter2D.
Definition Correlators.hh:941
void cnTwoInt(Scatter2DPtr h, ECorrPtr e2) const
Two-particle integrated cn.
Definition Correlators.hh:909
static pair< double, double > sampleVariance(T func)
Calculate the bootstrapped sample variance.
Definition Correlators.hh:862
ECorrPtr bookECorrelator(const string &name, const vector< double > &binIn)
Templated version of correlator booking which takes N desired harmonic and M number of particles,...
Definition Correlators.hh:677
void rawHookOut(const vector< MultiplexAOPtr > &raos, size_t iW) final
Transform RAW ECorrelator Profiles to have content before writing them. Overloaded method from Analys...
Definition Correlators.hh:958
static pair< double, double > sampleEnvelope(T func)
Calculate the bootstrapped sample envelope.
Definition Correlators.hh:880
static void nthPow(Scatter2DPtr hOut, const Scatter2DPtr hIn, const double &n, const double &k=1.0)
Take the n th power of all points in hIn and put the result in hOut.
Definition Correlators.hh:790
void vnTwoInt(Scatter2DPtr h, ECorrPtr e2) const
Two particle integrated vn.
Definition Correlators.hh:932
void vnEightInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4, ECorrPtr e6, ECorrPtr e8) const
Eight-particle integrated vn.
Definition Correlators.hh:1134
void fillScatter(Scatter2DPtr h, const vector< double > &binx, const F func, vector< pair< double, double > > &yErr) const
Helper method for turning correlators into Scatter2Ds with error estimates.
Definition Correlators.hh:758
shared_ptr< ECorrelator > ECorrPtr
Typedef of shared pointer to ECorrelator.
Definition Correlators.hh:607
void cnEightInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4, ECorrPtr e6, ECorrPtr e8) const
Eight-particle integrated cn.
Definition Correlators.hh:1086
ECorrPtr bookECorrelator(const string name, const vector< int > &h1, const vector< int > &h2, const vector< double > &binIn)
Book a gapped ECorrelator with two harmonic vectors.
Definition Correlators.hh:636
Representation of a HepMC event, and enabler of Projection caching.
Definition Event.hh:22
Definition RivetYODA.hh:1330
Base class for projections which return subsets of an event's particles.
Definition ParticleFinder.hh:11
Particle representation, either from a HepMC::GenEvent or reconstructed.
Definition Particle.hh:45
Base class for all Rivet projections.
Definition Projection.hh:29
Cmp< Projection > mkPCmp(const Projection &otherparent, const std::string &pname) const
CounterPtr & book(CounterPtr &, const std::string &name)
Book a counter.
double pT(const Vector3 &a, const Vector3 &b)
Calculate transverse momentum of pair of 3-vectors.
Definition Vector3.hh:639
Definition MC_CENT_PPB_Projections.hh:10
std::enable_if< std::is_arithmetic< NUM1 >::value &&std::is_arithmetic< NUM2 >::value, int >::type binIndex(NUM1 val, std::initializer_list< NUM2 > binedges, bool allow_overflow=false)
Return the bin index of the given value, val, given a vector of bin edges.
Definition MathUtils.hh:455
static constexpr double DBL_NAN
Convenient const for getting the double NaN value.
Definition Utils.hh:38
std::enable_if< std::is_arithmetic< NUM >::value, double >::type mean(const vector< NUM > &sample)
Definition MathUtils.hh:497
std::enable_if< std::is_arithmetic< NUM >::value, NUM >::type sqr(NUM a)
Named number-type squaring operation.
Definition MathUtils.hh:219