From cf19ad525fa3e1bffdb1316e8b11e26785a22382 Mon Sep 17 00:00:00 2001 From: Yuanjun Mei Date: Wed, 11 Nov 2026 10:28:44 +0100 Subject: [PATCH] Add new cuts and histograms with update with cuts ordering --- .../Tasks/multiparticleCorrelationsMei.cxx | 402 +++++++++++------- 1 file changed, 245 insertions(+), 157 deletions(-) diff --git a/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCorrelationsMei.cxx b/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCorrelationsMei.cxx index 4d4ce5047f1..1891910e581 100644 --- a/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCorrelationsMei.cxx +++ b/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCorrelationsMei.cxx @@ -128,7 +128,7 @@ enum EEventCuts { eOccupancyCut, eCorrCentCut, eMultNContribCut, - EEventCuts_N + eEventCuts_N }; enum EParticleCuts { @@ -138,7 +138,8 @@ enum EParticleCuts { eTPCNumClusterCut, eDCAxyCut, eDCAzCut, - EParticleCuts_N + eITSNClsCut, + eParticleCuts_N }; enum EBeforeAfterCuts { @@ -171,6 +172,10 @@ enum EParticleHistograms { eHistPhi, eHistEta, eHistCharge, + eHistTPCNumCluster, + eHistDCAxy, + eHistDCAz, + eHistITSNCls, eParticleHistograms_N }; @@ -205,8 +210,8 @@ static constexpr std::array CentralityEstim static constexpr std::array MultiplicityTablesNames = {"multTPC", "multFV0M", "multFT0C", "multFT0M", "multNTracksPV"}; static constexpr std::array OccupancyEstimatorNames = {"track", "FT0C"}; -static constexpr std::array ParticleHistTitles = {"p_{T}", "#varphi", "#eta", "Particle charge"}; -static constexpr std::array ParticleHistNames = {"eHistPt", "eHistPhi", "eHistEta", "eHistCharge"}; +static constexpr std::array ParticleHistTitles = {"p_{T}", "#varphi", "#eta", "Particle charge", "Number of found TPC clusters", "DCA XY (cm)", "DCA Z (cm)", "track.itsNCls()"}; +static constexpr std::array ParticleHistNames = {"eHistPt", "eHistPhi", "eHistEta", "eHistCharge", "eHistTPCNumCluster", "eHistDCAxy", "eHistDCAz", "eHistITSNCls"}; static constexpr std::array EventHistTitles = {"Centrality", "Multiplicity (trak.size())", "Reference Multiplicity", "Vertex X", "Vertex Y", "Vertex Z", "Impact Parameter", "Number of Contributor", "Occupancy"}; static constexpr std::array EventHistNames = {"eHistCentrality", "eHistMult", "eHistRefMult", "eHistVertexX", "eHistVertexY", "eHistVertexZ", "eHistImpactParameter", "eHistNumContrib", "eHistOccupancy"}; @@ -257,27 +262,32 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to Configurable> cfPtBins{"cfPtBins", {4000, 0., 10.}, "nPtBins, ptMin, ptMax"}; Configurable> cfPhiBins{"cfPhiBins", {180, 0., math::TwoPI}, "nPhiBins, phiMin, phiMax"}; Configurable> cfEtaBins{"cfEtaBins", {800, -3., 3.}, "nEtaBins, etaMin, etaMax"}; + Configurable> cfTPCNumClusterBins{"cfTPCNumClusterBins", {160, 0., 160.}, "nNumClusterBins, numClusterMin, numClusterMax"}; + Configurable> cfITSNClsBins{"cfITSNClsBins", {20, 0., 20.}, "nITSNClsBins, ITSNClsMin, ITSNClsMax"}; + Configurable> cfDCAxyBins{"cfDCAxyBins", {500, -7., 7.}, "nDCAxyBins, DCAxyMin, DCAxyMax"}; + Configurable> cfDCAzBins{"cfDCAzBins", {500, -7., 7.}, "nDCAzBins, DCAzMin, DCAzMax"}; - Configurable> cfQASelTrackBins{"cfQASelTrackBins", {200, 0., 15000., 400, 0., 6000.}, "bins for hist of selected tracks of QA: {nQASelTrackBinsBefore, minBefore, MaxBefore, nQASelTrackBinsAfter, minAfter, MaxAfter}"}; + Configurable> cfQASelTrackBins{"cfQASelTrackBins", {500, 0., 15000., 150, 0., 300.}, "bins for hist of selected tracks of QA: {nQASelTrackBinsBefore, minBefore, MaxBefore, nQASelTrackBinsAfter, minAfter, MaxAfter}"}; // *) Cuts Configurable cfMasterCutSwitch{"cfMasterCutSwitch", true, "switch on or off all cuts"}; // event level cuts - Configurable> cfEventCutSwitch{"cfEventCutSwitch", {"1Sel8", "1IsGoodITSLayersAll", "1NoCollInTimeRangeStandard", "1NoCollInRofStandard", "1IsGoodZvtxFT0vsPV", "1NoHighMultCollInPrevRof", "1NoSameBunchPileUp", "1IsVertexITSTPC", "1Cent", "1VertexZ", "1NContrib", "1Occupancy", "1CorrCent", "1MultNContribCut"}, "event cuts switch, on and off by the first number before name"}; + Configurable> cfEventCutSwitch{"cfEventCutSwitch", {"1Sel8", "1IsGoodITSLayersAll", "1NoCollInTimeRangeStandard", "1NoCollInRofStandard", "1IsGoodZvtxFT0vsPV", "1NoHighMultCollInPrevRof", "1Cent", "1VertexZ", "1Occupancy", "1CorrCent", "1MultNContribCut", "1NoSameBunchPileUp", "1IsVertexITSTPC", "1NContrib"}, "event cuts switch, on and off by the first number before name"}; Configurable> cfVertexZCutRange{"cfVertexZCutRange", {-10., 10.}, "vertex z position range: {min, max}[cm], with convention: min <= Vz <= max"}; Configurable> cfCentCutRange{"cfCentCutRange", {0., 80.}, "centrality range: {min, max}, with convention: min <= cent <= max"}; - Configurable cfNContribCutInf{"cfNContribCutInf", 2, "Cuts on number of tracks used for the vertex, only numContrib > cfNContribCutInf survives"}; - Configurable cfOccupancyCutSup{"cfOccupancyCutSup", 6000, "Cuts on number of tracks used for the vertex, only occupancy < cfOccupancyCutSup survives"}; + Configurable cfNContribCutInf{"cfNContribCutInf", 2, "cuts on number of tracks used for the vertex, only numContrib > cfNContribCutInf survives"}; + Configurable cfOccupancyCutSup{"cfOccupancyCutSup", 6000, "cuts on occupancy, only occupancy < cfOccupancyCutSup survives"}; Configurable> cfCorrCent{"cfCorrCent", {0, 2}, "choose the pair of centrality estimators for cut: 0=FT0C, 1=FT0M, 2=FV0A, 3=NTPV"}; - Configurable cfCorrCentCutRatio{"cfCorrCentCutRatio", 0.5, "correlation cut on centrality, only 2*|cent1 - cent2|/(cent1 + cent2) < cfCorrCentCutRatio survives"}; - Configurable> cfMultNContribBound{"cfMultNContribBound", {"0.160 * x + 153.563", "0.095 * x + -157.476"}, "geometric cuts on correlation of mult vs. nContrib switch, convention {upper bound, lower bound}"}; + Configurable cfCorrCentCutPercentage{"cfCorrCentCutPercentage", 30, "correlation cut on centrality, only |cent1 - cent2| < cfCorrCentCutPercentage survives"}; + Configurable> cfMultNContribBound{"cfMultNContribBound", {"0.19*x + 320", "0.056*x - 324"}, "geometric cuts on correlation of mult vs. nContrib switch, convention {upper bound, lower bound}"}; // particle level cuts - Configurable> cfParticleCutSwitch{"cfParticleCutSwitch", {"1Pt", "1Eta", "1Charge", "1TPCNumCluster", "1DCAxy", "1DCAz"}, "particle cuts switch, on and off by the first number before name"}; + Configurable> cfParticleCutSwitch{"cfParticleCutSwitch", {"1Pt", "1Eta", "1TPCNumCluster", "1DCAxy", "1DCAz", "1ITSNCls", "1Charge"}, "particle cuts switch, on and off by the first number before name"}; Configurable> cfPtCutRange{"cfPtCutRange", {0.2, 5.}, "pt cut range: {min, max}[GeV/c], with convention: min <= pt <= max"}; Configurable> cfEtaCutRange{"cfEtaCutRange", {-0.8, 0.8}, "eta cut range: {min, max}, with convention: min <= eta <= max"}; - Configurable> cfTPCNumClusterCutRange{"cfTPCNumClusterCutRange", {70, 160}, "range of cut on the number of clusters associated with a reconstructed track: {min, max}, with convention: min <= TPCNcls <= max"}; + Configurable cfTPCNumClusterCutInf{"cfTPCNumClusterCutInf", 70, "inf of cut on the number of found TPC clusters associated with a reconstructed track, with convention: cfTPCNumClusterCutInf <= TPCNumCluster"}; + Configurable cfITSNClsInf{"cfITSNClsInf", 5, "inf of cut on the number of found ITS clusters associated with a reconstructed track: {min, max}, with convention: cfITSNClsInf <= itsNCls"}; Configurable> cfDCAxyCutRange{"cfDCAxyCutRange", {-3.2, 3.2}, "range of cut on signed transverse distance of closest approach between the reconstructed track and the primary interaction vertex in x-y plane: {min, max}, with convention: min <= dcaXY <= max"}; Configurable> cfDCAzCutRange{"cfDCAzCutRange", {-2.4, 2.4}, "range of cut on signed transverse distance of closest approach between the reconstructed track and the primary interaction vertex in z direction: {min, max}, with convention: min <= dcaZ <= max"}; @@ -290,23 +300,27 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to Configurable> cfIVHarmonics{"cfIVHarmonics", {0.06, 0.07, 0.08}, "harmonics for internal validation: {v_2, v_3, v_4, ...}"}; Configurable cfRNGSeed{"cfRNGSeed", 0, "seed for random number generator"}; Configurable> cfQASelTrackCentRange{"cfQASelTrackCentRange", {60., 70.}, "centrality range for hist of selected tracks for QA"}; + Configurable cfChooseCharge{"cfChooseCharge", 0, "choose what charge to keep. 0 = +1 and -1, 1 = +1, -1 = -1"}; // *) Define and initialize all data members to be called in the main process* functions: // **) Task configuration: struct TaskConfiguration { float fSigmaInel = 7.71; float fEtaGap = 1.0; - float fCorrCentCutRatio = 0.5; + float fCorrCentCutPercentage = 50; std::vector fEventCutSwitch; + std::vector fEventCutOrdering; std::vector fVertexZCutRange = {-10., 10.}; std::vector fCentCutRange = {0., 80.}; std::vector fCorrCent = {0, 2}; std::vector fMultNContribBound; + std::vector fParticleCutSwitch; + std::vector fParticleCutOrdering; std::vector fPtCutRange = {0.2, 5.}; std::vector fEtaCutRange = {-0.8, 0.8}; - std::vector fTPCNumClusterCutRange = {70, 160}; + std::vector fDCAxyCutRange = {-3.2, 3.2}; std::vector fDCAzCutRange = {-2.4, 2.4}; @@ -320,7 +334,10 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to int fOccupancyEstimator = 0; int fNContribCutInf = 2; int fOccupancyCutSup = 6000; + int fTPCNumClusterCutInf = 70; + int fITSNClsInf = 5; int fRNGSeed = 0; + int fChooseCharge = 0; bool fDryRun = false; // book all histos and run without filling and calculating anything bool fExternalFileSwitch = false; @@ -683,88 +700,90 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to } } bool doNotPassCut = false; - for (int i = 0; i < EEventCuts_N; ++i) { - switch (i) { - default: - doNotPassCut = false; - break; - case eSel8Cut: - doNotPassCut = tc.fEventCutSwitch[eSel8Cut].starts_with("1") && !collision.sel8(); - break; - case eNoCollInTimeRangeStandard: - doNotPassCut = tc.fEventCutSwitch[eNoCollInTimeRangeStandard].starts_with("1") && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard); - break; - case eNoCollInRofStandard: - doNotPassCut = tc.fEventCutSwitch[eNoCollInRofStandard].starts_with("1") && !collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard); - break; - case eNoSameBunchPileUp: - doNotPassCut = tc.fEventCutSwitch[eNoSameBunchPileUp].starts_with("1") && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup); - break; - case eIsVertexITSTPC: - doNotPassCut = tc.fEventCutSwitch[eIsVertexITSTPC].starts_with("1") && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC); - break; - case eIsGoodITSLayersAll: - doNotPassCut = tc.fEventCutSwitch[eIsGoodITSLayersAll].starts_with("1") && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll); - break; - case eIsGoodZvtxFT0vsPV: - doNotPassCut = tc.fEventCutSwitch[eIsGoodZvtxFT0vsPV].starts_with("1") && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV); - break; - case eNoHighMultCollInPrevRof: - doNotPassCut = tc.fEventCutSwitch[eNoHighMultCollInPrevRof].starts_with("1") && !collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof); - break; - case eNContribCut: - doNotPassCut = tc.fEventCutSwitch[eNContribCut].starts_with("1") && (collision.numContrib() < tc.fNContribCutInf); - break; - case eOccupancyCut: { - auto thisOccu = chooseOccu(collision, tc.fOccupancyEstimator); - doNotPassCut = tc.fEventCutSwitch[eOccupancyCut].starts_with("1") && (thisOccu > tc.fOccupancyCutSup); - break; - } - case eVertexZCut: - if constexpr ((rs == eRec || rs == eRecAndSim) && rm == eReal) { - doNotPassCut = tc.fEventCutSwitch[eVertexZCut].starts_with("1") && (collision.posZ() > tc.fVertexZCutRange[1] || collision.posZ() < tc.fVertexZCutRange[0]); - } - if constexpr ((rs == eSim || rs == eRecAndSim) && rm == eMC) { - doNotPassCut = tc.fEventCutSwitch[eVertexZCut].starts_with("1") && (collision.mcCollision().posZ() > tc.fVertexZCutRange[1] || collision.mcCollision().posZ() < tc.fVertexZCutRange[0]); + for (int i = 0; i < eEventCuts_N; ++i) { + if (tc.fEventCutSwitch[i].starts_with("1")) { + switch (tc.fEventCutOrdering[i]) { + default: + doNotPassCut = false; + break; + case eSel8Cut: + doNotPassCut = !collision.sel8(); + break; + case eNoCollInTimeRangeStandard: + doNotPassCut = !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard); + break; + case eNoCollInRofStandard: + doNotPassCut = !collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard); + break; + case eNoSameBunchPileUp: + doNotPassCut = !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup); + break; + case eIsVertexITSTPC: + doNotPassCut = !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC); + break; + case eIsGoodITSLayersAll: + doNotPassCut = !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll); + break; + case eIsGoodZvtxFT0vsPV: + doNotPassCut = !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV); + break; + case eNoHighMultCollInPrevRof: + doNotPassCut = !collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof); + break; + case eNContribCut: + doNotPassCut = (collision.numContrib() < tc.fNContribCutInf); + break; + case eOccupancyCut: { + auto thisOccu = chooseOccu(collision, tc.fOccupancyEstimator); + doNotPassCut = (thisOccu > tc.fOccupancyCutSup); + break; } - break; - case eCentCut: - if constexpr ((rs == eRec || rs == eRecAndSim) && rm == eReal) { - auto thisCent = chooseCent(collision, tc.fCentralityEstimator); - doNotPassCut = tc.fEventCutSwitch[eCentCut].starts_with("1") && (thisCent > tc.fCentCutRange[1] || thisCent < tc.fCentCutRange[0]); + case eVertexZCut: + if constexpr ((rs == eRec || rs == eRecAndSim) && rm == eReal) { + doNotPassCut = (collision.posZ() > tc.fVertexZCutRange[1] || collision.posZ() < tc.fVertexZCutRange[0]); + } + if constexpr ((rs == eSim || rs == eRecAndSim) && rm == eMC) { + doNotPassCut = (collision.mcCollision().posZ() > tc.fVertexZCutRange[1] || collision.mcCollision().posZ() < tc.fVertexZCutRange[0]); + } + break; + case eCentCut: + if constexpr ((rs == eRec || rs == eRecAndSim) && rm == eReal) { + auto thisCent = chooseCent(collision, tc.fCentralityEstimator); + doNotPassCut = (thisCent > tc.fCentCutRange[1] || thisCent < tc.fCentCutRange[0]); + } + if constexpr ((rs == eSim || rs == eRecAndSim) && rm == eMC) { + auto impactParameter = collision.mcCollision().impactParameter(); + auto centralityMC = math::PI * impactParameter * impactParameter / tc.fSigmaInel; + doNotPassCut = (centralityMC > tc.fCentCutRange[1] || centralityMC < tc.fCentCutRange[0]); + } + break; + case eCorrCentCut: { + auto cent1 = chooseCent(collision, tc.fCorrCent[0]); + auto cent2 = chooseCent(collision, tc.fCorrCent[1]); + doNotPassCut = (std::fabsf(cent1 - cent2) > tc.fCorrCentCutPercentage); + break; } - if constexpr ((rs == eSim || rs == eRecAndSim) && rm == eMC) { - auto impactParameter = collision.mcCollision().impactParameter(); - auto centralityMC = math::PI * impactParameter * impactParameter / tc.fSigmaInel; - doNotPassCut = tc.fEventCutSwitch[eCentCut].starts_with("1") && (centralityMC > tc.fCentCutRange[1] || centralityMC < tc.fCentCutRange[0]); + case eMultNContribCut: { + auto thisRefMult = chooseMult(collision, tc.fMultiplicityTables); + auto thisNContrib = collision.numContrib(); + auto upper = geoCuts.fMultNContribBounds[eUpper]->Eval(thisRefMult); + auto lower = geoCuts.fMultNContribBounds[eLower]->Eval(thisRefMult); + doNotPassCut = (thisNContrib < lower || thisNContrib > upper); + break; } - break; - case eCorrCentCut: { - auto cent1 = chooseCent(collision, tc.fCorrCent[0]); - auto cent2 = chooseCent(collision, tc.fCorrCent[1]); - doNotPassCut = tc.fEventCutSwitch[eCorrCentCut].starts_with("1") && (2 * std::fabsf(cent1 - cent2) / (cent1 + cent2) > tc.fCorrCentCutRatio); - break; - } - case eMultNContribCut: { - auto thisRefMult = chooseMult(collision, tc.fMultiplicityTables); - auto thisNContrib = collision.numContrib(); - auto upper = geoCuts.fMultNContribBounds[eUpper]->Eval(thisRefMult); - auto lower = geoCuts.fMultNContribBounds[eLower]->Eval(thisRefMult); - doNotPassCut = tc.fEventCutSwitch[eMultNContribCut].starts_with("1") && (thisNContrib < lower || thisNContrib > upper); - break; - } - } // end of switch (i) { + } // end of switch (i) { - if (doNotPassCut) { - if constexpr (rm == eReal) { - misc.fMiscEventCutsCounter[eRec]->Fill(i + 0.5, 1); - } - if constexpr (rm == eMC) { - misc.fMiscEventCutsCounter[eSim]->Fill(i + 0.5, 1); + if (doNotPassCut) { + if constexpr (rm == eReal) { + misc.fMiscEventCutsCounter[eRec]->Fill(i + 0.5, 1); + } + if constexpr (rm == eMC) { + misc.fMiscEventCutsCounter[eSim]->Fill(i + 0.5, 1); + } + return false; } - return false; } - } // end of for (int i = 0; i < EEventCuts_N; ++i) { + } // end of for (int i = 0; i < eEventCuts_N; ++i) { return true; } // end of bool eventCuts(T1 const& collision) @@ -830,85 +849,89 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to template bool particleCuts(T const& track) { + if constexpr ((rs == eSim || rs == eRecAndSim) && rm == eMC) { + if (!track.has_mcParticle()) { + return false; + } + } bool doNotPassCut = false; - for (int i = 0; i < EParticleCuts_N; ++i) { - if constexpr ((rs == eRec || rs == eRecAndSim) && rm == eReal) { - switch (i) { + for (int i = 0; i < eParticleCuts_N; ++i) { + if (tc.fParticleCutSwitch[i].starts_with("1")) { + switch (tc.fParticleCutOrdering[i]) { default: doNotPassCut = false; break; case ePtCut: - doNotPassCut = tc.fParticleCutSwitch[ePtCut].starts_with("1") && (track.pt() < tc.fPtCutRange[0] || track.pt() > tc.fPtCutRange[1]); - break; - case eEtaCut: - doNotPassCut = tc.fParticleCutSwitch[eEtaCut].starts_with("1") && (track.eta() < tc.fEtaCutRange[0] || track.eta() > tc.fEtaCutRange[1]); - break; - case eChargeCut: - doNotPassCut = tc.fParticleCutSwitch[eChargeCut].starts_with("1") && (track.sign() != 1 && track.sign() != -1); - break; - case eTPCNumClusterCut: - doNotPassCut = tc.fParticleCutSwitch[eTPCNumClusterCut].starts_with("1") && (track.tpcNClsFound() < tc.fTPCNumClusterCutRange[0] || track.tpcNClsFound() > tc.fTPCNumClusterCutRange[1]); - break; - case eDCAxyCut: - doNotPassCut = tc.fParticleCutSwitch[eDCAxyCut].starts_with("1") && (track.dcaXY() < tc.fDCAxyCutRange[0] || track.dcaXY() > tc.fDCAxyCutRange[1]); - break; - case eDCAzCut: - doNotPassCut = tc.fParticleCutSwitch[eDCAzCut].starts_with("1") && (track.dcaZ() < tc.fDCAzCutRange[0] || track.dcaZ() > tc.fDCAzCutRange[1]); - break; - } - } // end of if constexpr ((rs == eRec || rs == eRecAndSim) && rm == eReal) { - - if constexpr ((rs == eSim || rs == eRecAndSim) && rm == eMC) { - if (!track.has_mcParticle()) { - return false; - } - auto thisMCParticle = track.mcParticle(); - switch (i) { - default: - doNotPassCut = false; - break; - case ePtCut: - doNotPassCut = tc.fParticleCutSwitch[ePtCut].starts_with("1") && (thisMCParticle.pt() < tc.fPtCutRange[0] || thisMCParticle.pt() > tc.fPtCutRange[1]); + if constexpr ((rs == eRec || rs == eRecAndSim) && rm == eReal) { + doNotPassCut = (track.pt() < tc.fPtCutRange[0] || track.pt() > tc.fPtCutRange[1]); + } + if constexpr ((rs == eSim || rs == eRecAndSim) && rm == eMC) { + auto thisMCParticle = track.mcParticle(); + doNotPassCut = (thisMCParticle.pt() < tc.fPtCutRange[0] || thisMCParticle.pt() > tc.fPtCutRange[1]); + } break; case eEtaCut: - doNotPassCut = tc.fParticleCutSwitch[eEtaCut].starts_with("1") && (thisMCParticle.eta() < tc.fEtaCutRange[0] || thisMCParticle.eta() > tc.fEtaCutRange[1]); + if constexpr ((rs == eRec || rs == eRecAndSim) && rm == eReal) { + doNotPassCut = (track.eta() < tc.fEtaCutRange[0] || track.eta() > tc.fEtaCutRange[1]); + } + if constexpr ((rs == eSim || rs == eRecAndSim) && rm == eMC) { + auto thisMCParticle = track.mcParticle(); + doNotPassCut = (thisMCParticle.eta() < tc.fEtaCutRange[0] || thisMCParticle.eta() > tc.fEtaCutRange[1]); + } break; case eChargeCut: - if (tc.fParticleCutSwitch[eChargeCut].starts_with("1")) { + if constexpr ((rs == eRec || rs == eRecAndSim) && rm == eReal) { + if (tc.fChooseCharge == 0) { + doNotPassCut = (track.sign() != 1 && track.sign() != -1); + } else if (tc.fChooseCharge == 1) { + doNotPassCut = (track.sign() != 1); + } else if (tc.fChooseCharge == -1) { + doNotPassCut = (track.sign() != -1); + } + } + if constexpr ((rs == eSim || rs == eRecAndSim) && rm == eMC) { + auto thisMCParticle = track.mcParticle(); TParticlePDG* particlePDG = pdg->GetParticle(thisMCParticle.pdgCode()); if (particlePDG) { const int chargeMCUnit = 3; auto chargeMC = particlePDG->Charge() / chargeMCUnit; - if (chargeMC != 1 && chargeMC != -1) { - doNotPassCut = true; + if (tc.fChooseCharge == 0) { + doNotPassCut = (chargeMC != 1 && chargeMC != -1); + } else if (tc.fChooseCharge == 1) { + doNotPassCut = (chargeMC != 1); + } else if (tc.fChooseCharge == -1) { + doNotPassCut = (chargeMC != -1); } } else { doNotPassCut = true; } - } // end of if (tc.fParticleCutSwitch[eChargeCut]) { + } break; case eTPCNumClusterCut: - doNotPassCut = tc.fParticleCutSwitch[eTPCNumClusterCut].starts_with("1") && (track.tpcNClsFound() < tc.fTPCNumClusterCutRange[0] || track.tpcNClsFound() > tc.fTPCNumClusterCutRange[1]); + doNotPassCut = (track.tpcNClsFound() < tc.fTPCNumClusterCutInf); break; case eDCAxyCut: - doNotPassCut = tc.fParticleCutSwitch[eDCAxyCut].starts_with("1") && (track.dcaXY() < tc.fDCAxyCutRange[0] || track.dcaXY() > tc.fDCAxyCutRange[1]); + doNotPassCut = (track.dcaXY() < tc.fDCAxyCutRange[0] || track.dcaXY() > tc.fDCAxyCutRange[1]); break; case eDCAzCut: - doNotPassCut = tc.fParticleCutSwitch[eDCAzCut].starts_with("1") && (track.dcaZ() < tc.fDCAzCutRange[0] || track.dcaZ() > tc.fDCAzCutRange[1]); + doNotPassCut = (track.dcaZ() < tc.fDCAzCutRange[0] || track.dcaZ() > tc.fDCAzCutRange[1]); + break; + case eITSNClsCut: + doNotPassCut = (track.itsNCls() < tc.fITSNClsInf); break; } - } // end of if constexpr ((rs == eSim || rs == eRecAndSim) && rm == eMC) { - if (doNotPassCut) { - if constexpr (rm == eReal) { - misc.fMiscParticleCutsCounter[eRec]->Fill(i + 0.5, 1); - } - if constexpr (rm == eMC) { - misc.fMiscParticleCutsCounter[eSim]->Fill(i + 0.5, 1); + if (doNotPassCut) { + if constexpr (rm == eReal) { + misc.fMiscParticleCutsCounter[eRec]->Fill(i + 0.5, 1); + } + if constexpr (rm == eMC) { + misc.fMiscParticleCutsCounter[eSim]->Fill(i + 0.5, 1); + } + return false; } - return false; - } - } // end of for (int i = 0; i < EParticleCuts_N; ++i) { + } // end of if (tc.fParticleCutSwitch[i].starts_with("1")) { + } // end of for (int i = 0; i < eParticleCuts_N; ++i) { return true; } // end of bool particleCuts(T const& track) @@ -922,12 +945,20 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to pc.fParticleHist[eHistPhi][eRec][eBefore]->Fill(track.phi()); pc.fParticleHist[eHistEta][eRec][eBefore]->Fill(track.eta()); pc.fParticleHist[eHistCharge][eRec][eBefore]->Fill(track.sign()); + pc.fParticleHist[eHistTPCNumCluster][eRec][eBefore]->Fill(track.tpcNClsFound()); + pc.fParticleHist[eHistDCAxy][eRec][eBefore]->Fill(track.dcaXY()); + pc.fParticleHist[eHistDCAz][eRec][eBefore]->Fill(track.dcaZ()); + pc.fParticleHist[eHistITSNCls][eRec][eBefore]->Fill(track.itsNCls()); } if constexpr (cuts == eAfter) { pc.fParticleHist[eHistPt][eRec][eAfter]->Fill(track.pt()); pc.fParticleHist[eHistPhi][eRec][eAfter]->Fill(track.phi()); pc.fParticleHist[eHistEta][eRec][eAfter]->Fill(track.eta()); pc.fParticleHist[eHistCharge][eRec][eAfter]->Fill(track.sign()); + pc.fParticleHist[eHistTPCNumCluster][eRec][eAfter]->Fill(track.tpcNClsFound()); + pc.fParticleHist[eHistDCAxy][eRec][eAfter]->Fill(track.dcaXY()); + pc.fParticleHist[eHistDCAz][eRec][eAfter]->Fill(track.dcaZ()); + pc.fParticleHist[eHistITSNCls][eRec][eAfter]->Fill(track.itsNCls()); } } @@ -1116,7 +1147,7 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to } } // end of if constexpr (rs == eRecAndSim) { } // end of for (auto const& track : tracks) { - if (tc.fMasterCutSwitch && passEventCutsReal && tc.fQualityAssuranceSwitch) { + if (tc.fMasterCutSwitch && passEventCutsReal && tc.fQualityAssuranceSwitch && (thisCollCentReal > tc.fQASelTrackCentRange[0] && thisCollCentReal < tc.fQASelTrackCentRange[1])) { qa.fHistSelectedTrack[eAfter]->Fill(nPassedTrack); } @@ -1210,13 +1241,14 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to tc.fNContribCutInf = cfNContribCutInf; tc.fOccupancyCutSup = cfOccupancyCutSup; tc.fCorrCent = cfCorrCent.value; - tc.fCorrCentCutRatio = cfCorrCentCutRatio; + tc.fCorrCentCutPercentage = cfCorrCentCutPercentage; tc.fMultNContribBound = cfMultNContribBound.value; tc.fParticleCutSwitch = cfParticleCutSwitch.value; tc.fPtCutRange = cfPtCutRange.value; tc.fEtaCutRange = cfEtaCutRange.value; - tc.fTPCNumClusterCutRange = cfTPCNumClusterCutRange.value; + tc.fTPCNumClusterCutInf = cfTPCNumClusterCutInf; + tc.fITSNClsInf = cfITSNClsInf; tc.fDCAxyCutRange = cfDCAxyCutRange.value; tc.fDCAzCutRange = cfDCAzCutRange.value; @@ -1228,6 +1260,57 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to tc.fInternalValidationSwitch = cfInternalValidationSwitch; tc.fIVHarmonics = cfIVHarmonics.value; tc.fQASelTrackCentRange = cfQASelTrackCentRange.value; + tc.fChooseCharge = cfChooseCharge; + + for (const auto& cutSwitch : tc.fEventCutSwitch) { + if (cutSwitch.ends_with("Sel8")) { + tc.fEventCutOrdering.push_back(static_cast(eSel8Cut)); + } else if (cutSwitch.ends_with("IsGoodITSLayersAll")) { + tc.fEventCutOrdering.push_back(static_cast(eIsGoodITSLayersAll)); + } else if (cutSwitch.ends_with("NoCollInTimeRangeStandard")) { + tc.fEventCutOrdering.push_back(static_cast(eNoCollInTimeRangeStandard)); + } else if (cutSwitch.ends_with("NoCollInRofStandard")) { + tc.fEventCutOrdering.push_back(static_cast(eNoCollInRofStandard)); + } else if (cutSwitch.ends_with("IsGoodZvtxFT0vsPV")) { + tc.fEventCutOrdering.push_back(static_cast(eIsGoodZvtxFT0vsPV)); + } else if (cutSwitch.ends_with("NoHighMultCollInPrevRof")) { + tc.fEventCutOrdering.push_back(static_cast(eNoHighMultCollInPrevRof)); + } else if (cutSwitch.ends_with("NoSameBunchPileUp")) { + tc.fEventCutOrdering.push_back(static_cast(eNoSameBunchPileUp)); + } else if (cutSwitch.ends_with("IsVertexITSTPC")) { + tc.fEventCutOrdering.push_back(static_cast(eIsVertexITSTPC)); + } else if (cutSwitch.ends_with("Cent")) { + tc.fEventCutOrdering.push_back(static_cast(eCentCut)); + } else if (cutSwitch.ends_with("VertexZ")) { + tc.fEventCutOrdering.push_back(static_cast(eVertexZCut)); + } else if (cutSwitch.ends_with("NContrib")) { + tc.fEventCutOrdering.push_back(static_cast(eNContribCut)); + } else if (cutSwitch.ends_with("Occupancy")) { + tc.fEventCutOrdering.push_back(static_cast(eOccupancyCut)); + } else if (cutSwitch.ends_with("CorrCent")) { + tc.fEventCutOrdering.push_back(static_cast(eCorrCentCut)); + } else if (cutSwitch.ends_with("MultNContribCut")) { + tc.fEventCutOrdering.push_back(static_cast(eMultNContribCut)); + } + } + + for (const auto& cutSwitch : tc.fParticleCutSwitch) { + if (cutSwitch.ends_with("Pt")) { + tc.fParticleCutOrdering.push_back(static_cast(ePtCut)); + } else if (cutSwitch.ends_with("Charge")) { + tc.fParticleCutOrdering.push_back(static_cast(eChargeCut)); + } else if (cutSwitch.ends_with("Eta")) { + tc.fParticleCutOrdering.push_back(static_cast(eEtaCut)); + } else if (cutSwitch.ends_with("TPCNumCluster")) { + tc.fParticleCutOrdering.push_back(static_cast(eTPCNumClusterCut)); + } else if (cutSwitch.ends_with("DCAxy")) { + tc.fParticleCutOrdering.push_back(static_cast(eDCAxyCut)); + } else if (cutSwitch.ends_with("DCAz")) { + tc.fParticleCutOrdering.push_back(static_cast(eDCAzCut)); + } else if (cutSwitch.ends_with("ITSNCls")) { + tc.fParticleCutOrdering.push_back(static_cast(eITSNClsCut)); + } + } // *) Book base list: auto* temp = new TList(); @@ -1255,32 +1338,32 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to misc.fMiscHistList->Add(misc.fMiscHistRunNumber); if (doprocessRec || doprocessRecSim) { - misc.fMiscEventCutsCounter[eRec] = new TH1D("[eEventCutsCounter][eRec]", "Event Level Cuts Counter for Reconstructed Data (Sequential)", EEventCuts_N, 0, EEventCuts_N); + misc.fMiscEventCutsCounter[eRec] = new TH1D("[eEventCutsCounter][eRec]", "Event Level Cuts Counter for Reconstructed Data (Sequential)", eEventCuts_N, 0, eEventCuts_N); misc.fMiscEventCutsCounter[eRec]->SetStats(false); - for (int i = 1; i <= EEventCuts_N; ++i) { + for (int i = 1; i <= eEventCuts_N; ++i) { misc.fMiscEventCutsCounter[eRec]->GetXaxis()->SetBinLabel(i, tc.fEventCutSwitch[i - 1].c_str()); } misc.fMiscHistList->Add(misc.fMiscEventCutsCounter[eRec]); - misc.fMiscParticleCutsCounter[eRec] = new TH1D("[eParticleCutsCounter][eRec]", "Particle Level Cuts Counter for Reconstructed Data (Sequential)", EParticleCuts_N, 0, EParticleCuts_N); + misc.fMiscParticleCutsCounter[eRec] = new TH1D("[eParticleCutsCounter][eRec]", "Particle Level Cuts Counter for Reconstructed Data (Sequential)", eParticleCuts_N, 0, eParticleCuts_N); misc.fMiscParticleCutsCounter[eRec]->SetStats(false); - for (int i = 1; i <= EParticleCuts_N; ++i) { + for (int i = 1; i <= eParticleCuts_N; ++i) { misc.fMiscParticleCutsCounter[eRec]->GetXaxis()->SetBinLabel(i, tc.fParticleCutSwitch[i - 1].c_str()); } misc.fMiscHistList->Add(misc.fMiscParticleCutsCounter[eRec]); } if (doprocessSim || doprocessRecSim) { - misc.fMiscEventCutsCounter[eSim] = new TH1D("[eEventCutsCounter][eSim]", "Event Level Cuts Counter for Simulated Data (Sequential)", EEventCuts_N, 0, EEventCuts_N); + misc.fMiscEventCutsCounter[eSim] = new TH1D("[eEventCutsCounter][eSim]", "Event Level Cuts Counter for Simulated Data (Sequential)", eEventCuts_N, 0, eEventCuts_N); misc.fMiscEventCutsCounter[eSim]->SetStats(false); - for (int i = 1; i <= EEventCuts_N; ++i) { + for (int i = 1; i <= eEventCuts_N; ++i) { misc.fMiscEventCutsCounter[eSim]->GetXaxis()->SetBinLabel(i, tc.fEventCutSwitch[i - 1].c_str()); } misc.fMiscHistList->Add(misc.fMiscEventCutsCounter[eSim]); - misc.fMiscParticleCutsCounter[eSim] = new TH1D("[eParticleCutsCounter][eSim]", "Particle Level Cuts Counter for Simulated Data (Sequential)", EParticleCuts_N, 0, EParticleCuts_N); + misc.fMiscParticleCutsCounter[eSim] = new TH1D("[eParticleCutsCounter][eSim]", "Particle Level Cuts Counter for Simulated Data (Sequential)", eParticleCuts_N, 0, eParticleCuts_N); misc.fMiscParticleCutsCounter[eSim]->SetStats(false); - for (int i = 1; i <= EParticleCuts_N; ++i) { + for (int i = 1; i <= eParticleCuts_N; ++i) { misc.fMiscParticleCutsCounter[eSim]->GetXaxis()->SetBinLabel(i, tc.fParticleCutSwitch[i - 1].c_str()); } misc.fMiscHistList->Add(misc.fMiscParticleCutsCounter[eSim]); @@ -1293,11 +1376,14 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to fBaseList->Add(pc.fParticleHistList); // any nested TList in the base TList appears as a subdir in the output ROOT file std::vector> lParticleHistBins(eParticleHistograms_N); - lParticleHistBins[eHistPt] = cfPtBins.value; lParticleHistBins[eHistPhi] = cfPhiBins.value; lParticleHistBins[eHistEta] = cfEtaBins.value; lParticleHistBins[eHistCharge] = {5, -2.5, 2.5}; + lParticleHistBins[eHistTPCNumCluster] = cfTPCNumClusterBins.value; + lParticleHistBins[eHistDCAxy] = cfDCAxyBins.value; + lParticleHistBins[eHistDCAz] = cfDCAzBins.value; + lParticleHistBins[eHistITSNCls] = cfITSNClsBins.value; if (doprocessRec || doprocessRecSim) { for (int i = 0; i < eParticleHistograms_N; ++i) { @@ -1316,6 +1402,9 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to if (doprocessSim || doprocessRecSim) { for (int i = 0; i < eParticleHistograms_N; ++i) { + if (i == eHistTPCNumCluster || i == eHistDCAxy || i == eHistDCAz || i == eHistITSNCls) { + continue; + } const auto& bins = lParticleHistBins[i]; for (int ba = 0; ba < eBeforeAfterCuts_N; ++ba) { if (!tc.fMasterCutSwitch && ba == eAfter) { @@ -1336,7 +1425,6 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to fBaseList->Add(ec.fEventHistList); std::vector> lEventHistBins(eEventHistograms_N); - lEventHistBins[eHistCentrality] = cfCentBins.value; lEventHistBins[eHistMult] = cfMultBins.value; lEventHistBins[eHistRefMult] = cfRefMultBins.value; @@ -1367,8 +1455,8 @@ struct MultiparticleCorrelationsMei // this name is used in lower-case format to ec.fEventHist[eHistCentrality][eRec][eAfter]->GetXaxis()->SetTitle(Form("%s (%s)", EventHistTitles[eHistCentrality], CentralityEstimatorNames[tc.fCentralityEstimator])); ec.fEventHist[eHistRefMult][eRec][eBefore]->GetXaxis()->SetTitle(Form("%s (%s)", EventHistTitles[eHistRefMult], MultiplicityTablesNames[tc.fMultiplicityTables])); ec.fEventHist[eHistRefMult][eRec][eAfter]->GetXaxis()->SetTitle(Form("%s (%s)", EventHistTitles[eHistRefMult], MultiplicityTablesNames[tc.fMultiplicityTables])); - ec.fEventHist[eHistOccupancy][eRec][eBefore]->GetXaxis()->SetTitle(Form("%s (%s)", EventHistTitles[eHistRefMult], OccupancyEstimatorNames[tc.fOccupancyEstimator])); - ec.fEventHist[eHistOccupancy][eRec][eAfter]->GetXaxis()->SetTitle(Form("%s (%s)", EventHistTitles[eHistRefMult], OccupancyEstimatorNames[tc.fOccupancyEstimator])); + ec.fEventHist[eHistOccupancy][eRec][eBefore]->GetXaxis()->SetTitle(Form("%s (%s)", EventHistTitles[eHistOccupancy], OccupancyEstimatorNames[tc.fOccupancyEstimator])); + ec.fEventHist[eHistOccupancy][eRec][eAfter]->GetXaxis()->SetTitle(Form("%s (%s)", EventHistTitles[eHistOccupancy], OccupancyEstimatorNames[tc.fOccupancyEstimator])); } // end of if (doprocessRec || doprocessRecSim) { if (doprocessSim || doprocessRecSim) {