diff --git a/PWGLF/Tasks/Nuspex/piKpRAA.cxx b/PWGLF/Tasks/Nuspex/piKpRAA.cxx index 33e06b49d56..9d046f974db 100644 --- a/PWGLF/Tasks/Nuspex/piKpRAA.cxx +++ b/PWGLF/Tasks/Nuspex/piKpRAA.cxx @@ -81,10 +81,10 @@ using TracksMC = soa::Join, kNEtaHists> dEdxPiV0{}; -std::array, kNEtaHists> dEdxPrV0{}; -std::array, kNEtaHists> dEdxElV0{}; -std::array, kNEtaHists> dEdxPiTOF{}; +std::array, kNEtaHists> dEdxPiV0{}; +std::array, kNEtaHists> dEdxPrV0{}; +std::array, kNEtaHists> dEdxElV0{}; +std::array, kNEtaHists> dEdxPiTOF{}; std::array, kNEtaHists> dEdx{}; std::array, kNEtaHists> nClVsdEdxPiV0{}; std::array, kNEtaHists> nClVsdEdxElV0{}; @@ -131,7 +131,6 @@ struct PiKpRAA { static constexpr float kLowEta[kNEtaHists] = {-0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6}; static constexpr float kHighEta[kNEtaHists] = {-0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8}; // static constexpr float kLowEta[kNEtaHists] = {0.0, 0.2, 0.4, 0.6}; - // static constexpr float kHighEta[kNEtaHists] = {0.2, 0.4, 0.6, 0.8}; static constexpr float DefaultLifetimeCuts[1][2] = {{30., 20.}}; Configurable> lifetimecut{"lifetimecut", {DefaultLifetimeCuts[0], 2, {"lifetimecutLambda", "lifetimecutK0S"}}, "lifetimecut"}; @@ -219,7 +218,6 @@ struct PiKpRAA { Configurable selINELgt0{"selINELgt0", true, "Select INEL > 0?"}; Configurable isApplyFT0CbasedOccupancy{"isApplyFT0CbasedOccupancy", false, "T0C Occu cut"}; Configurable applyNchSel{"applyNchSel", false, "Use mid-rapidity-based Nch selection"}; - Configurable skipRecoColGTOne{"skipRecoColGTOne", true, "Remove collisions if reconstructed more than once"}; // Event selection Configurable posZcut{"posZcut", +10.0, "z-vertex position cut"}; @@ -228,6 +226,7 @@ struct PiKpRAA { Configurable minOccCut{"minOccCut", 0., "min Occu cut"}; Configurable maxOccCut{"maxOccCut", 500., "max Occu cut"}; Configurable nSigmaNchCut{"nSigmaNchCut", 3., "nSigma Nch selection"}; + Configurable tpcNchAcceptance{"tpcNchAcceptance", 0.5, "Eta window to measure Nch MC for Nch vs Cent distribution"}; ConfigurableAxis binsPtPhiCut{"binsPtPhiCut", {VARIABLE_WIDTH, 0.0, 0.6, 0.8, 1.0, 1.4, 1.8, 2.2, 2.6, 3.0, 3.5, 4.0, 5.0, 7.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 45.0, 50.0}, "pT"}; ConfigurableAxis binsPtV0s{"binsPtV0s", {VARIABLE_WIDTH, 0, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.5, 3.0, 3.5, 4, 5, 7, 9, 12, 15, 20}, "pT"}; @@ -489,10 +488,10 @@ struct PiKpRAA { for (int i = 0; i < kNEtaHists; ++i) { dEdx[i] = registry.add(Form("dEdx_%s", endingEta[i]), Form("%s;Momentum (GeV/#it{c});dE/dx;", latexEta[i]), kTH3F, {axisPt, axisdEdx, axisCent}); pTVsP[i] = registry.add(Form("pTVsP_%s", endingEta[i]), Form("%s;Momentum (GeV/#it{c});#it{p}_{T} (GeV/#it{c});", latexEta[i]), kTH2F, {axisPt, axisPt}); - dEdxPiV0[i] = registry.add(Form("dEdxPiV0_%s", endingEta[i]), Form("#pi^{+} + #pi^{-}, %s;Momentum (GeV/#it{c});dE/dx;", latexEta[i]), kTH3F, {axisPtV0s, axisdEdx, axisCent}); - dEdxPrV0[i] = registry.add(Form("dEdxPrV0_%s", endingEta[i]), Form("p + #bar{p}, %s;Momentum (GeV/#it{c});dE/dx;", latexEta[i]), kTH3F, {axisPtV0s, axisdEdx, axisCent}); - dEdxElV0[i] = registry.add(Form("dEdxElV0_%s", endingEta[i]), Form("e^{+} + e^{-}, %s;Momentum (GeV/#it{c});dE/dx;", latexEta[i]), kTH3F, {axisPtV0s, axisdEdx, axisCent}); - dEdxPiTOF[i] = registry.add(Form("dEdxPiTOF_%s", endingEta[i]), Form("#pi^{+} + #pi^{-}, %s;Momentum (GeV/#it{c});dE/dx;", latexEta[i]), kTH3F, {axisPtV0s, axisdEdx, axisCent}); + dEdxPiV0[i] = registry.add(Form("dEdxPiV0_%s", endingEta[i]), Form("#pi^{+} + #pi^{-}, %s;Momentum (GeV/#it{c});dE/dx;", latexEta[i]), kTH2F, {axisPtV0s, axisdEdx}); + dEdxPrV0[i] = registry.add(Form("dEdxPrV0_%s", endingEta[i]), Form("p + #bar{p}, %s;Momentum (GeV/#it{c});dE/dx;", latexEta[i]), kTH2F, {axisPtV0s, axisdEdx}); + dEdxElV0[i] = registry.add(Form("dEdxElV0_%s", endingEta[i]), Form("e^{+} + e^{-}, %s;Momentum (GeV/#it{c});dE/dx;", latexEta[i]), kTH2F, {axisPtV0s, axisdEdx}); + dEdxPiTOF[i] = registry.add(Form("dEdxPiTOF_%s", endingEta[i]), Form("#pi^{+} + #pi^{-}, %s;Momentum (GeV/#it{c});dE/dx;", latexEta[i]), kTH2F, {axisPtV0s, axisdEdx}); nClVsdEdxPiV0[i] = registry.add(Form("NclVsdEdxPiV0_%s", endingEta[i]), Form("%s;#it{N}_{cl} used for PID;dE/dx;", latexEta[i]), kTH2F, {axisNcl, axisdEdx}); nClVsdEdxElV0[i] = registry.add(Form("NclVsdEdxElV0_%s", endingEta[i]), Form("%s;#it{N}_{cl} used for PID;dE/dx;", latexEta[i]), kTH2F, {axisNcl, axisdEdx}); nClVsdEdxPrV0[i] = registry.add(Form("NclVsdEdxPrV0_%s", endingEta[i]), Form("%s;#it{N}_{cl} used for PID;dE/dx;", latexEta[i]), kTH2F, {axisNcl, axisdEdx}); @@ -525,7 +524,7 @@ struct PiKpRAA { xtrkSel->SetBinLabel(12, "Passed all"); } - if (doprocessMC || doprocessSim) { + if (doprocessSim) { registry.add("zPosMC", "Generated Events With at least One Rec. Collision + Sel. criteria;;Entries;", kTH1F, {axisZpos}); registry.add("dcaVsPtPiDec", "Secondary pions from decays;#it{p}_{T} (GeV/#it{c});DCA_{xy} (cm);Centrality Perc.;", kTH3F, {axisPt, axisDCAxy, axisCent}); registry.add("dcaVsPtPrDec", "Secondary protons from decays;#it{p}_{T} (GeV/#it{c});DCA_{xy} (cm);Centrality Perc.;", kTH3F, {axisPt, axisDCAxy, axisCent}); @@ -534,19 +533,6 @@ struct PiKpRAA { registry.add("NclVsPhip", Form("#LTNcl#GT used for PID;%s (GeV/#it{c});#varphi", titlePorPt.data()), kTProfile2D, {{{axisXPhiCut}, {350, 0.0, 0.35}}}); } - if (doprocessMC) { - - registry.add("EventCounterMC", "Event counter", kTH1F, {axisEvent}); - - registry.add("PtPiVsCent", "", kTH2F, {axisPt, axisCent}); - registry.add("PtKaVsCent", "", kTH2F, {axisPt, axisCent}); - registry.add("PtPrVsCent", "", kTH2F, {axisPt, axisCent}); - - registry.add("PtPiVsCentMC", "", kTH2F, {axisPt, axisCent}); - registry.add("PtKaVsCentMC", "", kTH2F, {axisPt, axisCent}); - registry.add("PtPrVsCentMC", "", kTH2F, {axisPt, axisCent}); - } - if (doprocessSim) { // MC events passing the TVX requirement @@ -562,16 +548,20 @@ struct PiKpRAA { registry.add("PtKaVsCent_WithRecoEvt", "Generated Events With at least One Rec. Collision + Sel. criteria;;;", kTH2F, {axisPt, axisCent}); registry.add("PtPrVsCent_WithRecoEvt", "Generated Events With at least One Rec. Collision + Sel. criteria;;;", kTH2F, {axisPt, axisCent}); + registry.add("PtGenPiVsCent_WithRecoEvt", "Generated Events With at least One Rec. Collision + Sel. criteria;;;", kTH2F, {axisPt, axisCent}); + registry.add("PtGenKaVsCent_WithRecoEvt", "Generated Events With at least One Rec. Collision + Sel. criteria;;;", kTH2F, {axisPt, axisCent}); + registry.add("PtGenPrVsCent_WithRecoEvt", "Generated Events With at least One Rec. Collision + Sel. criteria;;;", kTH2F, {axisPt, axisCent}); + // Needed to calculate the denominator of the Acceptance X Efficiency registry.add("PtPiVsCentMC_WithRecoEvt", "Generated Events With at least One Rec. Collision;;;", kTH2F, {axisPt, axisCent}); registry.add("PtKaVsCentMC_WithRecoEvt", "Generated Events With at least One Rec. Collision;;;", kTH2F, {axisPt, axisCent}); registry.add("PtPrVsCentMC_WithRecoEvt", "Generated Events With at least One Rec. Collision;;;", kTH2F, {axisPt, axisCent}); // Needed for the Gen. Nch to Centrality conversion - registry.add("NchMCVsCent", "Generated Nch v.s. Centrality (At least Once Rec. Coll. + Sel. criteria);;Gen. Nch", kTH2F, {{axisCent, {nBinsNch, minNch, maxNch}}}); + registry.add("NchMCVsCent", "Generated Nch v.s. Centrality (At least Once Rec. Coll. + Sel. criteria);;Gen. Nch MC (|#eta|<0.8)", kTH2F, {{axisCent, {nBinsNch, minNch, maxNch}}}); // Needed to measure Event Loss - registry.add("NchMC_WithRecoEvt", "Generated Nch of Evts With at least one Rec. Coll. + Sel. criteria;Gen. Nch MC;Entries", kTH1F, {{nBinsNch, minNch, maxNch}}); + registry.add("NchMC_WithRecoEvt", "Generated Nch of Evts With at least one Rec. Coll. + Sel. criteria;Gen. Nch MC (|#eta|<0.8);Entries", kTH1F, {{nBinsNch, minNch, maxNch}}); registry.add("NchMC_AllGen", "Generated Nch of All Gen. Evts.;Gen. Nch;Entries", kTH1F, {{nBinsNch, minNch, maxNch}}); // Needed to measure Event Splitting @@ -580,22 +570,22 @@ struct PiKpRAA { registry.add("Centrality_AllRecoEvt", "Generated Events Irrespective of the number of times it was reconstructed + Evt. Selections;;Entries", kTH1F, {axisCent}); // Needed to calculate the numerator of the Signal Loss correction - registry.add("PtPiVsNchMC_WithRecoEvt", "Generated Events With at least One Rec. Collision;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); - registry.add("PtKaVsNchMC_WithRecoEvt", "Generated Events With at least One Rec. Collision;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); - registry.add("PtPrVsNchMC_WithRecoEvt", "Generated Events With at least One Rec. Collision;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); + registry.add("PtPiVsNchMC_WithRecoEvt", "Generated Events With at least One Rec. Collision;;Gen. Nch (|#eta|<0.8);", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); + registry.add("PtKaVsNchMC_WithRecoEvt", "Generated Events With at least One Rec. Collision;;Gen. Nch (|#eta|<0.8);", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); + registry.add("PtPrVsNchMC_WithRecoEvt", "Generated Events With at least One Rec. Collision;;Gen. Nch (|#eta|<0.8);", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); // Needed to calculate the denominator of the Signal Loss correction - registry.add("PtPiVsNchMC_AllGen", "All Generated Events;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); - registry.add("PtKaVsNchMC_AllGen", "All Generated Events;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); - registry.add("PtPrVsNchMC_AllGen", "All Generated Events;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); + registry.add("PtPiVsNchMC_AllGen", "All Generated Events;;Gen. Nch (|#eta|<0.8);", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); + registry.add("PtKaVsNchMC_AllGen", "All Generated Events;;Gen. Nch (|#eta|<0.8);", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); + registry.add("PtPrVsNchMC_AllGen", "All Generated Events;;Gen. Nch (|#eta|<0.8);", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); - registry.add("MCclosure_PtMCPiVsNchMC", "All Generated Events 4 MC closure;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); - registry.add("MCclosure_PtMCKaVsNchMC", "All Generated Events 4 MC closure;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); - registry.add("MCclosure_PtMCPrVsNchMC", "All Generated Events 4 MC closure;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); + registry.add("MCclosure_PtMCPiVsNchMC", "All Generated Events 4 MC closure;;Gen. Nch (|#eta|<0.8);", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); + registry.add("MCclosure_PtMCKaVsNchMC", "All Generated Events 4 MC closure;;Gen. Nch (|#eta|<0.8);", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); + registry.add("MCclosure_PtMCPrVsNchMC", "All Generated Events 4 MC closure;;Gen. Nch (|#eta|<0.8);", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); - registry.add("MCclosure_PtPiVsNchMC", "Gen Evts With at least one Rec. Coll. + Sel. criteria 4 MC closure;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); - registry.add("MCclosure_PtKaVsNchMC", "Gen Evts With at least one Rec. Coll. + Sel. criteria 4 MC closure;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); - registry.add("MCclosure_PtPrVsNchMC", "Gen Evts With at least one Rec. Coll. + Sel. criteria 4 MC closure;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); + registry.add("MCclosure_PtPiVsNchMC", "Gen Evts With at least one Rec. Coll. + Sel. criteria 4 MC closure;;Gen. Nch (|#eta|<0.8);", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); + registry.add("MCclosure_PtKaVsNchMC", "Gen Evts With at least one Rec. Coll. + Sel. criteria 4 MC closure;;Gen. Nch (|#eta|<0.8);", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); + registry.add("MCclosure_PtPrVsNchMC", "Gen Evts With at least one Rec. Coll. + Sel. criteria 4 MC closure;;Gen. Nch (|#eta|<0.8);", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}}); } LOG(info) << "\trequireGoodRct=" << requireGoodRct.value; @@ -902,7 +892,7 @@ struct PiKpRAA { // registry.fill(HIST("TOFExpEl2TOF"), momentum, tExpElTOF / tTOF); if (std::abs((tExpPiTOF / tTOF) - kOne) < v0Selections.maxExpTOFPi) { - dEdxPiTOF[indexEta]->Fill(momentum, dedx, centrality); + dEdxPiTOF[indexEta]->Fill(momentum, dedx); } } } @@ -1038,8 +1028,8 @@ struct PiKpRAA { nClVsdEdxPiV0[posIndexEta]->Fill(posNcl, posTrkdEdx); nClVsdEdxpPiV0[posIndexEta]->Fill(posNcl, posTrkdEdx); nClVsPpPiV0[negIndexEta]->Fill(negPorPt, negNcl); - dEdxPiV0[posIndexEta]->Fill(posTrkP, posTrkdEdx, centrality); - dEdxPiV0[negIndexEta]->Fill(negTrkP, negTrkdEdx, centrality); + dEdxPiV0[posIndexEta]->Fill(posTrkP, posTrkdEdx); + dEdxPiV0[negIndexEta]->Fill(negTrkP, negTrkdEdx); if (posTrkP > kMinPMIP && posTrkP < kMaxPMIP && posTrkdEdx > kMindEdxMIP && posTrkdEdx < kMaxdEdxMIP) { registry.fill(HIST("dEdxVsEtaPiMIPV0"), posTrkEta, posTrkdEdx); @@ -1064,7 +1054,7 @@ struct PiKpRAA { nClVsPpPrV0[posIndexEta]->Fill(posPorPt, posNcl); nClVsdEdxPrV0[posIndexEta]->Fill(posNcl, posTrkdEdx); nClVsdEdxpPrV0[posIndexEta]->Fill(posNcl, posTrkdEdx); - dEdxPrV0[posIndexEta]->Fill(posTrkP, posTrkdEdx, centrality); + dEdxPrV0[posIndexEta]->Fill(posTrkP, posTrkdEdx); } if (negTrackCharge < kZero) { registry.fill(HIST("nSigPiFromL"), negTrkPt, negTrack.tpcNSigmaPi()); @@ -1100,7 +1090,7 @@ struct PiKpRAA { nClVsPpPrV0[negIndexEta]->Fill(negPorPt, negNcl); nClVsdEdxPrV0[negIndexEta]->Fill(negNcl, negTrkdEdx); nClVsdEdxpPrV0[negIndexEta]->Fill(negNcl, negTrkdEdx); - dEdxPrV0[negIndexEta]->Fill(negTrkP, negTrkdEdx, centrality); + dEdxPrV0[negIndexEta]->Fill(negTrkP, negTrkdEdx); } } @@ -1135,8 +1125,8 @@ struct PiKpRAA { registry.fill(HIST("dEdxVsEtaElMIPV0p"), posTrkEta, posTrkdEdx); registry.fill(HIST("dEdxVsEtaElMIPV0"), negTrkEta, negTrkdEdx); registry.fill(HIST("dEdxVsEtaElMIPV0p"), negTrkEta, negTrkdEdx); - dEdxElV0[posIndexEta]->Fill(posTrkP, posTrkdEdx, centrality); - dEdxElV0[negIndexEta]->Fill(negTrkP, negTrkdEdx, centrality); + dEdxElV0[posIndexEta]->Fill(posTrkP, posTrkdEdx); + dEdxElV0[negIndexEta]->Fill(negTrkP, negTrkdEdx); } } } // v0s @@ -1145,264 +1135,6 @@ struct PiKpRAA { Preslice perCollision = aod::track::collisionId; Service pdg; - void processMC(aod::McCollisions::iterator const& mccollision, soa::SmallGroups const& collisions, BCsRun3 const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::McParticles const& mcParticles, TracksMC const& tracksMC) - { - for (const auto& collision : collisions) { - // Event selection - if (!isEventSelected(collision)) { - continue; - } - // MC collision? - if (!collision.has_mcCollision()) { - continue; - } - - registry.fill(HIST("EventCounterMC"), EvCutLabel::All); - - if (std::fabs(mccollision.posZ()) > posZcut) - continue; - - registry.fill(HIST("zPos"), collision.posZ()); - registry.fill(HIST("zPosMC"), mccollision.posZ()); - registry.fill(HIST("EventCounterMC"), EvCutLabel::VtxZ); - - const auto& foundBC = collision.foundBC_as(); - uint64_t timeStamp{foundBC.timestamp()}; - const int magField{getMagneticField(timeStamp)}; - - if (v0Selections.applyPhiCut) { - const int nextRunNumber{foundBC.runNumber()}; - if (currentRunNumberPhiSel != nextRunNumber) { - loadPhiCutSelections(timeStamp); - currentRunNumberPhiSel = nextRunNumber; - LOG(info) << "\tcurrentRunNumberPhiSel= " << currentRunNumberPhiSel << " timeStamp = " << timeStamp; - } - - // return if phi cut objects are nullptr - if (!(phiCut.hPhiCutHigh && phiCut.hPhiCutLow)) - return; - } - - // Remove collisions if reconstructed more than once - if (skipRecoColGTOne && (collisions.size() > kOne)) - continue; - - const float centrality{isT0Ccent ? collision.centFT0C() : collision.centFT0M()}; - registry.fill(HIST("T0Ccent"), centrality); - - const auto& groupedTracks{tracksMC.sliceBy(perCollision, collision.globalIndex())}; - - // Track selection with Open DCAxy - for (const auto& track : groupedTracks) { - // Track Selection - if (track.eta() < v0Selections.minEtaDaughter || track.eta() > v0Selections.maxEtaDaughter) - continue; - - if (track.pt() < v0Selections.minPt || track.pt() > v0Selections.maxPt) - continue; - - if (!trkSelGlobalOpenDCAxy.IsSelected(track)) - continue; - - if (!track.has_mcParticle()) - continue; - - // Get the MC particle - const auto& particle{track.mcParticle()}; - auto charge{0.}; - auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (pdgParticle != nullptr) - charge = pdgParticle->Charge(); - else - continue; - - // Is it a charged particle? - if (std::abs(charge) < kMinCharge) - continue; - - float phiPrime{track.phi()}; - phiPrimeFunc(phiPrime, magField, charge); - - const float pOrPt{v0Selections.usePinPhiSelection ? track.p() : track.pt()}; - if (v0Selections.applyPhiCut) { - if (!passesPhiSelection(pOrPt, phiPrime)) - continue; - } - - const int16_t nclFound{track.tpcNClsFound()}; - const int16_t nclPID{track.tpcNClsPID()}; - const int16_t ncl = v0Selections.useNclsPID ? nclPID : nclFound; - if (v0Selections.applyNclSel && ncl < v0Selections.minNcl) - continue; - - bool isPrimary{false}; - bool isDecay{false}; - bool isMaterial{false}; - if (particle.isPhysicalPrimary()) - isPrimary = true; - else if (particle.getProcess() == TMCProcess::kPDecay) - isDecay = true; - else - isMaterial = true; - - bool isPi{false}; - bool isPr{false}; - if (particle.pdgCode() == PDG_t::kPiPlus || particle.pdgCode() == PDG_t::kPiMinus) - isPi = true; - else if (particle.pdgCode() == PDG_t::kProton || particle.pdgCode() == PDG_t::kProtonBar) - isPr = true; - else - continue; - - if (isPrimary && !isDecay && !isMaterial) { - if (isPi && !isPr) - registry.fill(HIST("dcaVsPtPi"), track.pt(), track.dcaXY(), centrality); - if (isPr && !isPi) - registry.fill(HIST("dcaVsPtPr"), track.pt(), track.dcaXY(), centrality); - } - - if (isDecay && !isPrimary && !isMaterial) { - if (isPi && !isPr) - registry.fill(HIST("dcaVsPtPiDec"), track.pt(), track.dcaXY(), centrality); - if (isPr && !isPi) - registry.fill(HIST("dcaVsPtPrDec"), track.pt(), track.dcaXY(), centrality); - } - - if (isMaterial && !isPrimary && !isDecay) { - if (isPi && !isPr) - registry.fill(HIST("dcaVsPtPiMat"), track.pt(), track.dcaXY(), centrality); - if (isPr && !isPi) - registry.fill(HIST("dcaVsPtPrMat"), track.pt(), track.dcaXY(), centrality); - } - } - - // Global track + DCAxy selections - for (const auto& track : groupedTracks) { - // Track Selection - if (track.eta() < v0Selections.minEtaDaughter || track.eta() > v0Selections.maxEtaDaughter) - continue; - - if (track.pt() < v0Selections.minPt || track.pt() > v0Selections.maxPt) - continue; - - if (!trkSelGlobal.IsSelected(track)) - continue; - - // Has MC particle? - if (!track.has_mcParticle()) - continue; - - // Get the MC particle - const auto& particle{track.mcParticle()}; - auto charge{0.}; - auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (pdgParticle != nullptr) - charge = pdgParticle->Charge(); - else - continue; - - // Is it a charged particle? - if (std::abs(charge) < kMinCharge) - continue; - - float phiPrime{track.phi()}; - phiPrimeFunc(phiPrime, magField, charge); - - const float pOrPt{v0Selections.usePinPhiSelection ? track.p() : track.pt()}; - if (v0Selections.applyPhiCut) { - if (!passesPhiSelection(pOrPt, phiPrime)) - continue; - } - - const int16_t nclFound{track.tpcNClsFound()}; - const int16_t nclPID{track.tpcNClsPID()}; - const int16_t ncl = v0Selections.useNclsPID ? nclPID : nclFound; - if (v0Selections.applyNclSel && ncl < v0Selections.minNcl) - continue; - - int indexEta{-999}; - const float eta{track.eta()}; - for (int i = 0; i < kNEtaHists; ++i) { - if (eta >= kLowEta[i] && eta < kHighEta[i]) { - indexEta = i; - break; - } - } - - if (indexEta < kZeroInt || indexEta > kSevenInt) - continue; - - registry.fill(HIST("NclVsPhip"), pOrPt, phiPrime, ncl); - registry.fill(HIST("NclVsEtaPID"), eta, ncl); - registry.fill(HIST("NclVsEtaPIDp"), eta, ncl); - - bool isPrimary{false}; - if (particle.isPhysicalPrimary()) - isPrimary = true; - - if (!isPrimary) - continue; - - bool isPi{false}; - bool isKa{false}; - bool isPr{false}; - - if (particle.pdgCode() == PDG_t::kPiPlus || particle.pdgCode() == PDG_t::kPiMinus) - isPi = true; - else if (particle.pdgCode() == PDG_t::kKPlus || particle.pdgCode() == PDG_t::kKMinus) - isKa = true; - else if (particle.pdgCode() == PDG_t::kProton || particle.pdgCode() == PDG_t::kProtonBar) - isPr = true; - else - continue; - - if (isPi && !isKa && !isPr) - registry.fill(HIST("PtPiVsCent"), track.pt(), centrality); - if (isKa && !isPi && !isPr) - registry.fill(HIST("PtKaVsCent"), track.pt(), centrality); - if (isPr && !isPi && !isKa) - registry.fill(HIST("PtPrVsCent"), track.pt(), centrality); - } - - // Generated MC - for (const auto& particle : mcParticles) { - if (particle.eta() < v0Selections.minEtaDaughter || particle.eta() > v0Selections.maxEtaDaughter) - continue; - - if (particle.pt() < v0Selections.minPt || particle.pt() > v0Selections.maxPt) - continue; - - auto charge{0.}; - // Get the MC particle - auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (pdgParticle != nullptr) - charge = pdgParticle->Charge(); - else - continue; - - // Is it a charged particle? - if (std::abs(charge) < kMinCharge) - continue; - - // Is it a primary particle? - bool isPrimary{true}; - if (!particle.isPhysicalPrimary()) - isPrimary = false; - - if (isPrimary) { - if (particle.pdgCode() == PDG_t::kPiPlus || particle.pdgCode() == PDG_t::kPiMinus) // pion - registry.fill(HIST("PtPiVsCentMC"), particle.pt(), centrality); - else if (particle.pdgCode() == PDG_t::kKPlus || particle.pdgCode() == PDG_t::kKMinus) // kaon - registry.fill(HIST("PtKaVsCentMC"), particle.pt(), centrality); - else if (particle.pdgCode() == PDG_t::kProton || particle.pdgCode() == PDG_t::kProtonBar) // proton - registry.fill(HIST("PtPrVsCentMC"), particle.pt(), centrality); - else - continue; - } - } - } // Collisions - } - PROCESS_SWITCH(PiKpRAA, processMC, "Process MC closure", false); void processSim(aod::McCollisions::iterator const& mccollision, soa::SmallGroups const& collisions, BCsRun3 const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::McParticles const& mcParticles, TracksMC const& tracksMC) { @@ -1412,6 +1144,7 @@ struct PiKpRAA { // By counting number of primary charged particles in |eta| < 1 //--------------------------- int nChMC{0}; + int nChMCEta08{0}; int nChFT0A{0}; int nChFT0C{0}; for (const auto& particle : mcParticles) { @@ -1444,6 +1177,10 @@ struct PiKpRAA { nChFT0C++; } + if (std::abs(eta) < tpcNchAcceptance) { + nChMCEta08++; + } + // INEL > 0 if (std::abs(eta) > kOne) continue; @@ -1557,9 +1294,24 @@ struct PiKpRAA { registry.fill(HIST("Centrality_WRecoEvt"), centrality); registry.fill(HIST("zPosMC"), mccollision.posZ()); + //--------------------------- + // Event selection + // for reconstructed collisions + //--------------------------- + if (!isEventSelected(collision)) { + continue; + } + + registry.fill(HIST("Centrality_WRecoEvtWSelCri"), centrality); + registry.fill(HIST("NchMCVsCent"), centrality, nChMCEta08); + registry.fill(HIST("NchMC_WithRecoEvt"), nChMCEta08); // Numerator of event loss correction + registry.fill(HIST("zPos"), collision.posZ()); + registry.fill(HIST("T0Ccent"), centrality); + //--------------------------- // All Generated events with at least one associated reconstructed collision // The Generated events are not subjected to any selection criteria + // However, the associated reconstructed collisions pass the selection criteria // This histograms are used for the denominator of the tracking efficiency //--------------------------- for (const auto& particle : mcParticles) { @@ -1590,34 +1342,19 @@ struct PiKpRAA { if (isPrimary) { if (particle.pdgCode() == PDG_t::kPiPlus || particle.pdgCode() == PDG_t::kPiMinus) { registry.fill(HIST("PtPiVsCentMC_WithRecoEvt"), particle.pt(), centrality); // Denominator of tracking efficiency - registry.fill(HIST("PtPiVsNchMC_WithRecoEvt"), particle.pt(), nChMC); // Numerator of signal loss + registry.fill(HIST("PtPiVsNchMC_WithRecoEvt"), particle.pt(), nChMCEta08); // Numerator of signal loss } else if (particle.pdgCode() == PDG_t::kKPlus || particle.pdgCode() == PDG_t::kKMinus) { registry.fill(HIST("PtKaVsCentMC_WithRecoEvt"), particle.pt(), centrality); - registry.fill(HIST("PtKaVsNchMC_WithRecoEvt"), particle.pt(), nChMC); + registry.fill(HIST("PtKaVsNchMC_WithRecoEvt"), particle.pt(), nChMCEta08); } else if (particle.pdgCode() == PDG_t::kProton || particle.pdgCode() == PDG_t::kProtonBar) { registry.fill(HIST("PtPrVsCentMC_WithRecoEvt"), particle.pt(), centrality); - registry.fill(HIST("PtPrVsNchMC_WithRecoEvt"), particle.pt(), nChMC); + registry.fill(HIST("PtPrVsNchMC_WithRecoEvt"), particle.pt(), nChMCEta08); } else { continue; } } } // Loop over generated particles per generated collision - //--------------------------- - // Reconstructed collisions subjected to selection criteria - //--------------------------- - - // Event selection - if (!isEventSelected(collision)) { - continue; - } - - registry.fill(HIST("Centrality_WRecoEvtWSelCri"), centrality); - registry.fill(HIST("NchMCVsCent"), centrality, nChMC); - registry.fill(HIST("NchMC_WithRecoEvt"), nChMC); - registry.fill(HIST("T0Ccent"), centrality); - registry.fill(HIST("zPos"), collision.posZ()); - const auto& groupedTracks{tracksMC.sliceBy(perCollision, collision.globalIndex())}; //--------------------------- @@ -1711,9 +1448,8 @@ struct PiKpRAA { } //--------------------------- + // Needed for the number of the Tracking Efficiency // Global track + DCAxy selections - // This is needed for the number of the Tracking Efficiency - // and the spectra to be corrected //--------------------------- int nCh{0}; for (const auto& track : groupedTracks) { @@ -1799,16 +1535,19 @@ struct PiKpRAA { } if (isPi && !isKa && !isPr) { - registry.fill(HIST("PtPiVsCent_WithRecoEvt"), track.pt(), centrality); // Numerator of tracking efficiency - registry.fill(HIST("MCclosure_PtPiVsNchMC"), track.pt(), nChMC); + registry.fill(HIST("PtPiVsCent_WithRecoEvt"), track.pt(), centrality); // Numerator of tracking efficiency + registry.fill(HIST("PtGenPiVsCent_WithRecoEvt"), particle.pt(), centrality); // Numerator of tracking efficiency + registry.fill(HIST("MCclosure_PtPiVsNchMC"), track.pt(), nChMCEta08); } if (isKa && !isPi && !isPr) { registry.fill(HIST("PtKaVsCent_WithRecoEvt"), track.pt(), centrality); - registry.fill(HIST("MCclosure_PtKaVsNchMC"), track.pt(), nChMC); + registry.fill(HIST("PtGenKaVsCent_WithRecoEvt"), particle.pt(), centrality); + registry.fill(HIST("MCclosure_PtKaVsNchMC"), track.pt(), nChMCEta08); } if (isPr && !isPi && !isKa) { registry.fill(HIST("PtPrVsCent_WithRecoEvt"), track.pt(), centrality); - registry.fill(HIST("MCclosure_PtPrVsNchMC"), track.pt(), nChMC); + registry.fill(HIST("PtGenPrVsCent_WithRecoEvt"), particle.pt(), centrality); + registry.fill(HIST("MCclosure_PtPrVsNchMC"), track.pt(), nChMCEta08); } registry.fill(HIST("PtResolution"), particle.pt(), (track.pt() - particle.pt()) / particle.pt()); } // Loop over reconstructed tracks @@ -1820,11 +1559,6 @@ struct PiKpRAA { // All Generated events irrespective of whether there is an associated reconstructed collision // Consequently, the centrality being a reconstructed quantity, might not always be available // Therefore it is expressed as a function of the generated pT and the generated Nch in ∣eta∣ < 0.8 - //--------------------------- - - //--------------------------- - // Generated Pt spectra of all INEL > 0 Generated evets - // irrespective of whether there is a reconstructed collision // This is used for the denominator of the signal loss correction // Also for MC closure: True Pt vs Generated Nch //--------------------------- @@ -1855,14 +1589,14 @@ struct PiKpRAA { if (isPrimary) { if (particle.pdgCode() == PDG_t::kPiPlus || particle.pdgCode() == PDG_t::kPiMinus) { - registry.fill(HIST("PtPiVsNchMC_AllGen"), particle.pt(), nChMC); - registry.fill(HIST("MCclosure_PtMCPiVsNchMC"), particle.pt(), nChMC); + registry.fill(HIST("PtPiVsNchMC_AllGen"), particle.pt(), nChMCEta08); + registry.fill(HIST("MCclosure_PtMCPiVsNchMC"), particle.pt(), nChMCEta08); } else if (particle.pdgCode() == PDG_t::kKPlus || particle.pdgCode() == PDG_t::kKMinus) { - registry.fill(HIST("PtKaVsNchMC_AllGen"), particle.pt(), nChMC); - registry.fill(HIST("MCclosure_PtMCKaVsNchMC"), particle.pt(), nChMC); + registry.fill(HIST("PtKaVsNchMC_AllGen"), particle.pt(), nChMCEta08); + registry.fill(HIST("MCclosure_PtMCKaVsNchMC"), particle.pt(), nChMCEta08); } else if (particle.pdgCode() == PDG_t::kProton || particle.pdgCode() == PDG_t::kProtonBar) { - registry.fill(HIST("PtPrVsNchMC_AllGen"), particle.pt(), nChMC); - registry.fill(HIST("MCclosure_PtMCPrVsNchMC"), particle.pt(), nChMC); + registry.fill(HIST("PtPrVsNchMC_AllGen"), particle.pt(), nChMCEta08); + registry.fill(HIST("MCclosure_PtMCPrVsNchMC"), particle.pt(), nChMCEta08); } else { continue; } @@ -1872,7 +1606,7 @@ struct PiKpRAA { //--------------------------- // This is used for the denominator of the event loss correction //--------------------------- - registry.fill(HIST("NchMC_AllGen"), nChMC); + registry.fill(HIST("NchMC_AllGen"), nChMCEta08); } PROCESS_SWITCH(PiKpRAA, processSim, "Process Sim", false);