|
|
Rivet 4.0.0
|
2#ifndef RIVET_Correlators_HH
3#define RIVET_Correlators_HH
22#include "Rivet/Analysis.hh"
23#include "Rivet/Projection.hh"
24#include "Rivet/Projections/ParticleFinder.hh"
25#include "YODA/Scatter.h"
54 int pMaxIn = 0, vector<double> pTbinEdgesIn = {});
58 int pMaxIn, const YODA::Estimate1D hIn);
61 using Projection::operator =;
78 bool overflow = false) const;
89 vector<int> n1, vector<int> n2) const;
102 const vector<pair<double,double>>
109 static vector<int> hVec( int n, int m) {
111 cout << "Harmonic Vector: Number of particles must be an even number." << endl;
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);
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) {
132 if (tmpN > nMax) nMax = tmpN;
133 if (tmpP > pMax) pMax = tmpP;
135 return make_pair(nMax,pMax);
150 if (nMax != other->nMax) return CmpState::NEQ;
151 if (pMax != other->pMax) return CmpState::NEQ;
152 if (pTbinEdges != other->pTbinEdges) return CmpState::NEQ;
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];
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];
181 const complex<double> recCorr( int order, vector<int> n,
182 vector<int> p, bool useP, double pT = 0.) const;
188 const complex<double> twoPartCorr( int n1, int n2, int p1 = 1,
189 int p2 = 1, double pT = 0., bool useP = false) const;
195 const complex<double> _ZERO = {0., 0.};
196 const double _TINY = 1e-10;
199 typedef vector< vector<complex<double>> > Vec2D;
203 map<double, Vec2D> pVec;
209 vector<double> pTbinEdges;
230 static const int BOOT_BINS = 9;
247 virtual ~CorBinBase() {};
249 virtual void fill( const pair<double, double>& cor, const double& weight = 1.0) = 0;
250 virtual double mean() const = 0;
259 class CorSingleBin : public CorBinBase {
264 : _sumWX(0.), _sumW(0.), _sumW2(0.), _numEntries(0.)
273 void fill( const pair<double, double>& cor, const double& weight = 1.0) {
275 if (cor.second < 1e-10) return;
278 _sumWX += cor.first * weight;
279 _sumW += weight * cor.second;
280 _sumW2 += weight * weight * cor.second * cor.second;
285 double mean() const {
286 if (_sumW < 1e-10) return 0;
287 return _sumWX / _sumW;
291 double sumW() const {
296 double sumW2() const {
301 double sumWX() const {
306 double numEntries() const {
311 void addContent( double ne, double sw, double sw2, double swx) {
321 double _sumWX, _sumW, _sumW2, _numEntries;
330 class CorBin : public CorBinBase {
338 CorBin() : binIndex(0), nBins(BOOT_BINS) {
339 for( size_t i = 0; i < nBins; ++i) bins.push_back(CorSingleBin());
346 void fill( const pair<double, double>& cor, const double& weight = 1.0) {
348 if (cor.second < 1e-10) return;
356 double mean() const {
359 for ( auto b : bins) {
360 if (b.sumW() < 1e-10) continue;
368 const vector<CorSingleBin>& getBins() const {
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;});
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;});
390 vector<CorSingleBin> bins;
421 : h1(h), h2({}), binX(binIn), binContent(binIn.size() - 1), reference()
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()
434 int index = getBinIndex(obs);
435 if (index < 0) return;
444 cout << "Trying to fill gapped correlator, but harmonics behind the gap (h2) are not given!" << endl;
447 int index = getBinIndex(obs);
448 if (index < 0) return;
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);
475 cout << "Trying to fill gapped correlator, but harmonics behind "
476 "the gap (h2) are not given!" << endl;
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);
498 vector<const CorBinBase*> ret(binContent.size());
499 transform(binContent.begin(), binContent.end(), ret.begin(), []( const CorBin& b) {return &b;});
525 if (reference.mean() < 1e-10)
526 cout << "Warning: ECorrelator, reference bin is zero." << endl;
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(),
549 const YODA::Dbn2D& uBin = pPtr->bin(0);
550 refs[i]->addContent(uBin.numEntries(), uBin.sumW(), uBin.sumW2(),
561 int getBinIndex( const double& obs) const {
564 if (obs >= binX.back()) return -1;
566 if (obs < binX[0]) return -1;
568 for ( int i = 0, N = binX.size() - 1; i < N; ++i, ++index)
569 if (obs >= binX[i] && obs < binX[i + 1]) break;
579 vector<CorBin> binContent;
584 vector <string> profs;
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);
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>();
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());
623 vector<string> eCorrProfs;
624 for ( int i = 0; i < BOOT_BINS; ++i) {
626 book(tmp, "TMP/"+ name+ "-"+to_string(i),binIn);
627 eCorrProfs.push_back( name+ "-"+to_string(i));
629 ecPtr->setProfs(eCorrProfs);
630 eCorrPtrs.push_back(ecPtr);
637 const vector<int>& h2, const vector<double>& binIn) {
639 vector<string> eCorrProfs;
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));
645 ecPtr->setProfs(eCorrProfs);
646 eCorrPtrs.push_back(ecPtr);
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());
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());
676 template< unsigned int N, unsigned int M>
685 template< unsigned int N, unsigned int M>
694 template< unsigned int N, unsigned int 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());
706 list<ECorrPtr> eCorrPtrs;
716 : Analysis(n), errorMethod(VARIANCE)
728 vector<YODA::Point2D> points;
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]);
736 xMid = h->points()[i].x();
737 xeMin = h->points()[i].xErrMinus();
738 xeMax = h->points()[i].xErrPlus();
740 double yVal = func(i);
741 if (std::isnan(yVal)) yVal = 0.;
743 points.push_back(YODA::Point2D(xMid, yVal, xeMin, xeMax, yErr, yErr));
747 for ( int i = 0, N = points.size(); i < N; ++i) h->addPoint(points[i]);
759 vector<pair<double, double> >& yErr) const {
760 vector<YODA::Point2D> points;
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]);
768 xMid = h->points()[i].x();
769 xeMin = h->points()[i].xErrMinus();
770 xeMax = h->points()[i].xErrPlus();
772 const double yVal = func(i);
773 if (std::isnan(yVal))
774 points.push_back(YODA::Point2D(xMid, 0., xeMin, xeMax,0., 0.));
776 points.push_back(YODA::Point2D(xMid, yVal, xeMin, xeMax,
777 yErr[i].first, yErr[i].second));
781 for ( int i = 0, N = points.size(); i < N; ++i) {
782 h->addPoint(points[i]);
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;
796 if (hIn->points().size() != hOut->points().size()) {
797 cout << "nthRoot: Scatterplots: " << hIn->name() << " and " <<
798 hOut->name() << " not initialized with same length." << endl;
801 vector<YODA::Point2D> points;
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 ));
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 ));
819 for ( int i = 0, N = points.size(); i < N; ++i)
820 hOut->addPoint(points[i]);
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;
833 vector<YODA::Point2D> points;
834 vector<YODA::Point2D> pIn = h->points();
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 ));
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 ));
853 for ( int i = 0, N = points.size(); i < N; ++i) h->addPoint(points[i]);
865 for ( int i = 0; i < BOOT_BINS; ++i) avg += func(i);
866 avg /= double(BOOT_BINS);
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);
883 for ( int i = 0; i < BOOT_BINS; ++i) avg += func(i);
884 avg /= double(BOOT_BINS);
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;
893 return pair<double, double>(fabs(avg - yMin), fabs(yMax - avg));
903 cout << "Error: Error method not found!" << endl;
904 return pair<double, double>(0.,0.);
910 const vector<CorBin>& bins = e2->getBins();
911 const vector<double>& binx = e2->getBinX();
913 if (binx.size() - 1 != bins.size()){
914 cout << "cnTwoInt: Bin size (x,y) differs!" << endl;
917 vector<const CorBinBase*> binPtrs;
919 auto cn = [&]( const int i) { return binPtrs[i]->mean(); };
921 vector<pair<double,double>> yErr;
922 for ( int j = 0, N = bins.size(); j < N; ++j) {
923 binPtrs = bins[j].getBinPtrs();
926 binPtrs = e2->getBinPtrs();
949 void rawHookIn(YODA::AnalysisObjectPtr yao) final {
951 for ( auto ec : eCorrPtrs) if(ec->fillFromProfile(yao, name())) break;;
958 void rawHookOut( const vector<MultiplexAOPtr>& raos, size_t iW) final {
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>();
966 if (binx.size() - 1 != corBins.size()){
967 cout << "corrPlot: Bin size (x,y) differs!" << endl;
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;
975 rao.get()->setActiveWeightIdx(iW);
976 YODA::Profile1DPtr pPtr = dynamic_pointer_cast<YODA::Profile1D>(rao.get()->activeAO());
978 vector<YODA::Dbn2D> profBins;
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>();
988 pPtr->bin(j).set( YODA::Dbn2D(binPtrs[i]->numEntries(), binPtrs[i]->sumW(),
989 binPtrs[i]->sumW2(), 0., 0., binPtrs[i]->sumWX(), 0, 0) );
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;
1005 if (binx != e4->getBinX()) {
1006 cout << "Error in cnFourInt: Correlator x-binning differs!" << endl;
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;
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();
1023 e2binPtrs = e2->getBinPtrs();
1024 e4binPtrs = e4->getBinPtrs();
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;
1047 if (binx != e4->getBinX() || binx != e6->getBinX()) {
1048 cout << "Error in cnSixInt: Correlator x-binning differs!" << endl;
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;
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();
1070 e2binPtrs = e2->getBinPtrs();
1071 e4binPtrs = e4->getBinPtrs();
1072 e6binPtrs = e6->getBinPtrs();
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;
1097 if (binx != e4->getBinX() || binx != e6->getBinX() || binx != e8->getBinX()) {
1098 cout << "Error in cnEightInt: Correlator x-binning differs!" << endl;
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;
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();
1125 e2binPtrs = e2->getBinPtrs();
1126 e4binPtrs = e4->getBinPtrs();
1127 e6binPtrs = e6->getBinPtrs();
1128 e8binPtrs = e8->getBinPtrs();
1136 nthPow(h, 1./8., -1./33.);
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;
1149 vector<const CorBinBase*> e2binPtrs;
1150 vector<const CorBinBase*> refPtrs;
1151 auto vn = [&] ( const int i) {
1153 if (ref.mean() <= 0) return 0.;
1154 return e2binPtrs[i]->mean() / sqrt(ref.mean());
1157 auto vnerr = [&] ( const int i) {
1159 if (refPtrs[i]-> mean() <=0) return 0.;
1160 return e2binPtrs[i]->mean() / sqrt(refPtrs[i]-> mean());
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();
1170 e2binPtrs = e2Dif->getBinPtrs();
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;
1186 if (binx != e4Dif->getBinX()) {
1187 cout << "Error in vnFourDif: Correlator x-binning differs!" << endl;
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) {
1197 if (denom <= 0 ) return 0.;
1198 return ((2 * ref2.mean() * e2bins[i].mean() - e4bins[i].mean()) / pow(denom, 0.75));
1201 auto vnerr = [&] ( const int i) {
1202 double denom2 = 2 * ref2Ptrs[i]->mean() * ref2Ptrs[i]->mean() -
1203 ref4Ptrs[i]->mean();
1205 if (denom2 <= 0) return 0.;
1206 return ((2 * ref2Ptrs[i]-> mean() * e2binPtrs[i]->mean() - e4binPtrs[i]->mean()) / pow(denom2, 0.75));
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();
1218 e2binPtrs = e2Dif->getBinPtrs();
1219 e4binPtrs = e4Dif->getBinPtrs();
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
|