diff --git a/PWGCF/GenericFramework/Core/FlowPtContainer.cxx b/PWGCF/GenericFramework/Core/FlowPtContainer.cxx index ad5f5880d53..eb142b9ca31 100644 --- a/PWGCF/GenericFramework/Core/FlowPtContainer.cxx +++ b/PWGCF/GenericFramework/Core/FlowPtContainer.cxx @@ -14,9 +14,9 @@ /// \author Emil Gorm Nielsen, NBI, emil.gorm.nielsen@cern.ch #include "FlowPtContainer.h" + #include #include -#include FlowPtContainer::FlowPtContainer() : fCMTermList(0), fCorrList(0), @@ -26,27 +26,24 @@ FlowPtContainer::FlowPtContainer() : fCMTermList(0), fCumulantList(0), fCentralMomentList(0), mpar(0), + nSubevents(0), fillCounter(0), fEventWeight(EventWeight::UnityWeight), fUseCentralMoments(true), fUseGap(false), sumP(), - insub1(), - insub2(), + insub(), corrNum(), - corrNum1(), - corrNum2(), + corrNumSub(), corrDen(), - corrDen1(), - corrDen2(), + corrDenSub(), cmVal(), - cmVal1(), - cmVal2(), + cmValSub(), cmDen(), - cmDen1(), - cmDen2(), + cmDenSub(), arr(), - warr() {} + warr(), + subevents() {} FlowPtContainer::~FlowPtContainer() { delete fCMTermList; @@ -66,27 +63,24 @@ FlowPtContainer::FlowPtContainer(const char* name) : TNamed(name, name), fCumulantList(0), fCentralMomentList(0), mpar(0), + nSubevents(0), fillCounter(0), fEventWeight(EventWeight::UnityWeight), fUseCentralMoments(true), fUseGap(false), sumP(), - insub1(), - insub2(), + insub(), corrNum(), - corrNum1(), - corrNum2(), + corrNumSub(), corrDen(), - corrDen1(), - corrDen2(), + corrDenSub(), cmVal(), - cmVal1(), - cmVal2(), + cmValSub(), cmDen(), - cmDen1(), - cmDen2(), + cmDenSub(), arr(), - warr() {} + warr(), + subevents() {} FlowPtContainer::FlowPtContainer(const char* name, const char* title) : TNamed(name, title), fCMTermList(0), fCorrList(0), @@ -96,27 +90,24 @@ FlowPtContainer::FlowPtContainer(const char* name, const char* title) : TNamed(n fCumulantList(0), fCentralMomentList(0), mpar(0), + nSubevents(0), fillCounter(0), fEventWeight(EventWeight::UnityWeight), fUseCentralMoments(true), fUseGap(false), sumP(), - insub1(), - insub2(), + insub(), corrNum(), - corrNum1(), - corrNum2(), + corrNumSub(), corrDen(), - corrDen1(), - corrDen2(), + corrDenSub(), cmVal(), - cmVal1(), - cmVal2(), + cmValSub(), cmDen(), - cmDen1(), - cmDen2(), + cmDenSub(), arr(), - warr() {} + warr(), + subevents() {} void FlowPtContainer::initialise(const o2::framework::AxisSpec axis, const int& m, const GFWCorrConfigs& configs, const int& nsub) { arr.resize(3 * 3 * 3 * 3); @@ -391,8 +382,13 @@ void FlowPtContainer::initialise(int nbinsx, double xlow, double xhigh, const in } LOGF(info, "Container %s initialized with m = %i\n", this->GetName(), mpar); }; -void FlowPtContainer::initialiseSubevent(const o2::framework::AxisSpec axis, const int& m, const int& nsub) +void FlowPtContainer::initialiseSubevent(const o2::framework::AxisSpec axis, const int& m, const int& nsubev, const int& nsub) { + if (nsubev < 1) { + LOGF(fatal, "Need at least one subevent"); + return; + } + nSubevents = nsubev; if (!mpar) mpar = m; std::vector multiBins = axis.binEdges; @@ -408,16 +404,28 @@ void FlowPtContainer::initialiseSubevent(const o2::framework::AxisSpec axis, con delete fSubList; fSubList = new TList(); fSubList->SetOwner(kTRUE); - for (int subEv = 0; subEv < 2; ++subEv) { - for (int m = 0; m < mpar; ++m) { - fSubList->Add(new BootstrapProfile(Form("mpt_sub%i_%ipar", subEv + 1, m + 1), this->GetTitle(), nMultiBins, &multiBins[0])); - } - } - for (int m = 2; m <= mpar; ++m) { - for (int k = 0; k < m - 1; ++k) { - fSubList->Add(new BootstrapProfile(Form("mpt_%isub1_%isub2_%ipar", m - k - 1, k + 1, m), this->GetTitle(), nMultiBins, &multiBins[0])); + + // Get all possible subevent combinations given m particles and nsubev subevents - also considering not using all m particles, e.g. all lower orders + std::vector current; + getSubevents(m, nsubev + 1, current, subevents); + // remove unused "extra" subevent + for (auto& subevent : subevents) + subevent.pop_back(); + subevents.erase(subevents.begin(), subevents.begin() + 1); + + std::vector histnames; + for (auto subevent : subevents) { + std::string outstr = "ptptsub"; + int index = 0; + for (auto particles : subevent) { + outstr += "_" + std::to_string(particles) + + "sub" + std::to_string(index + 1); + ++index; } + histnames.push_back(outstr); } + for (const auto& name : histnames) + fSubList->Add(new BootstrapProfile(name.c_str(), this->GetTitle(), nMultiBins, &multiBins[0])); if (fSubCMList) delete fSubCMList; @@ -453,8 +461,13 @@ void FlowPtContainer::initialiseSubevent(const o2::framework::AxisSpec axis, con } LOGF(info, "Container %s initialized Subevents and %i subsamples", this->GetName(), nsub); } -void FlowPtContainer::initialiseSubevent(int nbinsx, double* xbins, const int& m, const int& nsub) +void FlowPtContainer::initialiseSubevent(int nbinsx, double* xbins, const int& m, const int& nsubev, const int& nsub) { + if (nsubev < 1) { + LOGF(fatal, "Need at least one subevent"); + return; + } + nSubevents = nsubev; if (!mpar) mpar = m; @@ -462,16 +475,28 @@ void FlowPtContainer::initialiseSubevent(int nbinsx, double* xbins, const int& m delete fSubList; fSubList = new TList(); fSubList->SetOwner(kTRUE); - for (int subEv = 0; subEv < 2; ++subEv) { - for (int m = 0; m < mpar; ++m) { - fSubList->Add(new BootstrapProfile(Form("mpt_sub%i_%ipar", subEv + 1, m + 1), this->GetTitle(), nbinsx, xbins)); - } - } - for (int m = 2; m <= mpar; ++m) { - for (int k = 0; k < m - 1; ++k) { - fSubList->Add(new BootstrapProfile(Form("mpt_%isub1_%isub2_%ipar", m - k - 1, k + 1, m), this->GetTitle(), nbinsx, xbins)); + + // Get all possible subevent combinations given m particles and nsubev subevents - also considering not using all m particles, e.g. all lower orders + std::vector current; + getSubevents(m, nsubev + 1, current, subevents); + // remove unused "extra" subevent + for (auto& subevent : subevents) + subevent.pop_back(); + subevents.erase(subevents.begin(), subevents.begin() + 1); + + std::vector histnames; + for (auto subevent : subevents) { + std::string outstr = "ptptsub"; + int index = 0; + for (auto particles : subevent) { + outstr += "_" + std::to_string(particles) + + "sub" + std::to_string(index + 1); + ++index; } + histnames.push_back(outstr); } + for (const auto& name : histnames) + fSubList->Add(new BootstrapProfile(name.c_str(), this->GetTitle(), nbinsx, xbins)); if (fSubCMList) delete fSubCMList; @@ -507,24 +532,41 @@ void FlowPtContainer::initialiseSubevent(int nbinsx, double* xbins, const int& m } LOGF(info, "Container %s initialized Subevents and %i subsamples", this->GetName(), nsub); } -void FlowPtContainer::initialiseSubevent(int nbinsx, double xlow, double xhigh, const int& m, const int& nsub) +void FlowPtContainer::initialiseSubevent(int nbinsx, double xlow, double xhigh, const int& m, const int& nsubev, const int& nsub) { + if (nsubev < 1) { + LOGF(fatal, "Need at least one subevent"); + return; + } + nSubevents = nsubev; if (!mpar) mpar = m; if (fSubList) delete fSubList; fSubList = new TList(); fSubList->SetOwner(kTRUE); - for (int subEv = 0; subEv < 2; ++subEv) { - for (int m = 0; m < mpar; ++m) { - fSubList->Add(new BootstrapProfile(Form("mpt_sub%i_%ipar", subEv + 1, m + 1), this->GetTitle(), nbinsx, xlow, xhigh)); - } - } - for (int m = 2; m <= mpar; ++m) { - for (int k = 0; k < m - 1; ++k) { - fSubList->Add(new BootstrapProfile(Form("mpt_%isub1_%isub2_%ipar", m - k - 1, k + 1, m), this->GetTitle(), nbinsx, xlow, xhigh)); + + // Get all possible subevent combinations given m particles and nsubev subevents - also considering not using all m particles, e.g. all lower orders + std::vector current; + getSubevents(mpar, nsubev + 1, current, subevents); + // remove unused "extra" subevent + for (auto& subevent : subevents) + subevent.pop_back(); + subevents.erase(subevents.begin(), subevents.begin() + 1); + + std::vector histnames; + for (auto subevent : subevents) { + std::string outstr = "ptptsub"; + int index = 0; + for (auto particles : subevent) { + outstr += "_" + std::to_string(particles) + + "sub" + std::to_string(index + 1); + ++index; } + histnames.push_back(outstr); } + for (const auto& name : histnames) + fSubList->Add(new BootstrapProfile(name.c_str(), this->GetTitle(), nbinsx, xlow, xhigh)); if (fSubCMList) delete fSubCMList; @@ -566,19 +608,11 @@ void FlowPtContainer::fill(const double& w, const double& pt) } return; } -void FlowPtContainer::fillSub1(const double& w, const double& pt) -{ - for (size_t i = 0; i < insub1.size(); ++i) { - insub1[i] += std::pow(w, i % (mpar + 1)) * std::pow(pt, i / (mpar + 1)); - } - return; -} -void FlowPtContainer::fillSub2(const double& w, const double& pt) +void FlowPtContainer::fillSub(const double& w, const double& pt, int subIndex) { - for (size_t i = 0; i < insub2.size(); ++i) { - insub2[i] += std::pow(w, i % (mpar + 1)) * std::pow(pt, i / (mpar + 1)); + for (size_t i = 0; i < insub[subIndex].size(); ++i) { + insub[subIndex][i] += std::pow(w, i % (mpar + 1)) * std::pow(pt, i / (mpar + 1)); } - return; } void FlowPtContainer::calculateCorrelations() { @@ -609,49 +643,31 @@ void FlowPtContainer::calculateCorrelations() } void FlowPtContainer::calculateSubeventCorrelations() { - corrNum1.clear(); - corrNum1.resize(mpar + 1, 0); - corrNum1[0] = 1.0; - corrDen1.clear(); - corrDen1.resize(mpar + 1, 0); - corrDen1[0] = 1.0; - corrNum2.clear(); - corrNum2.resize(mpar + 1, 0); - corrNum2[0] = 1.0; - corrDen2.clear(); - corrDen2.resize(mpar + 1, 0); - corrDen2[0] = 1.0; + corrNumSub.clear(); + corrNumSub.resize(nSubevents, std::vector(mpar + 1, 0)); + for (auto& corrnum : corrNumSub) + corrnum[0] = 1.0; + corrDenSub.resize(nSubevents, std::vector(mpar + 1, 0)); + for (auto& corrden : corrDenSub) + corrden[0] = 1.0; - double sumNum1 = 0; - double sumDenum1 = 0; - std::vector valNum1; - std::vector valDenum1; - double sumNum2 = 0; - double sumDenum2 = 0; - std::vector valNum2; - std::vector valDenum2; - - for (int m(1); m <= mpar; ++m) { - for (int k(1); k <= m; ++k) { - // correlations in subevent 1 - valNum1.push_back(SignArray[k - 1] * corrNum1[m - k] * (FactorialArray[m - 1] / FactorialArray[m - k]) * insub1[getVectorIndex(k, k)]); - valDenum1.push_back(SignArray[k - 1] * corrDen1[m - k] * (FactorialArray[m - 1] / FactorialArray[m - k]) * insub1[getVectorIndex(k, 0)]); - // correlations in subevent 2 - valNum2.push_back(SignArray[k - 1] * corrNum2[m - k] * (FactorialArray[m - 1] / FactorialArray[m - k]) * insub2[getVectorIndex(k, k)]); - valDenum2.push_back(SignArray[k - 1] * corrDen2[m - k] * (FactorialArray[m - 1] / FactorialArray[m - k]) * insub2[getVectorIndex(k, 0)]); + for (int subIndex = 0; subIndex < nSubevents; ++subIndex) { + double sumNum = 0.0; + double sumDenum = 0.0; + std::vector valNum; + std::vector valDenum; + for (int m(1); m <= mpar; ++m) { + for (int k(1); k <= m; ++k) { + valNum.push_back(SignArray[k - 1] * corrNumSub[subIndex][m - k] * (FactorialArray[m - 1] / FactorialArray[m - k]) * insub[subIndex][getVectorIndex(k, k)]); + valDenum.push_back(SignArray[k - 1] * corrDenSub[subIndex][m - k] * (FactorialArray[m - 1] / FactorialArray[m - k]) * insub[subIndex][getVectorIndex(k, 0)]); + } + sumNum = orderedAddition(valNum); + sumDenum = orderedAddition(valDenum); + valNum.clear(); + valDenum.clear(); + corrNumSub[subIndex][m] = sumNum; + corrDenSub[subIndex][m] = sumDenum; } - sumNum1 = orderedAddition(valNum1); - sumDenum1 = orderedAddition(valDenum1); - sumNum2 = orderedAddition(valNum2); - sumDenum2 = orderedAddition(valDenum2); - valNum1.clear(); - valDenum1.clear(); - valNum2.clear(); - valDenum2.clear(); - corrNum1[m] = sumNum1; - corrDen1[m] = sumDenum1; - corrNum2[m] = sumNum2; - corrDen2[m] = sumDenum2; } return; } @@ -666,21 +682,25 @@ void FlowPtContainer::fillPtProfiles(const double& centmult, const double& rn) } void FlowPtContainer::fillSubeventPtProfiles(const double& centmult, const double& rn) { - // Fill the correlations within subevents, requires that the CalculateSubeventCorrelations have been called with the correct input vectors right before - for (int m = 1; m <= mpar; ++m) { - if (corrDen1[m] != 0) - dynamic_cast(fSubList->At(m - 1))->FillProfile(centmult, corrNum1[m] / corrDen1[m], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : corrDen1[m], rn); - // subevent 2 profiles offset by mpar positions - if (corrDen2[m] != 0) - dynamic_cast(fSubList->At(mpar + m - 1))->FillProfile(centmult, corrNum2[m] / corrDen2[m], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : corrDen2[m], rn); - } - - // Fill the cross-subevent correlations - for (int m = 2; m <= mpar; ++m) { - for (int k = 0; k < m - 1; ++k) { - if (corrDen1[m - k - 1] != 0 && corrDen2[k + 1] != 0) - dynamic_cast(fSubList->FindObject(Form("mpt_%isub1_%isub2_%ipar", m - k - 1, k + 1, m)))->FillProfile(centmult, corrNum1[m - k - 1] / corrDen1[m - k - 1] * corrNum2[k + 1] / corrDen2[k + 1], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : corrDen1[m - k - 1] * corrDen2[k + 1], rn); + int histCounter = 0; + for (const auto& subevent : subevents) { + double val = 1.0; + double dn = 1.0; + int subIndex = 0; + bool valid = true; + for (auto m : subevent) { + if (corrDenSub[subIndex][m] == 0) { + valid = false; + break; + } else { + val *= corrNumSub[subIndex][m] / corrDenSub[subIndex][m]; + dn *= corrDenSub[subIndex][m]; + } + ++subIndex; } + if (valid) + dynamic_cast(fSubList->At(histCounter))->FillProfile(centmult, val, (fEventWeight == EventWeight::UnityWeight) ? 1.0 : dn, rn); + ++histCounter; } return; } @@ -879,91 +899,93 @@ void FlowPtContainer::fillCMSubeventProfiles(const double& centmult, const doubl // do I need to add an extra return statement here to match fillCMProfiles? if (mpar < 1) return; + if (nSubevents < 2) + return; int indOffset = 0; for (int im = 1; im <= mpar; im++) { indOffset += im; } // 0th order correlation - cmDen1.push_back(1.); - cmVal1.push_back(1.); - cmDen2.push_back(1.); - cmVal2.push_back(1.); + cmDenSub[0].push_back(1.); + cmValSub[0].push_back(1.); + cmDenSub[nSubevents - 1].push_back(1.); + cmValSub[nSubevents - 1].push_back(1.); - cmDen1.push_back(insub1[getVectorIndex(1, 0)]); - cmDen1.push_back(insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] - insub1[getVectorIndex(2, 0)]); - cmDen1.push_back(insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] - 3 * insub1[getVectorIndex(2, 0)] * insub1[getVectorIndex(1, 0)] + 2 * insub1[getVectorIndex(3, 0)]); - cmDen1.push_back(insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] - 6 * insub1[getVectorIndex(2, 0)] * insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] + 8 * insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(3, 0)] + 3 * insub1[getVectorIndex(2, 0)] * insub1[getVectorIndex(2, 0)] - 6 * insub1[getVectorIndex(4, 0)]); + cmDenSub[0].push_back(insub[0][getVectorIndex(1, 0)]); + cmDenSub[0].push_back(insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] - insub[0][getVectorIndex(2, 0)]); + cmDenSub[0].push_back(insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] - 3 * insub[0][getVectorIndex(2, 0)] * insub[0][getVectorIndex(1, 0)] + 2 * insub[0][getVectorIndex(3, 0)]); + cmDenSub[0].push_back(insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] - 6 * insub[0][getVectorIndex(2, 0)] * insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] + 8 * insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(3, 0)] + 3 * insub[0][getVectorIndex(2, 0)] * insub[0][getVectorIndex(2, 0)] - 6 * insub[0][getVectorIndex(4, 0)]); - cmDen2.push_back(insub2[getVectorIndex(1, 0)]); - cmDen2.push_back(insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] - insub2[getVectorIndex(2, 0)]); - cmDen2.push_back(insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] - 3 * insub2[getVectorIndex(2, 0)] * insub2[getVectorIndex(1, 0)] + 2 * insub2[getVectorIndex(3, 0)]); - cmDen2.push_back(insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] - 6 * insub2[getVectorIndex(2, 0)] * insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] + 8 * insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(3, 0)] + 3 * insub2[getVectorIndex(2, 0)] * insub2[getVectorIndex(2, 0)] - 6 * insub2[getVectorIndex(4, 0)]); + cmDenSub[nSubevents - 1].push_back(insub[nSubevents - 1][getVectorIndex(1, 0)]); + cmDenSub[nSubevents - 1].push_back(insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] - insub[nSubevents - 1][getVectorIndex(2, 0)]); + cmDenSub[nSubevents - 1].push_back(insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] - 3 * insub[nSubevents - 1][getVectorIndex(2, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] + 2 * insub[nSubevents - 1][getVectorIndex(3, 0)]); + cmDenSub[nSubevents - 1].push_back(insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] - 6 * insub[nSubevents - 1][getVectorIndex(2, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] + 8 * insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(3, 0)] + 3 * insub[nSubevents - 1][getVectorIndex(2, 0)] * insub[nSubevents - 1][getVectorIndex(2, 0)] - 6 * insub[nSubevents - 1][getVectorIndex(4, 0)]); - if (cmDen1[1] != 0) { - cmVal1.push_back(insub1[getVectorIndex(1, 1)] / cmDen1[1]); - dynamic_cast(fSubCMList->At(0))->FillProfile(centmult, cmVal1[1], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen1[1], rn); + if (cmDenSub[0][1] != 0) { + cmValSub[0].push_back(insub[0][getVectorIndex(1, 1)] / cmDenSub[0][1]); + dynamic_cast(fSubCMList->At(0))->FillProfile(centmult, cmValSub[0][1], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[0][1], rn); } - if (cmDen2[1] != 0) { - cmVal2.push_back(insub2[getVectorIndex(1, 1)] / cmDen2[1]); - dynamic_cast(fSubCMList->At(indOffset + 0))->FillProfile(centmult, cmVal2[1], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen2[1], rn); + if (cmDenSub[nSubevents - 1][1] != 0) { + cmValSub[nSubevents - 1].push_back(insub[nSubevents - 1][getVectorIndex(1, 1)] / cmDenSub[nSubevents - 1][1]); + dynamic_cast(fSubCMList->At(indOffset + 0))->FillProfile(centmult, cmValSub[nSubevents - 1][1], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[nSubevents - 1][1], rn); } if (mpar < 2) return; - if (insub1[getVectorIndex(2, 0)] != 0 && cmDen1[2] != 0) { - cmVal1.push_back(1 / cmDen1[2] * (insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] - insub1[getVectorIndex(2, 2)])); - dynamic_cast(fSubCMList->At(1))->FillProfile(centmult, cmVal1[2], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen1[2], rn); - cmVal1.push_back(1 / cmDen1[2] * (insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 1)] - insub1[getVectorIndex(2, 1)])); - dynamic_cast(fSubCMList->At(2))->FillProfile(centmult, cmVal1[3], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen1[2], rn); - } - if (insub2[getVectorIndex(2, 0)] != 0 && cmDen2[2] != 0) { - cmVal2.push_back(1 / cmDen2[2] * (insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] - insub2[getVectorIndex(2, 2)])); - dynamic_cast(fSubCMList->At(indOffset + 1))->FillProfile(centmult, cmVal2[2], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen2[2], rn); - cmVal2.push_back(1 / cmDen2[2] * (insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 1)] - insub2[getVectorIndex(2, 1)])); - dynamic_cast(fSubCMList->At(indOffset + 2))->FillProfile(centmult, cmVal2[3], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen2[2], rn); + if (insub[0][getVectorIndex(2, 0)] != 0 && cmDenSub[0][2] != 0) { + cmValSub[0].push_back(1 / cmDenSub[0][2] * (insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] - insub[0][getVectorIndex(2, 2)])); + dynamic_cast(fSubCMList->At(1))->FillProfile(centmult, cmValSub[0][2], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[0][2], rn); + cmValSub[0].push_back(1 / cmDenSub[0][2] * (insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 1)] - insub[0][getVectorIndex(2, 1)])); + dynamic_cast(fSubCMList->At(2))->FillProfile(centmult, cmValSub[0][3], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[0][2], rn); + } + if (insub[nSubevents - 1][getVectorIndex(2, 0)] != 0 && cmDenSub[nSubevents - 1][2] != 0) { + cmValSub[nSubevents - 1].push_back(1 / cmDenSub[nSubevents - 1][2] * (insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] - insub[nSubevents - 1][getVectorIndex(2, 2)])); + dynamic_cast(fSubCMList->At(indOffset + 1))->FillProfile(centmult, cmValSub[nSubevents - 1][2], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[nSubevents - 1][2], rn); + cmValSub[nSubevents - 1].push_back(1 / cmDenSub[nSubevents - 1][2] * (insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 1)] - insub[nSubevents - 1][getVectorIndex(2, 1)])); + dynamic_cast(fSubCMList->At(indOffset + 2))->FillProfile(centmult, cmValSub[nSubevents - 1][3], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[nSubevents - 1][2], rn); } if (mpar < 3) return; - if (insub1[getVectorIndex(3, 0)] != 0 && cmDen1[3] != 0) { - cmVal1.push_back(1 / cmDen1[3] * (insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] - 3 * insub1[getVectorIndex(2, 2)] * insub1[getVectorIndex(1, 1)] + 2 * insub1[getVectorIndex(3, 3)])); - dynamic_cast(fSubCMList->At(3))->FillProfile(centmult, cmVal1[4], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen1[3], rn); - cmVal1.push_back(1 / cmDen1[3] * (insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 0)] - 2 * insub1[getVectorIndex(2, 1)] * insub1[getVectorIndex(1, 1)] + 2 * insub1[getVectorIndex(3, 2)] - insub1[getVectorIndex(2, 2)] * insub1[getVectorIndex(1, 0)])); - dynamic_cast(fSubCMList->At(4))->FillProfile(centmult, cmVal1[5], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen1[3], rn); - cmVal1.push_back(1 / cmDen1[3] * (insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] - 2 * insub1[getVectorIndex(2, 1)] * insub1[getVectorIndex(1, 0)] + 2 * insub1[getVectorIndex(3, 1)] - insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(2, 0)])); - dynamic_cast(fSubCMList->At(5))->FillProfile(centmult, cmVal1[6], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen1[3], rn); - } - if (insub2[getVectorIndex(3, 0)] != 0 && cmDen2[3] != 0) { - cmVal2.push_back(1 / cmDen2[3] * (insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] - 3 * insub2[getVectorIndex(2, 2)] * insub2[getVectorIndex(1, 1)] + 2 * insub2[getVectorIndex(3, 3)])); - dynamic_cast(fSubCMList->At(indOffset + 3))->FillProfile(centmult, cmVal2[4], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen2[3], rn); - cmVal2.push_back(1 / cmDen2[3] * (insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 0)] - 2 * insub2[getVectorIndex(2, 1)] * insub2[getVectorIndex(1, 1)] + 2 * insub2[getVectorIndex(3, 2)] - insub2[getVectorIndex(2, 2)] * insub2[getVectorIndex(1, 0)])); - dynamic_cast(fSubCMList->At(indOffset + 4))->FillProfile(centmult, cmVal2[5], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen2[3], rn); - cmVal2.push_back(1 / cmDen2[3] * (insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] - 2 * insub2[getVectorIndex(2, 1)] * insub2[getVectorIndex(1, 0)] + 2 * insub2[getVectorIndex(3, 1)] - insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(2, 0)])); - dynamic_cast(fSubCMList->At(indOffset + 5))->FillProfile(centmult, cmVal2[6], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen2[3], rn); + if (insub[0][getVectorIndex(3, 0)] != 0 && cmDenSub[0][3] != 0) { + cmValSub[0].push_back(1 / cmDenSub[0][3] * (insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] - 3 * insub[0][getVectorIndex(2, 2)] * insub[0][getVectorIndex(1, 1)] + 2 * insub[0][getVectorIndex(3, 3)])); + dynamic_cast(fSubCMList->At(3))->FillProfile(centmult, cmValSub[0][4], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[0][3], rn); + cmValSub[0].push_back(1 / cmDenSub[0][3] * (insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 0)] - 2 * insub[0][getVectorIndex(2, 1)] * insub[0][getVectorIndex(1, 1)] + 2 * insub[0][getVectorIndex(3, 2)] - insub[0][getVectorIndex(2, 2)] * insub[0][getVectorIndex(1, 0)])); + dynamic_cast(fSubCMList->At(4))->FillProfile(centmult, cmValSub[0][5], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[0][3], rn); + cmValSub[0].push_back(1 / cmDenSub[0][3] * (insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] - 2 * insub[0][getVectorIndex(2, 1)] * insub[0][getVectorIndex(1, 0)] + 2 * insub[0][getVectorIndex(3, 1)] - insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(2, 0)])); + dynamic_cast(fSubCMList->At(5))->FillProfile(centmult, cmValSub[0][6], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[0][3], rn); + } + if (insub[nSubevents - 1][getVectorIndex(3, 0)] != 0 && cmDenSub[nSubevents - 1][3] != 0) { + cmValSub[nSubevents - 1].push_back(1 / cmDenSub[nSubevents - 1][3] * (insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] - 3 * insub[nSubevents - 1][getVectorIndex(2, 2)] * insub[nSubevents - 1][getVectorIndex(1, 1)] + 2 * insub[nSubevents - 1][getVectorIndex(3, 3)])); + dynamic_cast(fSubCMList->At(indOffset + 3))->FillProfile(centmult, cmValSub[nSubevents - 1][4], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[nSubevents - 1][3], rn); + cmValSub[nSubevents - 1].push_back(1 / cmDenSub[nSubevents - 1][3] * (insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 0)] - 2 * insub[nSubevents - 1][getVectorIndex(2, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] + 2 * insub[nSubevents - 1][getVectorIndex(3, 2)] - insub[nSubevents - 1][getVectorIndex(2, 2)] * insub[nSubevents - 1][getVectorIndex(1, 0)])); + dynamic_cast(fSubCMList->At(indOffset + 4))->FillProfile(centmult, cmValSub[nSubevents - 1][5], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[nSubevents - 1][3], rn); + cmValSub[nSubevents - 1].push_back(1 / cmDenSub[nSubevents - 1][3] * (insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] - 2 * insub[nSubevents - 1][getVectorIndex(2, 1)] * insub[nSubevents - 1][getVectorIndex(1, 0)] + 2 * insub[nSubevents - 1][getVectorIndex(3, 1)] - insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(2, 0)])); + dynamic_cast(fSubCMList->At(indOffset + 5))->FillProfile(centmult, cmValSub[nSubevents - 1][6], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[nSubevents - 1][3], rn); } if (mpar < 4) return; - if (insub1[getVectorIndex(4, 0)] != 0 && cmDen1[4] != 0) { - cmVal1.push_back(1 / cmDen1[4] * (insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] - 6 * insub1[getVectorIndex(2, 2)] * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] + 3 * insub1[getVectorIndex(2, 2)] * insub1[getVectorIndex(2, 2)] + 8 * insub1[getVectorIndex(3, 3)] * insub1[getVectorIndex(1, 1)] - 6 * insub1[getVectorIndex(4, 4)])); - dynamic_cast(fSubCMList->At(6))->FillProfile(centmult, cmVal1[7], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen1[4], rn); - cmVal1.push_back(1 / cmDen1[4] * (insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 0)] - 3 * insub1[getVectorIndex(2, 2)] * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 0)] - 3 * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(2, 1)] + 3 * insub1[getVectorIndex(2, 2)] * insub1[getVectorIndex(2, 1)] + 6 * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(3, 2)] - 6 * insub1[getVectorIndex(4, 3)])); - dynamic_cast(fSubCMList->At(7))->FillProfile(centmult, cmVal1[8], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen1[4], rn); - cmVal1.push_back(1 / cmDen1[4] * (insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] - insub1[getVectorIndex(2, 2)] * insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] - insub1[getVectorIndex(2, 0)] * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 1)] + insub1[getVectorIndex(2, 0)] * insub1[getVectorIndex(2, 2)] - 4 * insub1[getVectorIndex(2, 1)] * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 0)] + 4 * insub1[getVectorIndex(3, 2)] * insub1[getVectorIndex(1, 0)] + 4 * insub1[getVectorIndex(3, 1)] * insub1[getVectorIndex(1, 1)] + 2 * insub1[getVectorIndex(2, 1)] * insub1[getVectorIndex(2, 1)] - 6 * insub1[getVectorIndex(4, 2)])); - dynamic_cast(fSubCMList->At(8))->FillProfile(centmult, cmVal1[9], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen1[4], rn); - cmVal1.push_back(1 / cmDen1[4] * (insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] - 3 * insub1[getVectorIndex(2, 1)] * insub1[getVectorIndex(1, 0)] * insub1[getVectorIndex(1, 0)] - 3 * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(2, 0)] * insub1[getVectorIndex(1, 0)] + 3 * insub1[getVectorIndex(2, 1)] * insub1[getVectorIndex(2, 0)] + 2 * insub1[getVectorIndex(1, 1)] * insub1[getVectorIndex(3, 0)] + 6 * insub1[getVectorIndex(3, 1)] * insub1[getVectorIndex(1, 0)] - 6 * insub1[getVectorIndex(4, 1)])); - dynamic_cast(fSubCMList->At(9))->FillProfile(centmult, cmVal1[10], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen1[4], rn); - } - if (insub2[getVectorIndex(4, 0)] != 0 && cmDen2[4] != 0) { - cmVal2.push_back(1 / cmDen2[4] * (insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] - 6 * insub2[getVectorIndex(2, 2)] * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] + 3 * insub2[getVectorIndex(2, 2)] * insub2[getVectorIndex(2, 2)] + 8 * insub2[getVectorIndex(3, 3)] * insub2[getVectorIndex(1, 1)] - 6 * insub2[getVectorIndex(4, 4)])); - dynamic_cast(fSubCMList->At(indOffset + 6))->FillProfile(centmult, cmVal2[7], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen2[4], rn); - cmVal2.push_back(1 / cmDen2[4] * (insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 0)] - 3 * insub2[getVectorIndex(2, 2)] * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 0)] - 3 * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(2, 1)] + 3 * insub2[getVectorIndex(2, 2)] * insub2[getVectorIndex(2, 1)] + 6 * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(3, 2)] - 6 * insub2[getVectorIndex(4, 3)])); - dynamic_cast(fSubCMList->At(indOffset + 7))->FillProfile(centmult, cmVal2[8], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen2[4], rn); - cmVal2.push_back(1 / cmDen2[4] * (insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] - insub2[getVectorIndex(2, 2)] * insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] - insub2[getVectorIndex(2, 0)] * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 1)] + insub2[getVectorIndex(2, 0)] * insub2[getVectorIndex(2, 2)] - 4 * insub2[getVectorIndex(2, 1)] * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 0)] + 4 * insub2[getVectorIndex(3, 2)] * insub2[getVectorIndex(1, 0)] + 4 * insub2[getVectorIndex(3, 1)] * insub2[getVectorIndex(1, 1)] + 2 * insub2[getVectorIndex(2, 1)] * insub2[getVectorIndex(2, 1)] - 6 * insub2[getVectorIndex(4, 2)])); - dynamic_cast(fSubCMList->At(indOffset + 8))->FillProfile(centmult, cmVal2[9], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen2[4], rn); - cmVal2.push_back(1 / cmDen2[4] * (insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] - 3 * insub2[getVectorIndex(2, 1)] * insub2[getVectorIndex(1, 0)] * insub2[getVectorIndex(1, 0)] - 3 * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(2, 0)] * insub2[getVectorIndex(1, 0)] + 3 * insub2[getVectorIndex(2, 1)] * insub2[getVectorIndex(2, 0)] + 2 * insub2[getVectorIndex(1, 1)] * insub2[getVectorIndex(3, 0)] + 6 * insub2[getVectorIndex(3, 1)] * insub2[getVectorIndex(1, 0)] - 6 * insub2[getVectorIndex(4, 1)])); - dynamic_cast(fSubCMList->At(indOffset + 9))->FillProfile(centmult, cmVal2[10], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen2[4], rn); + if (insub[0][getVectorIndex(4, 0)] != 0 && cmDenSub[0][4] != 0) { + cmValSub[0].push_back(1 / cmDenSub[0][4] * (insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] - 6 * insub[0][getVectorIndex(2, 2)] * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] + 3 * insub[0][getVectorIndex(2, 2)] * insub[0][getVectorIndex(2, 2)] + 8 * insub[0][getVectorIndex(3, 3)] * insub[0][getVectorIndex(1, 1)] - 6 * insub[0][getVectorIndex(4, 4)])); + dynamic_cast(fSubCMList->At(6))->FillProfile(centmult, cmValSub[0][7], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[0][4], rn); + cmValSub[0].push_back(1 / cmDenSub[0][4] * (insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 0)] - 3 * insub[0][getVectorIndex(2, 2)] * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 0)] - 3 * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(2, 1)] + 3 * insub[0][getVectorIndex(2, 2)] * insub[0][getVectorIndex(2, 1)] + 6 * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(3, 2)] - 6 * insub[0][getVectorIndex(4, 3)])); + dynamic_cast(fSubCMList->At(7))->FillProfile(centmult, cmValSub[0][8], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[0][4], rn); + cmValSub[0].push_back(1 / cmDenSub[0][4] * (insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] - insub[0][getVectorIndex(2, 2)] * insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] - insub[0][getVectorIndex(2, 0)] * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 1)] + insub[0][getVectorIndex(2, 0)] * insub[0][getVectorIndex(2, 2)] - 4 * insub[0][getVectorIndex(2, 1)] * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 0)] + 4 * insub[0][getVectorIndex(3, 2)] * insub[0][getVectorIndex(1, 0)] + 4 * insub[0][getVectorIndex(3, 1)] * insub[0][getVectorIndex(1, 1)] + 2 * insub[0][getVectorIndex(2, 1)] * insub[0][getVectorIndex(2, 1)] - 6 * insub[0][getVectorIndex(4, 2)])); + dynamic_cast(fSubCMList->At(8))->FillProfile(centmult, cmValSub[0][9], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[0][4], rn); + cmValSub[0].push_back(1 / cmDenSub[0][4] * (insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] - 3 * insub[0][getVectorIndex(2, 1)] * insub[0][getVectorIndex(1, 0)] * insub[0][getVectorIndex(1, 0)] - 3 * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(2, 0)] * insub[0][getVectorIndex(1, 0)] + 3 * insub[0][getVectorIndex(2, 1)] * insub[0][getVectorIndex(2, 0)] + 2 * insub[0][getVectorIndex(1, 1)] * insub[0][getVectorIndex(3, 0)] + 6 * insub[0][getVectorIndex(3, 1)] * insub[0][getVectorIndex(1, 0)] - 6 * insub[0][getVectorIndex(4, 1)])); + dynamic_cast(fSubCMList->At(9))->FillProfile(centmult, cmValSub[0][10], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[0][4], rn); + } + if (insub[nSubevents - 1][getVectorIndex(4, 0)] != 0 && cmDenSub[nSubevents - 1][4] != 0) { + cmValSub[nSubevents - 1].push_back(1 / cmDenSub[nSubevents - 1][4] * (insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] - 6 * insub[nSubevents - 1][getVectorIndex(2, 2)] * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] + 3 * insub[nSubevents - 1][getVectorIndex(2, 2)] * insub[nSubevents - 1][getVectorIndex(2, 2)] + 8 * insub[nSubevents - 1][getVectorIndex(3, 3)] * insub[nSubevents - 1][getVectorIndex(1, 1)] - 6 * insub[nSubevents - 1][getVectorIndex(4, 4)])); + dynamic_cast(fSubCMList->At(indOffset + 6))->FillProfile(centmult, cmValSub[nSubevents - 1][7], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[nSubevents - 1][4], rn); + cmValSub[nSubevents - 1].push_back(1 / cmDenSub[nSubevents - 1][4] * (insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 0)] - 3 * insub[nSubevents - 1][getVectorIndex(2, 2)] * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 0)] - 3 * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(2, 1)] + 3 * insub[nSubevents - 1][getVectorIndex(2, 2)] * insub[nSubevents - 1][getVectorIndex(2, 1)] + 6 * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(3, 2)] - 6 * insub[nSubevents - 1][getVectorIndex(4, 3)])); + dynamic_cast(fSubCMList->At(indOffset + 7))->FillProfile(centmult, cmValSub[nSubevents - 1][8], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[nSubevents - 1][4], rn); + cmValSub[nSubevents - 1].push_back(1 / cmDenSub[nSubevents - 1][4] * (insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] - insub[nSubevents - 1][getVectorIndex(2, 2)] * insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] - insub[nSubevents - 1][getVectorIndex(2, 0)] * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] + insub[nSubevents - 1][getVectorIndex(2, 0)] * insub[nSubevents - 1][getVectorIndex(2, 2)] - 4 * insub[nSubevents - 1][getVectorIndex(2, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 0)] + 4 * insub[nSubevents - 1][getVectorIndex(3, 2)] * insub[nSubevents - 1][getVectorIndex(1, 0)] + 4 * insub[nSubevents - 1][getVectorIndex(3, 1)] * insub[nSubevents - 1][getVectorIndex(1, 1)] + 2 * insub[nSubevents - 1][getVectorIndex(2, 1)] * insub[nSubevents - 1][getVectorIndex(2, 1)] - 6 * insub[nSubevents - 1][getVectorIndex(4, 2)])); + dynamic_cast(fSubCMList->At(indOffset + 8))->FillProfile(centmult, cmValSub[nSubevents - 1][9], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[nSubevents - 1][4], rn); + cmValSub[nSubevents - 1].push_back(1 / cmDenSub[nSubevents - 1][4] * (insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] - 3 * insub[nSubevents - 1][getVectorIndex(2, 1)] * insub[nSubevents - 1][getVectorIndex(1, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] - 3 * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(2, 0)] * insub[nSubevents - 1][getVectorIndex(1, 0)] + 3 * insub[nSubevents - 1][getVectorIndex(2, 1)] * insub[nSubevents - 1][getVectorIndex(2, 0)] + 2 * insub[nSubevents - 1][getVectorIndex(1, 1)] * insub[nSubevents - 1][getVectorIndex(3, 0)] + 6 * insub[nSubevents - 1][getVectorIndex(3, 1)] * insub[nSubevents - 1][getVectorIndex(1, 0)] - 6 * insub[nSubevents - 1][getVectorIndex(4, 1)])); + dynamic_cast(fSubCMList->At(indOffset + 9))->FillProfile(centmult, cmValSub[nSubevents - 1][10], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[nSubevents - 1][4], rn); } // Fill cross terms @@ -976,8 +998,8 @@ void FlowPtContainer::fillCMSubeventProfiles(const double& centmult, const doubl for (int third = 1; third < m; ++third) { if (third > fourth) continue; - if (insub1[getVectorIndex(m, 0)] != 0 && insub2[getVectorIndex(m, 0)] != 0 && cmDen1[m] * cmDen2[m] != 0) - dynamic_cast(fSubCMList->FindObject(Form("cm%i_%i%isub1_%i%isub2", m, first, second, third, fourth)))->FillProfile(centmult, cmVal1[second * (second - 1) / 2 + second - first + 1] * cmVal2[fourth * (fourth - 1) / 2 + fourth - third + 1], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDen1[m] * cmDen2[m], rn); + if (insub[0][getVectorIndex(m, 0)] != 0 && insub[nSubevents - 1][getVectorIndex(m, 0)] != 0 && cmDenSub[0][m] * cmDenSub[nSubevents - 1][m] != 0) + dynamic_cast(fSubCMList->FindObject(Form("cm%i_%i%isub1_%i%isub2", m, first, second, third, fourth)))->FillProfile(centmult, cmValSub[0][second * (second - 1) / 2 + second - first + 1] * cmValSub[nSubevents - 1][fourth * (fourth - 1) / 2 + fourth - third + 1], (fEventWeight == EventWeight::UnityWeight) ? 1.0 : cmDenSub[0][m] * cmDenSub[nSubevents - 1][m], rn); } } } @@ -1852,3 +1874,19 @@ TH1* FlowPtContainer::raiseHistToPower(TH1* inh, double p) } return reth; } +void FlowPtContainer::getSubevents(int k, int n, std::vector& current, std::vector>& subevents) +{ + if (n == 1) { + // Last box gets all remaining objects + current.push_back(k); + subevents.push_back(current); + current.pop_back(); + return; + } + + for (int i = 0; i <= k; ++i) { + current.push_back(i); + getSubevents(k - i, n - 1, current, subevents); + current.pop_back(); + } +} diff --git a/PWGCF/GenericFramework/Core/FlowPtContainer.h b/PWGCF/GenericFramework/Core/FlowPtContainer.h index 86bd676d7d8..67bf59d4cf5 100644 --- a/PWGCF/GenericFramework/Core/FlowPtContainer.h +++ b/PWGCF/GenericFramework/Core/FlowPtContainer.h @@ -51,12 +51,13 @@ class FlowPtContainer : public TNamed void initialise(int nbinsx, double* xbins, const int& m, const GFWCorrConfigs& configs, const int& nsub = 10); void initialise(int nbinsx, double xlow, double xhigh, const int& m, const GFWCorrConfigs& configs, const int& nsub = 10); // initial pt-pt correlations with two subevents - void initialiseSubevent(const o2::framework::AxisSpec axis, const int& m, const int& nsub = 10); - void initialiseSubevent(int nbinsx, double* xbins, const int& m, const int& nsub = 10); - void initialiseSubevent(int nbinsx, double xlow, double xhigh, const int& m, const int& nsub = 10); + void initialiseSubevent(const o2::framework::AxisSpec axis, const int& m, const int& nsubev = 2, const int& nsub = 10); + void initialiseSubevent(int nbinsx, double* xbins, const int& m, const int& nsubev = 2, const int& nsub = 10); + void initialiseSubevent(int nbinsx, double xlow, double xhigh, const int& m, const int& nsubev = 2, const int& nsub = 10); void fill(const double& w, const double& pt); - void fillSub1(const double& w, const double& pt); - void fillSub2(const double& w, const double& pt); + void fillSub(const double& w, const double& pt, int subIndex); + void fillSub1(const double& w, const double& pt) { fillSub(w, pt, 0); } + void fillSub2(const double& w, const double& pt) { fillSub(w, pt, nSubevents - 1); } void fillArray(FillType a, FillType b, double c, double d); int getVectorIndex(const int i, const int j) { return j * (mpar + 1) + i; } // index for 2d array for storing pt correlations int getVectorIndex(const int i, const int j, const int k, const int l) { return i + j * 3 + k * 3 * 3 + l * 3 * 3 * 5; } // index for 4d array for std vnpt correlation - size 3x3x3x3 @@ -134,16 +135,14 @@ class FlowPtContainer : public TNamed { sumP.clear(); sumP.resize((mpar + 1) * (mpar + 1)); - insub1.clear(); - insub1.resize((mpar + 1) * (mpar + 1)); - insub2.clear(); - insub2.resize((mpar + 1) * (mpar + 1)); + insub.clear(); + insub.resize(nSubevents, std::vector((mpar + 1) * (mpar + 1))); cmVal.clear(); - cmVal1.clear(); - cmVal2.clear(); + cmValSub.clear(); + cmValSub.resize(nSubevents); cmDen.clear(); - cmDen1.clear(); - cmDen2.clear(); + cmDenSub.clear(); + cmDenSub.resize(nSubevents); fillCounter = 0; arr.clear(); arr.resize(3 * 3 * 5 * 5, {0.0, 0.0}); @@ -160,6 +159,7 @@ class FlowPtContainer : public TNamed TList* fCentralMomentList; int mpar; + int nSubevents; //! int fillCounter; //! unsigned int fEventWeight; //! bool fUseCentralMoments; @@ -167,20 +167,15 @@ class FlowPtContainer : public TNamed void mergeBSLists(TList* source, TList* target); TH1* raiseHistToPower(TH1* inh, double p); std::vector sumP; //! - std::vector insub1; //! - std::vector insub2; //! + std::vector> insub; //! std::vector corrNum; //! - std::vector corrNum1; //! - std::vector corrNum2; //! + std::vector> corrNumSub; //! std::vector corrDen; //! - std::vector corrDen1; //! - std::vector corrDen2; //! + std::vector> corrDenSub; //! std::vector cmVal; //! - std::vector cmVal1; //! - std::vector cmVal2; //! + std::vector> cmValSub; //! std::vector cmDen; //! - std::vector cmDen1; //! - std::vector cmDen2; //! + std::vector> cmDenSub; //! std::vector> arr; //! std::vector warr; //! std::vector fCovFirstIndex; //! @@ -224,6 +219,8 @@ class FlowPtContainer : public TNamed double getStdABDDD(T& inarr); private: + std::vector> subevents; + void getSubevents(int k, int n, std::vector& current, std::vector>& subevents); static constexpr float FactorialArray[9] = {1., 1., 2., 6., 24., 120., 720., 5040., 40320.}; static constexpr int SignArray[9] = {1, -1, 1, -1, 1, -1, 1, -1, 1}; ClassDef(FlowPtContainer, 2); diff --git a/PWGCF/GenericFramework/Core/GFWConfig.h b/PWGCF/GenericFramework/Core/GFWConfig.h index 779d06b604e..4bfeb10d98e 100644 --- a/PWGCF/GenericFramework/Core/GFWConfig.h +++ b/PWGCF/GenericFramework/Core/GFWConfig.h @@ -12,13 +12,17 @@ #ifndef PWGCF_GENERICFRAMEWORK_CORE_GFWCONFIG_H_ #define PWGCF_GENERICFRAMEWORK_CORE_GFWCONFIG_H_ +#include "GFW.h" + +#include "Framework/Logger.h" + +#include +#include + #include +#include #include #include -#include -#include -#include -#include "GFW.h" namespace o2 { diff --git a/PWGCF/GenericFramework/Tasks/flowGfwLightIons.cxx b/PWGCF/GenericFramework/Tasks/flowGfwLightIons.cxx index 99ed3ee1a41..e644b6f8692 100644 --- a/PWGCF/GenericFramework/Tasks/flowGfwLightIons.cxx +++ b/PWGCF/GenericFramework/Tasks/flowGfwLightIons.cxx @@ -74,6 +74,7 @@ float nchlow = 0; float nchup = 3000; std::vector centbinning(90); int nBootstrap = 10; +std::vector> etagapsPtPt; GFWRegions regions; GFWCorrConfigs configs; std::vector multGlobalCorrCutPars; @@ -84,8 +85,10 @@ std::vector multGlobalT0ACutPars; std::vector firstRunsOfFill; } // namespace o2::analysis::gfw -struct FlowGfwLightIons { +// Let's restrict it to a maximum of four subevents for pt-pt correlations, otherwise output overhead starts to become too much +static constexpr double LongArrayDouble[4][2] = {{-2.0, -2.0}, {-2.0, -2.0}, {-2.0, -2.0}, {-2.0, -2.0}}; +struct FlowGfwLightIons { O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples") O2_DEFINE_CONFIGURABLE(cfgMpar, int, 4, "Highest order of pt-pt correlations") O2_DEFINE_CONFIGURABLE(cfgCentEstimator, int, 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FT0A") @@ -114,6 +117,7 @@ struct FlowGfwLightIons { O2_DEFINE_CONFIGURABLE(cfgEtaPtPt, float, 0.4, "eta cut for pt-pt correlations used in subevent vn-pt"); O2_DEFINE_CONFIGURABLE(cfgEtaPtPtGap, float, 0.4, "eta gap for subevent pt-pt correlations"); O2_DEFINE_CONFIGURABLE(cfgEtaPtPtFull, float, 0.8, "eta cut for pure pt-pt correlations"); + Configurable> cfgPtPtGaps{"cfgPtPtGaps", {LongArrayDouble[0], 4, 2, {"subevent 1", "subevent 2", "subevent 3", "subevent 4"}, {"etamin", "etamax"}}, "{etamin,etamax} for all ptpt-subevents"}; O2_DEFINE_CONFIGURABLE(cfgVtxZ, float, 10, "vertex cut (cm)"); O2_DEFINE_CONFIGURABLE(cfgOccupancySelection, int, 2000, "Max occupancy selection, -999 to disable"); struct : ConfigurableGroup { @@ -329,6 +333,16 @@ struct FlowGfwLightIons { } firstRunOfCurrentFill = o2::analysis::gfw::firstRunsOfFill.begin(); + for (int i = 0; i < 4; ++i) { // o2-linter: disable=magic-number (maximum of 4 subevents) + if (cfgPtPtGaps->getData()[i][0] < -1. || cfgPtPtGaps->getData()[i][1] < -1.) + continue; + o2::analysis::gfw::etagapsPtPt.push_back(std::make_pair(cfgPtPtGaps->getData()[i][0], cfgPtPtGaps->getData()[i][1])); + } + + for (const auto& [etamin, etamax] : o2::analysis::gfw::etagapsPtPt) { + LOGF(info, "pt-pt subevent: {%.1f,%.1f}", etamin, etamax); + } + AxisSpec phiAxis = {o2::analysis::gfw::phibins, o2::analysis::gfw::philow, o2::analysis::gfw::phiup, "#phi"}; AxisSpec etaAxis = {o2::analysis::gfw::etabins, -cfgEta, cfgEta, "#eta"}; AxisSpec vtxAxis = {o2::analysis::gfw::vtxZbins, -cfgVtxZ, cfgVtxZ, "Vtx_{z} (cm)"}; @@ -351,8 +365,8 @@ struct FlowGfwLightIons { AxisSpec t0aAxis = {300, 0, 30000, "N_{ch} (T0A)"}; AxisSpec v0aAxis = {800, 0, 80000, "N_{ch} (V0A)"}; AxisSpec multpvAxis = {600, 0, 600, "N_{ch} (PV)"}; - AxisSpec dcaZAXis = {200, -2, 2, "DCA_{z} (cm)"}; - AxisSpec dcaXYAXis = {200, -0.5, 0.5, "DCA_{xy} (cm)"}; + AxisSpec dcaZAxis = {200, -2, 2, "DCA_{z} (cm)"}; + AxisSpec dcaXYAxis = {200, -0.5, 0.5, "DCA_{xy} (cm)"}; std::vector timebinning(289); std::generate(timebinning.begin(), timebinning.end(), [n = -24 / 288., step = 24 / 288.]() mutable { n += step; @@ -387,7 +401,7 @@ struct FlowGfwLightIons { } if (doprocessMCReco || doprocessData) { registry.add("trackQA/before/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); - registry.add("trackQA/before/pt_dcaXY_dcaZ", "", {HistType::kTH3D, {ptAxis, dcaXYAXis, dcaZAXis}}); + registry.add("trackQA/before/pt_dcaXY_dcaZ", "", {HistType::kTH3D, {ptAxis, dcaXYAxis, dcaZAxis}}); registry.add("trackQA/before/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}}); registry.add("trackQA/before/chi2prTPCcls", "#chi^{2}/cluster for the TPC track segment; #chi^{2}/TPC cluster", {HistType::kTH1D, {{100, 0., 5.}}}); registry.add("trackQA/before/chi2prITScls", "#chi^{2}/cluster for the ITS track; #chi^{2}/ITS cluster", {HistType::kTH1D, {{100, 0., 50.}}}); @@ -442,7 +456,6 @@ struct FlowGfwLightIons { registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kIsGoodITSLayersAll, "kIsGoodITSLayersAll"); registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kMultCuts, "after Mult cuts"); registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kTrackCent, "has track + within cent"); - LOGF(info, "eventsel N bins = %d", registry.get(HIST("eventQA/eventSel"))->GetNbinsX()); if (!cfgRunByRun && cfgFillWeights) { registry.add("phi_eta_vtxz_ref", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); } @@ -478,14 +491,14 @@ struct FlowGfwLightIons { fFCptFull->setUseCentralMoments(cfgUseCentralMoments); fFCptFull->setUseGapMethod(true); fFCptFull->initialise(multAxis, cfgMpar, o2::analysis::gfw::configs, cfgNbootstrap); - fFCptFull->initialiseSubevent(multAxis, cfgMpar, cfgNbootstrap); + fFCptFull->initialiseSubevent(multAxis, cfgMpar, o2::analysis::gfw::etagapsPtPt.size(), cfgNbootstrap); fFCptgen->setUseCentralMoments(cfgUseCentralMoments); fFCptgen->setUseGapMethod(true); fFCptgen->initialise(multAxis, cfgMpar, o2::analysis::gfw::configs, cfgNbootstrap); fFCptgenFull->setUseCentralMoments(cfgUseCentralMoments); fFCptgenFull->setUseGapMethod(true); fFCptgenFull->initialise(multAxis, cfgMpar, o2::analysis::gfw::configs, cfgNbootstrap); - fFCptgenFull->initialiseSubevent(multAxis, cfgMpar, cfgNbootstrap); + fFCptgenFull->initialiseSubevent(multAxis, cfgMpar, o2::analysis::gfw::etagapsPtPt.size(), cfgNbootstrap); fPtDepDCAxy = new TF1("ptDepDCAxy", Form("[0]*%s", cfgDCAxy->c_str()), 0.001, 100); fPtDepDCAxy->SetParameter(0, cfgDCAxyNSigma); @@ -864,6 +877,7 @@ struct FlowGfwLightIons { template void fillOutputContainers(const float& centmult, const double& rndm, const int& run = 0) { + std::cout << "1" << std::endl; if (dt == kGen) { fFCptgen->calculateCorrelations(); fFCptgenFull->calculateCorrelations(); @@ -873,6 +887,7 @@ struct FlowGfwLightIons { fFCptFull->calculateCorrelations(); fFCptFull->calculateSubeventCorrelations(); } + std::cout << "2" << std::endl; if (dt == kGen) { fFCptgen->fillPtProfiles(centmult, rndm); fFCptgen->fillCMProfiles(centmult, rndm); @@ -888,7 +903,9 @@ struct FlowGfwLightIons { fFCptFull->fillCMProfiles(centmult, rndm); fFCptFull->fillCMSubeventProfiles(centmult, rndm); } + std::cout << "3" << std::endl; for (uint l_ind = 0; l_ind < corrconfigs.size(); ++l_ind) { + std::cout << "l_ind = " << l_ind << std::endl; if (!corrconfigs.at(l_ind).pTDif) { uint8_t vnptmask = o2::analysis::gfw::configs.GetpTCorrMasks()[l_ind]; auto dnx = fGFW->Calculate(corrconfigs.at(l_ind), 0, kTRUE).real(); @@ -898,8 +915,17 @@ struct FlowGfwLightIons { } auto val = fGFW->Calculate(corrconfigs.at(l_ind), 0, kFALSE).real() / dnx; if (std::abs(val) < 1) { + std::cout << "Filling flow profiles" << std::endl; (dt == kGen) ? fFCgen->FillProfile(corrconfigs.at(l_ind).Head.c_str(), centmult, val, (cfgUseMultiplicityFlowWeights) ? dnx : 1.0, rndm) : fFC->FillProfile(corrconfigs.at(l_ind).Head.c_str(), centmult, val, (cfgUseMultiplicityFlowWeights) ? dnx : 1.0, rndm); + std::cout << "Filling flow-pt profiles:" << std::endl; + std::cout << "centmult = " << centmult << std::endl; + std::cout << "val = " << val << std::endl; + std::cout << "weight = " << ((cfgUseMultiplicityFlowWeights) ? dnx : 1.0) << std::endl; + std::cout << "rndm = " << rndm << std::endl; + std::cout << "mask = " << static_cast(vnptmask) << std::endl; + std::cout << "fFCpt = " << &(*fFCpt) << std::endl; (dt == kGen) ? fFCptgen->fillVnPtProfiles(centmult, val, (cfgUseMultiplicityFlowWeights) ? dnx : 1.0, rndm, vnptmask) : fFCpt->fillVnPtProfiles(centmult, val, (cfgUseMultiplicityFlowWeights) ? dnx : 1.0, rndm, vnptmask); + std::cout << "flow-pt profiles filled" << std::endl; if (cfgRunByRun && cfgFillFlowRunByRun && dt != kGen && l_ind == 0) { tpfsList[run][pfCorr22]->Fill(centmult, val, (cfgUseMultiplicityFlowWeights) ? dnx : 1.0); } @@ -917,6 +943,7 @@ struct FlowGfwLightIons { (dt == kGen) ? fFCgen->FillProfile(Form("%s_pt_%i", corrconfigs.at(l_ind).Head.c_str(), i), centmult, val, (cfgUseMultiplicityFlowWeights) ? dnx : 1.0, rndm) : fFC->FillProfile(Form("%s_pt_%i", corrconfigs.at(l_ind).Head.c_str(), i), centmult, val, (cfgUseMultiplicityFlowWeights) ? dnx : 1.0, rndm); } } + std::cout << "4" << std::endl; return; } @@ -1182,10 +1209,13 @@ struct FlowGfwLightIons { (dt == kGen) ? fFCptgenFull->fill(1., track.pt()) : fFCptFull->fill(weff, track.pt()); // Fill pt-pt correlations in subevents - if (track.eta() < -cfgEtaPtPtGap && track.eta() > -cfgEtaPtPtFull) - (dt == kGen) ? fFCptgenFull->fillSub1(weff, track.pt()) : fFCptFull->fillSub1(weff, track.pt()); - if (track.eta() > cfgEtaPtPtGap && track.eta() < cfgEtaPtPtFull) - (dt == kGen) ? fFCptgenFull->fillSub2(weff, track.pt()) : fFCptFull->fillSub2(weff, track.pt()); + std::size_t index = 0; + for (const auto& [etamin, etamax] : o2::analysis::gfw::etagapsPtPt) { + if (etamin < track.eta() && track.eta() < etamax) { + (dt == kGen) ? fFCptgenFull->fillSub(weff, track.pt(), index) : fFCptFull->fillSub(weff, track.pt(), index); + } + ++index; + } return; }