diff --git a/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCumulants.cxx b/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCumulants.cxx index 85e7b44cf06..c880ec3d9a0 100644 --- a/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCumulants.cxx +++ b/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCumulants.cxx @@ -199,6 +199,7 @@ static constexpr int NumTwoPCorrBins = NumHarmonics; static constexpr int NumFourPCorrBins = (NumHarmonics * (NumHarmonics + 1)) / 2; // , , , , , , with , in this case 6 static constexpr int NumSixPCorrBins = (NumHarmonics * (NumHarmonics - 1) * (NumHarmonics - 2)) / 6; // , without , in this case 1 static constexpr std::array, NumFourPCorrBins> FourPHarmonicIndex{{{2, 2}, {2, 3}, {2, 4}, {3, 3}, {3, 4}, {4, 4}}}; +static constexpr int NumEtaGap = 3; // *) Main task: struct MultiparticleCumulants { // this name is used in lower-case format to name the TDirectoryFile in AnalysisResults.root @@ -237,7 +238,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam Configurable cfTpcNClsFoundCutSwitch{"cfTpcNClsFoundCutSwitch", true, ""}; Configurable cfDCAXYCutSwitch{"cfDCAXYCutSwitch", true, ""}; Configurable cfDCAZCutSwitch{"cfDCAZCutSwitch", true, ""}; - Configurable cfEtaGapSwitch{"cfEtaGapSwitch", true, "Eta gap switch"}; + Configurable cfHasTpcItsCutSwitch{"cfHasTpcItsCutSwitch", false, ""}; // *) Event cut Configurable> cfVertexZCut{"cfVertexZCut", {-10., 10.}, "vertex z position range: {min, max}[cm]"}; @@ -259,7 +260,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam Configurable> cfTpcNClsFoundCut{"cfTpcNClsFoundCut", {70., 160.}, "range of found TPC clusters for this track geometry: {min, max}"}; Configurable> cfDCAXYCut{"cfDCAXYCut", {-3.2, 3.2}, "range of distance-of-closest-approach (DCA) of the extrapolated track to the primary position in XY-direction: {min, max}[cm]"}; Configurable> cfDCAZCut{"cfDCAZCut", {-2.4, 2.4}, "range of distance-of-closest-approach (DCA) of the extrapolated track to the primary position in Z-direction: {min, max}[cm]"}; - Configurable cfEtaGap{"cfEtaGap", 1., "|dEta| > gap"}; + Configurable> cfHasTpcItsCut{"cfHasTpcItsCut", {1, 1}, "1 to keep tracks that have {TPC, ITS} information"}; + Configurable> cfEtaGap{"cfEtaGap", {0.4, 0.8, 1.0}, "|dEta| > gap"}; // *) Others Configurable cfFileWithWeights{"cfFileWithWeights", "/scratch3/go52dab/O2tutorial/tutorial3-6/weights.root", "path to external ROOT file which holds all particle weights in O2 format"}; @@ -305,7 +307,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam bool fTpcNClsFoundCutSwitch = true; bool fDCAXYCutSwitch = true; bool fDCAZCutSwitch = true; - bool fEtaGapSwitch = true; + bool fHasTpcItsCutSwitch = true; std::vector fVertexZCut = {-10., 10.}; std::vector fCentCut = {10., 20.}; @@ -319,7 +321,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam std::vector fTpcNClsFoundCut = {70., 160.}; std::vector fDCAXYCut = {-3.2, 3.2}; std::vector fDCAZCut = {-2.4, 2.4}; - float fEtaGap = 1.; + std::vector fHasTpcItsCut = {1, 1}; + std::vector fEtaGap = {0.4, 0.8, 1.0}; std::vector fPtBins = {1000, 0., 100.}; std::vector fPhiBins = {1000, 0., o2::constants::math::TwoPI}; @@ -396,10 +399,10 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam std::array, MaxHarmonic> fQvectorBefore; std::array, MaxHarmonic> fQvectorAfter; - std::array, MaxHarmonic> fQvectorBeforeA; // Q-vector with eta gap - std::array, MaxHarmonic> fQvectorBeforeB; // Q-vector with eta gap - std::array, MaxHarmonic> fQvectorAfterA; // Q-vector with eta gap - std::array, MaxHarmonic> fQvectorAfterB; // Q-vector with eta gap + std::array, MaxHarmonic>, NumEtaGap> fQvectorBeforeA; + std::array, MaxHarmonic>, NumEtaGap> fQvectorBeforeB; + std::array, MaxHarmonic>, NumEtaGap> fQvectorAfterA; + std::array, MaxHarmonic>, NumEtaGap> fQvectorAfterB; } mcc; struct MultiparticleCorrelationProfile { @@ -409,11 +412,11 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam std::array fTwoParticleCorrelationProfiles{}; // [cut] std::array fFourParticleCorrelationProfiles{}; std::array fSixParticleCorrelationProfiles{}; - std::array fTwoParticleCorrelationGapProfiles{}; + std::array, NumEtaGap> fTwoParticleCorrelationGapProfiles{}; std::array, eBeforeAfter_N>, eBeforeAfter_N> fTwoParticleCorrelationHistograms{}; // [cut][event weight][n] std::array, eBeforeAfter_N>, eBeforeAfter_N> fFourParticleCorrelationHistograms{}; std::array, eBeforeAfter_N>, eBeforeAfter_N> fSixParticleCorrelationHistograms{}; - std::array fTwoParticleCorrelationGapHistograms{}; + std::array, NumEtaGap> fTwoParticleCorrelationGapHistograms{}; } mc; struct EventByEventQuantities { @@ -435,7 +438,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam std::array, eBeforeAfter_N> fSixParticleCorrelationMaxEbye{}; } ebye; - template // rm should be rec and sim + template bool ctEventCuts(T1 const& collision, T2 const& rlCollisionCentAll, T3 const& rlCollisionMultAll, T4 const& rlCollisionNumContrib) { bool pass = true; @@ -552,6 +555,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam bool bTpcNClsFoundCut = true; bool bDCAXYCut = true; bool bDCAZCut = true; + bool bHasTpcCut = true; + bool bHasItsCut = true; // *) For rec event and sim event bPtCut = track.pt() < tc.fPtCut[1] && track.pt() > tc.fPtCut[0]; @@ -568,6 +573,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam track.dcaXY() > tc.fDCAXYCut[0]; bDCAZCut = track.dcaZ() < tc.fDCAZCut[1] && track.dcaZ() > tc.fDCAZCut[0]; + bHasTpcCut = track.hasTPC() && tc.fHasTpcItsCut[0]; + bHasItsCut = track.hasITS() && tc.fHasTpcItsCut[1]; } // *) For sim event only @@ -603,6 +610,10 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam if (tc.fDCAZCutSwitch) { pass &= bDCAZCut; } + if (tc.fHasTpcItsCutSwitch) { + pass &= bHasTpcCut; + pass &= bHasItsCut; + } return pass; } @@ -618,15 +629,15 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam return TComplex::Conjugate(mcc.fQvectorBefore[-n][p]); } if (n >= 0) { - return mcc.fQvectorAfterB[n][p]; + return mcc.fQvectorAfter[n][p]; } - return TComplex::Conjugate(mcc.fQvectorAfterB[-n][p]); + return TComplex::Conjugate(mcc.fQvectorAfter[-n][p]); } - TComplex mccTwo(int n1, int n2, EnBeforeAfter eba) - { - return mccQ(n1, 1, eba) * mccQ(n2, 1, eba) - mccQ(n1 + n2, 2, eba); - } + // TComplex mccTwo(int n1, int n2, EnBeforeAfter eba) + // { + // return mccQ(n1, 1, eba) * mccQ(n2, 1, eba) - mccQ(n1 + n2, 2, eba); + // } template TComplex mccRecursion(int n, std::array harmonic, EnBeforeAfter eba, int mult = 1, int skip = 0) @@ -865,10 +876,13 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam for (int p = 0; p < mcc.MaxPower; p++) { mcc.fQvectorBefore[h][p] = TComplex(0., 0.); mcc.fQvectorAfter[h][p] = TComplex(0., 0.); - mcc.fQvectorBeforeA[h][p] = TComplex(0., 0.); - mcc.fQvectorBeforeB[h][p] = TComplex(0., 0.); - mcc.fQvectorAfterA[h][p] = TComplex(0., 0.); - mcc.fQvectorAfterB[h][p] = TComplex(0., 0.); + + for (int n = 0; n < NumEtaGap; n++) { + mcc.fQvectorBeforeA[n][h][p] = TComplex(0., 0.); + mcc.fQvectorBeforeB[n][h][p] = TComplex(0., 0.); + mcc.fQvectorAfterA[n][h][p] = TComplex(0., 0.); + mcc.fQvectorAfterB[n][h][p] = TComplex(0., 0.); + } } } @@ -959,7 +973,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam mc.fTwoParticleCorrelationProfiles[i] = nullptr; mc.fFourParticleCorrelationProfiles[i] = nullptr; mc.fSixParticleCorrelationProfiles[i] = nullptr; - mc.fTwoParticleCorrelationGapProfiles[i] = nullptr; + mc.fTwoParticleCorrelationGapProfiles[i].fill(nullptr); for (int j = 0; j < eBeforeAfter_N; j++) { // before/after weight mc.fTwoParticleCorrelationHistograms[i][j].fill(nullptr); @@ -972,18 +986,22 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam for (int i = 0; i < eBeforeAfter_N; i++) { // before/after cut mc.fTwoParticleCorrelationProfiles[i] = new TProfile(Form("prof2%sCut", BeforeAfterNames[i + 2]), Form("2-p correlation %s cut", BeforeAfterNames[i]), NumTwoPCorrBins, 0., 1.); mc.fTwoParticleCorrelationProfiles[i]->Sumw2(); - mc.fTwoParticleCorrelationGapProfiles[i] = new TProfile(Form("prof2%sCutGapped", BeforeAfterNames[i + 2]), Form("2-p correlation %s cut, gapped", BeforeAfterNames[i]), NumTwoPCorrBins, 0., 1.); - mc.fTwoParticleCorrelationProfiles[i]->Sumw2(); - mc.fTwoParticleCorrelationGapProfiles[i]->Sumw2(); + for (int n = 0; n < NumEtaGap; n++) { + mc.fTwoParticleCorrelationGapProfiles[n][i] = new TProfile(Form("prof2%sCutGapped%.1f", BeforeAfterNames[i + 2], tc.fEtaGap[n]), Form("2-p correlation %s cut, %.1f gapped", BeforeAfterNames[i], tc.fEtaGap[n]), NumTwoPCorrBins, 0., 1.); + mc.fTwoParticleCorrelationGapProfiles[n][i]->Sumw2(); + } for (int k = 0; k < NumTwoPCorrBins; k++) { // v2, v3, v4 mc.fTwoParticleCorrelationProfiles[i]->GetXaxis()->SetBinLabel(k + 1, Form("", k + 2)); - mc.fTwoParticleCorrelationGapProfiles[i]->GetXaxis()->SetBinLabel(k + 1, Form("", k + 2)); - if (static_cast(i)) { - mc.fTwoParticleCorrelationGapHistograms[k] = new TH1D(Form("hist2v%dAfterCutWithGap", k + 2), Form("2-p correlation v%d^2 after cut with gap", k + 2), static_cast(tc.fTwoParticleCorrBins[0]), tc.fTwoParticleCorrBins[1], tc.fTwoParticleCorrBins[2]); - mc.fTwoParticleCorrelationGapHistograms[k]->Sumw2(); - mc.fTwoParticleCorrelationGapHistograms[k]->SetOption("HIST"); + for (int n = 0; n < NumEtaGap; n++) { + mc.fTwoParticleCorrelationGapProfiles[n][i]->GetXaxis()->SetBinLabel(k + 1, Form("", k + 2)); + if (static_cast(i)) { + mc.fTwoParticleCorrelationGapHistograms[n][k] = new TH1D(Form("hist2v%dAfterCutWithGap%.1f", k + 2, tc.fEtaGap[n]), Form("2-p correlation v%d^2 after cut with gap %.1f", k + 2, tc.fEtaGap[n]), static_cast(tc.fTwoParticleCorrBins[0]), tc.fTwoParticleCorrBins[1], tc.fTwoParticleCorrBins[2]); + mc.fTwoParticleCorrelationGapHistograms[n][k]->Sumw2(); + mc.fTwoParticleCorrelationGapHistograms[n][k]->SetOption("HIST"); + } } + for (int j = 0; j < eBeforeAfter_N; j++) { // before/after weight mc.fTwoParticleCorrelationHistograms[i][j][k] = new TH1D(Form("hist2v%d%sCut%sWeight", k + 2, BeforeAfterNames[i + 2], BeforeAfterNames[j + 2]), Form("2-p correlation v%d^2 %s cut %s weight", k + 2, BeforeAfterNames[i], BeforeAfterNames[j]), static_cast(tc.fTwoParticleCorrBins[0]), tc.fTwoParticleCorrBins[1], tc.fTwoParticleCorrBins[2]); if (static_cast(j)) { @@ -993,11 +1011,15 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam mc.fMultiparticleCorrelationBackupList->Add(mc.fTwoParticleCorrelationHistograms[i][j][k]); } if (static_cast(i)) { - mc.fMultiparticleCorrelationBackupList->Add(mc.fTwoParticleCorrelationGapHistograms[k]); + for (int n = 0; n < NumEtaGap; n++) { + mc.fMultiparticleCorrelationBackupList->Add(mc.fTwoParticleCorrelationGapHistograms[n][k]); + } } } mc.fMultiparticleCorrelationByRunMap[ebye.fRunNumber]->Add(mc.fTwoParticleCorrelationProfiles[i]); - mc.fMultiparticleCorrelationByRunMap[ebye.fRunNumber]->Add(mc.fTwoParticleCorrelationGapProfiles[i]); + for (int n = 0; n < NumEtaGap; n++) { + mc.fMultiparticleCorrelationByRunMap[ebye.fRunNumber]->Add(mc.fTwoParticleCorrelationGapProfiles[n][i]); + } } // Define 4p profiles and histograms: @@ -1117,31 +1139,27 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam ev.fEventHistograms[eNumContrib][eRec][eBefore]->Fill(rlCollisionNumContrib); // Fill centrality correlation histograms before cut: - if (tc.fCentCorrCutSwitch) { - for (int i = 0; i < eCentEstm_N; i++) { - for (int j = i + 1; j < eCentEstm_N; j++) { - auto* h = cr.fCorrHistograms[eCorrCent][i][j][eBefore]; - if (!h) { - LOGF(fatal, "Missing histogram cr.fCorrHistograms[eCorrCent][%d][%d][eBefore]", i, j); - } - if (rlCollisionCentAll[i] >= 0. && rlCollisionCentAll[j] >= 0.) { - h->Fill(rlCollisionCentAll[i], rlCollisionCentAll[j]); - } + for (int i = 0; i < eCentEstm_N; i++) { + for (int j = i + 1; j < eCentEstm_N; j++) { + auto* h = cr.fCorrHistograms[eCorrCent][i][j][eBefore]; + if (!h) { + LOGF(fatal, "Missing histogram cr.fCorrHistograms[eCorrCent][%d][%d][eBefore]", i, j); + } + if (rlCollisionCentAll[i] >= 0. && rlCollisionCentAll[j] >= 0.) { + h->Fill(rlCollisionCentAll[i], rlCollisionCentAll[j]); } } } // Fill multiplicity correlation histograms before cut: - if (tc.fMultCorrCutSwitch) { - for (int i = 0; i < eMultEstm_N; i++) { - for (int j = i + 1; j < eMultEstm_N; j++) { - auto* h = cr.fCorrHistograms[eCorrMult][i][j][eBefore]; - if (!h) { - LOGF(fatal, "Missing histogram cr.fCorrHistograms[eCorrMult][%d][%d][eBefore]", i, j); - } - if (rlCollisionMultAll[i] >= 0. && rlCollisionMultAll[j] >= 0.) { - h->Fill(rlCollisionMultAll[i], rlCollisionMultAll[j]); - } + for (int i = 0; i < eMultEstm_N; i++) { + for (int j = i + 1; j < eMultEstm_N; j++) { + auto* h = cr.fCorrHistograms[eCorrMult][i][j][eBefore]; + if (!h) { + LOGF(fatal, "Missing histogram cr.fCorrHistograms[eCorrMult][%d][%d][eBefore]", i, j); + } + if (rlCollisionMultAll[i] >= 0. && rlCollisionMultAll[j] >= 0.) { + h->Fill(rlCollisionMultAll[i], rlCollisionMultAll[j]); } } } @@ -1214,28 +1232,24 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam ev.fEventHistograms[eNumContrib][eRec][eAfter]->Fill(rlCollisionNumContrib); // Fill centrality correlation histograms after cut: - if (tc.fCentCorrCutSwitch) { - for (int i = 0; i < eCentEstm_N; i++) { - for (int j = i + 1; j < eCentEstm_N; j++) { - auto* h = cr.fCorrHistograms[eCorrCent][i][j][eAfter]; - if (!h) { - LOGF(fatal, "Missing histogram cr.fCorrHistograms[eCorrCent][%d][%d][eAfter]", i, j); - } - h->Fill(rlCollisionCentAll[i], rlCollisionCentAll[j]); + for (int i = 0; i < eCentEstm_N; i++) { + for (int j = i + 1; j < eCentEstm_N; j++) { + auto* h = cr.fCorrHistograms[eCorrCent][i][j][eAfter]; + if (!h) { + LOGF(fatal, "Missing histogram cr.fCorrHistograms[eCorrCent][%d][%d][eAfter]", i, j); } + h->Fill(rlCollisionCentAll[i], rlCollisionCentAll[j]); } } // Fill multiplicity correlation histograms after cut: - if (tc.fMultCorrCutSwitch) { - for (int i = 0; i < eMultEstm_N; i++) { - for (int j = i + 1; j < eMultEstm_N; j++) { - auto* h = cr.fCorrHistograms[eCorrMult][i][j][eAfter]; - if (!h) { - LOGF(fatal, "Missing histogram cr.fCorrHistograms[eCorrMult][%d][%d][eAfter]", i, j); - } - h->Fill(rlCollisionMultAll[i], rlCollisionMultAll[j]); + for (int i = 0; i < eMultEstm_N; i++) { + for (int j = i + 1; j < eMultEstm_N; j++) { + auto* h = cr.fCorrHistograms[eCorrMult][i][j][eAfter]; + if (!h) { + LOGF(fatal, "Missing histogram cr.fCorrHistograms[eCorrMult][%d][%d][eAfter]", i, j); } + h->Fill(rlCollisionMultAll[i], rlCollisionMultAll[j]); } } @@ -1252,10 +1266,10 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam int nTracksBefore = tracks.size(); int nTracksAfter = 0; - int nTracksBeforeA = 0; - int nTracksBeforeB = 0; - int nTracksAfterA = 0; - int nTracksAfterB = 0; + std::vector nTracksBeforeA = {0, 0, 0}; + std::vector nTracksBeforeB = {0, 0, 0}; + std::vector nTracksAfterA = {0, 0, 0}; + std::vector nTracksAfterB = {0, 0, 0}; // Calculate Q-vectors for available angles and weights: double dEta = 0.; @@ -1272,10 +1286,6 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam if constexpr (rs == eRec || rs == eRecAndSim) { - // Fill phi/pt real histogram with this run number: - wt.fPhiByRunMap.at(ebye.fRunNumber)->Fill(track.phi()); - wt.fPtRealByRunMap.at(ebye.fRunNumber)->Fill(track.pt()); - // Fill track histograms before cut: pc.fParticleHistograms[ePt][eRec][eBefore]->Fill(track.pt()); pc.fParticleHistograms[ePhi][eRec][eBefore]->Fill(track.phi()); @@ -1306,26 +1316,30 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam wPhiToPowerP = std::pow(wPhi, p); } mcc.fQvectorBefore[h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); - if (tc.fEtaGapSwitch) { - if (dEta < -tc.fEtaGap / 2.) { - mcc.fQvectorBeforeA[h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); - } else if (dEta > tc.fEtaGap / 2.) { - mcc.fQvectorBeforeB[h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); + for (int n = 0; n < NumEtaGap; n++) { + if (dEta < -tc.fEtaGap[n] / 2.) { + mcc.fQvectorBeforeA[n][h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); + } else if (dEta > tc.fEtaGap[n] / 2.) { + mcc.fQvectorBeforeB[n][h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); } } } } - if (tc.fEtaGapSwitch) { - if (dEta < -tc.fEtaGap / 2.) { - nTracksBeforeA += 1; - } else if (dEta > tc.fEtaGap / 2.) { - nTracksBeforeB += 1; + for (int n = 0; n < NumEtaGap; n++) { + if (dEta < -tc.fEtaGap[n] / 2.) { + nTracksBeforeA[n] += 1; + } else if (dEta > tc.fEtaGap[n] / 2.) { + nTracksBeforeB[n] += 1; } } if (ctParticleCuts(track)) { + // Fill phi/pt real histogram with this run number: + wt.fPhiByRunMap.at(ebye.fRunNumber)->Fill(track.phi()); + wt.fPtRealByRunMap.at(ebye.fRunNumber)->Fill(track.pt()); + // Fill particle histograms after cut: pc.fParticleHistograms[ePt][eRec][eAfter]->Fill(track.pt()); pc.fParticleHistograms[ePhi][eRec][eAfter]->Fill(track.phi()); @@ -1337,22 +1351,24 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam wPhiToPowerP = std::pow(wPhi, p); } mcc.fQvectorAfter[h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); - if (tc.fEtaGapSwitch) { - if (dEta < -tc.fEtaGap / 2.) { - mcc.fQvectorAfterA[h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); - } else if (dEta > tc.fEtaGap / 2.) { - mcc.fQvectorAfterB[h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); + + for (int n = 0; n < NumEtaGap; n++) { + if (dEta < -tc.fEtaGap[n] / 2.) { + mcc.fQvectorAfterA[n][h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); + } else if (dEta > tc.fEtaGap[n] / 2.) { + mcc.fQvectorAfterB[n][h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); } } } } nTracksAfter += 1; - if (tc.fEtaGapSwitch) { - if (dEta < -tc.fEtaGap / 2.) { - nTracksAfterA += 1; - } else if (dEta > tc.fEtaGap / 2.) { - nTracksAfterB += 1; + + for (int n = 0; n < NumEtaGap; n++) { + if (dEta < -tc.fEtaGap[n] / 2.) { + nTracksAfterA[n] += 1; + } else if (dEta > tc.fEtaGap[n] / 2.) { + nTracksAfterB[n] += 1; } } } @@ -1372,14 +1388,19 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // Corresponding MC truth simulated particle auto mcparticle = track.mcParticle(); - // Fill pt MC sim histogram with this run number: - wt.fPtMCByRunMap.at(ebye.fRunNumber)->Fill(mcparticle.pt()); - // Fill MC particle histograms before cut: pc.fParticleHistograms[ePt][eSim][eBefore]->Fill(mcparticle.pt()); pc.fParticleHistograms[ePhi][eSim][eBefore]->Fill(mcparticle.phi()); if (ctParticleCuts(mcparticle)) { + // Fill pt MC sim histogram with this run number: + double diffPt = std::abs(track.pt() - mcparticle.pt()); + double diffPtMax = (tc.fPtBins[2] - tc.fPtBins[1]) / tc.fPtBins[0]; + if (diffPt < diffPtMax) { + wt.fPtMCByRunMap.at(ebye.fRunNumber)->Fill(mcparticle.pt()); + } else { + LOGF(info, "|RecPt - SimPt| = %e > %e", diffPt, diffPtMax); + } // Fill MC particle histograms after cut: pc.fParticleHistograms[ePt][eSim][eAfter]->Fill(mcparticle.pt()); pc.fParticleHistograms[ePhi][eSim][eAfter]->Fill(mcparticle.phi()); @@ -1445,28 +1466,34 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam LOGF(warning, "cent=%f, nTracksAfter = %d, wTwoRecursionAfter = %e", rlCollisionCent, nTracksAfter, wTwoRecursionAfter); } - // Before cut, with gap: - TComplex qvba = mcc.fQvectorBeforeA[mcc.h2][0]; - TComplex qvbb = mcc.fQvectorBeforeB[mcc.h2][0]; - double gapb = 0.; - if (nTracksBeforeA * nTracksBeforeB > 0) { - gapb = (qvba * TComplex::Conjugate(qvbb)).Re() / (nTracksBeforeA * nTracksBeforeB); - mc.fTwoParticleCorrelationGapProfiles[eBefore]->Fill(mc.fTwoParticleCorrelationGapProfiles[eBefore]->GetXaxis()->GetBinCenter(i + 1), gapb, nTracksBeforeA * nTracksBeforeB); - // mc.fTwoParticleCorrelationGapHistograms[i]->Fill(gapb, nTracksBeforeA * nTracksBeforeB); - } else { - LOGF(warning, "cent=%f, nTracksBeforeA = %d, nTracksBeforeB = %d", rlCollisionCent, nTracksBeforeA, nTracksBeforeB); - } + for (int n = 0; n < NumEtaGap; n++) { + + // Before cut, with gap: + TComplex qvba = mcc.fQvectorBeforeA[n][mcc.h2][1]; + TComplex qvbb = mcc.fQvectorBeforeB[n][mcc.h2][1]; + double wba = mcc.fQvectorBeforeA[n][0][1].Re(); + double wbb = mcc.fQvectorBeforeB[n][0][1].Re(); + double gapb = 0.; + if (wba * wbb > 0) { // nTracksBeforeA[n] * nTracksBeforeB[n] > 0 + gapb = (qvba * TComplex::Conjugate(qvbb)).Re() / (wba * wbb); + mc.fTwoParticleCorrelationGapProfiles[n][eBefore]->Fill(mc.fTwoParticleCorrelationGapProfiles[n][eBefore]->GetXaxis()->GetBinCenter(i + 1), gapb, wba * wbb); + } else { + LOGF(warning, "etagap=%f, cent=%f, nTracksBeforeA = %d, nTracksBeforeB = %d", tc.fEtaGap[n], rlCollisionCent, nTracksBeforeA[n], nTracksBeforeB[n]); + } - // After cut, with gap: - TComplex qvaa = mcc.fQvectorAfterA[mcc.h2][0]; - TComplex qvab = mcc.fQvectorAfterB[mcc.h2][0]; - double gapa = 0.; - if (nTracksAfterA * nTracksAfterB > 0) { - gapa = (qvaa * TComplex::Conjugate(qvab)).Re() / (nTracksAfterA * nTracksAfterB); - mc.fTwoParticleCorrelationGapProfiles[eAfter]->Fill(mc.fTwoParticleCorrelationGapProfiles[eAfter]->GetXaxis()->GetBinCenter(i + 1), gapa, nTracksAfterA * nTracksAfterB); - mc.fTwoParticleCorrelationGapHistograms[i]->Fill(gapa, nTracksAfterA * nTracksAfterB); - } else { - LOGF(warning, "cent=%f, nTracksAfterA = %d, nTracksAfterB = %d", rlCollisionCent, nTracksAfterA, nTracksAfterB); + // After cut, with gap: + TComplex qvaa = mcc.fQvectorAfterA[n][mcc.h2][1]; + TComplex qvab = mcc.fQvectorAfterB[n][mcc.h2][1]; + double waa = mcc.fQvectorAfterA[n][0][1].Re(); + double wab = mcc.fQvectorAfterB[n][0][1].Re(); + double gapa = 0.; + if (waa * wab > 0) { + gapa = (qvaa * TComplex::Conjugate(qvab)).Re() / (waa * wab); + mc.fTwoParticleCorrelationGapProfiles[n][eAfter]->Fill(mc.fTwoParticleCorrelationGapProfiles[n][eAfter]->GetXaxis()->GetBinCenter(i + 1), gapa, waa * wab); + mc.fTwoParticleCorrelationGapHistograms[n][i]->Fill(gapa, waa * wab); + } else { + LOGF(warning, "etagap=%f, cent=%f, nTracksAfterA = %d, nTracksAfterB = %d", tc.fEtaGap[n], rlCollisionCent, nTracksAfterA[n], nTracksAfterB[n]); + } } } @@ -1922,7 +1949,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam tc.fTpcNClsFoundCutSwitch = cfTpcNClsFoundCutSwitch; tc.fDCAXYCutSwitch = cfDCAXYCutSwitch; tc.fDCAZCutSwitch = cfDCAZCutSwitch; - tc.fEtaGapSwitch = cfEtaGapSwitch; + tc.fHasTpcItsCutSwitch = cfHasTpcItsCutSwitch; tc.fVertexZCut = cfVertexZCut; tc.fCentCut = cfCentCut; @@ -1941,6 +1968,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam tc.fTpcNClsFoundCut = cfTpcNClsFoundCut; tc.fDCAXYCut = cfDCAXYCut; tc.fDCAZCut = cfDCAZCut; + tc.fHasTpcItsCut = cfHasTpcItsCut; tc.fEtaGap = cfEtaGap; tc.fPtBins = cfPtBins;