From e3854c6ed5040c7cd3a6ea574dd18a27f9bedde1 Mon Sep 17 00:00:00 2001 From: Jaideep Tanwar <141036812+jtanwar2212@users.noreply.github.com> Date: Fri, 5 Sep 2025 17:25:45 +0530 Subject: [PATCH 1/8] Update nucleitpcpbpb.cxx --- PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx | 519 ++++++++++++++++++++++----- 1 file changed, 427 insertions(+), 92 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx index fc99ebcd30b..afcc0a596ee 100644 --- a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx +++ b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx @@ -59,6 +59,12 @@ static const std::vector particlePdgCodes{211, 2212, o2::constants::physics static const std::vector particleMasses{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton, o2::constants::physics::MassDeuteron, o2::constants::physics::MassTriton, o2::constants::physics::MassHelium3, o2::constants::physics::MassAlpha}; static const std::vector particleCharge{1, 1, 1, 1, 2, 2}; const int nBetheParams = 6; +std::vector hfMothCodes = {511, 521, 531, 541, 5122}; // b-mesons + Lambda_b +enum NucleiFlags { + kIsPhysicalPrimary = BIT(0), + kIsSecondaryFromMaterial = BIT(1), + kIsSecondaryFromWeakDecay = BIT(2) +}; static const std::vector betheBlochParNames{"p0", "p1", "p2", "p3", "p4", "resolution"}; constexpr double kBetheBlochDefault[nParticles][nBetheParams]{ {13.611469, 3.598765, -0.021138, 2.039562, 0.651040, 0.09}, // pion @@ -103,7 +109,7 @@ struct PrimParticles { } }; // struct PrimParticles //---------------------------------------------------------------------------------------------------------------- -std::vector> hmass; +std::vector> hmass; std::vector> hmassnsigma; } // namespace //---------------------------------------------------------------------------------------------------------------- @@ -144,6 +150,7 @@ struct NucleitpcPbPb { Configurable cfgmaxGetMeanItsClsSizeRequire{"cfgmaxGetMeanItsClsSizeRequire", true, "Require maxGetMeanItsClsSize Cut"}; Configurable cfgDCAwithptRequire{"cfgDCAwithptRequire", true, "Require DCA cuts with pt dependance"}; Configurable cfgRequirebetaplot{"cfgRequirebetaplot", true, "Require beta plot"}; + Configurable cfgIncludeMaterialInEfficiency{"cfgIncludeMaterialInEfficiency", false, "Require from material in efficiency"}; Configurable> cfgBetheBlochParams{"cfgBetheBlochParams", {kBetheBlochDefault[0], nParticles, nBetheParams, particleNames, betheBlochParNames}, "TPC Bethe-Bloch parameterisation for light nuclei"}; Configurable> cfgTrackPIDsettings{"cfgTrackPIDsettings", {kTrackPIDSettings[0], nParticles, nTrkSettings, particleNames, trackPIDsettingsNames}, "track selection and PID criteria"}; @@ -165,14 +172,18 @@ struct NucleitpcPbPb { ConfigurableAxis axisdEdx{"axisdEdx", {4000, 0, 4000}, "d#it{E}/d#it{x}"}; ConfigurableAxis axisCent{"axisCent", {100, 0, 100}, "centrality"}; ConfigurableAxis axisVtxZ{"axisVtxZ", {120, -20, 20}, "z"}; - ConfigurableAxis axisImpt{"axisImpt", {100, 0, 20}, "impact parameter"}; + ConfigurableAxis ptAxis{"ptAxis", {200, 0, 10}, "#it{p}_{T} (GeV/#it{c})"}; ConfigurableAxis axiseta{"axiseta", {100, -1, 1}, "eta"}; ConfigurableAxis axisrapidity{"axisrapidity", {100, -2, 2}, "rapidity"}; - ConfigurableAxis axismass{"axismass", {100, -10, 10}, "mass^{2}"}; + ConfigurableAxis axismass{"axismass", {100, -10, 10}, "mass"}; + ConfigurableAxis axismassnsigma{"axismassnsigma", {100, 0, 20}, "mass"}; ConfigurableAxis nsigmaAxis{"nsigmaAxis", {160, -10, 10}, "n#sigma_{#pi^{+}}"}; ConfigurableAxis speciesBitAxis{"speciesBitAxis", {8, -0.5, 7.5}, "particle type 0: pion, 1: proton, 2: deuteron, 3: triton, 4:He3, 5:He4"}; + ConfigurableAxis speciesTrackingAxis{"speciesTrackingAxis", {11, -0.5, 10.5}, "particle type 0: pion, 1: proton, 2: deuteron, 3: triton, 4:He3, 5:He4"}; ConfigurableAxis axisDCA{"axisDCA", {400, -10., 10.}, "DCA axis"}; + ConfigurableAxis particleAntiAxis{"particleAntiAxis", {2, 0, 2}, "Particle/Anti-particle"}; // 0 = particle, 1 = anti-particle + ConfigurableAxis decayTypeAxis{"decayTypeAxis", {3, 0, 3}, "Decay type"}; // 0 = primary, 1 = from decay // CCDB Service ccdb; @@ -213,6 +224,7 @@ struct NucleitpcPbPb { histos.add("histCentFT0C", "histCentFT0C", kTH1F, {axisCent}); histos.add("histCentFT0M", "histCentFT0M", kTH1F, {axisCent}); histos.add("histCentFTOC_cut", "histCentFTOC_cut", kTH1F, {axisCent}); + histos.add("hSpectra", " ", HistType::kTHnSparseF, {speciesBitAxis, ptAxis, nsigmaAxis, {5, -2.5, 2.5}, axisCent, axisDCA, axisDCA}); } histos.add("histeta", "histeta", kTH1F, {axiseta}); histos.add("Tofsignal", "Tofsignal", kTH2F, {axisRigidity, {4000, 0.2, 1.2, "#beta"}}); @@ -224,43 +236,71 @@ struct NucleitpcPbPb { for (int i = 0; i < nParticles; i++) { TString histName = primaryParticles[i].name; if (cfgFillmass) { - hmass[2 * i] = histos.add(Form("histmass_pt/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}", HistType::kTH2F, {ptAxis, axismass}); - hmass[2 * i + 1] = histos.add(Form("histmass_ptanti/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}", HistType::kTH2F, {ptAxis, axismass}); + hmass[2 * i] = histos.add(Form("histmass_pt/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}; centrality(%)", HistType::kTH3F, {ptAxis, axismass, axisCent}); + hmass[2 * i + 1] = histos.add(Form("histmass_ptanti/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}; centrality(%)", HistType::kTH3F, {ptAxis, axismass, axisCent}); } } for (int i = 0; i < nParticles; i++) { TString histName = primaryParticles[i].name; if (cfgFillmassnsigma) { - hmassnsigma[2 * i] = histos.add(Form("histmass_nsigma/histmass_%s", histName.Data()), ";nsigma; mass^{2}", HistType::kTH2F, {nsigmaAxis, axismass}); - hmassnsigma[2 * i + 1] = histos.add(Form("histmass_nsigmaanti/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}", HistType::kTH2F, {nsigmaAxis, axismass}); + hmassnsigma[2 * i] = histos.add(Form("histmass_nsigma/histmass_%s", histName.Data()), ";nsigma; mass^{2}", HistType::kTH2F, {nsigmaAxis, axismassnsigma}); + hmassnsigma[2 * i + 1] = histos.add(Form("histmass_nsigmaanti/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}", HistType::kTH2F, {nsigmaAxis, axismassnsigma}); } } - histos.add("hSpectra", " ", HistType::kTHnSparseF, {speciesBitAxis, ptAxis, nsigmaAxis, {5, -2.5, 2.5}, axisCent, axisDCA, axisDCA, axisrapidity, axiseta}); - if (doprocessMC) { + histomc.add("hSpectramc", " ", HistType::kTHnSparseF, {speciesBitAxis, {5, -2.5, 2.5}, axisCent, ptAxis, ptAxis}); + + histomc.add("hDenomEffAcc", "Denominator for Efficiency x Acceptance", + {HistType::kTHnSparseF, {speciesBitAxis, ptAxis, axisrapidity, axisCent, particleAntiAxis, decayTypeAxis}}); + histomc.add("hNumerEffAcc", "Numerator for Efficiency x Acceptance", + {HistType::kTHnSparseF, {speciesBitAxis, ptAxis, axisrapidity, axisCent, particleAntiAxis, decayTypeAxis}}); + histomc.add("hDenomSignalLoss", "Denominator for Signal Loss", + {HistType::kTHnSparseF, {speciesBitAxis, ptAxis, axisrapidity, axisCent, particleAntiAxis, decayTypeAxis}}); + histomc.add("hNumerSignalLoss", "Numerator for Signal Loss", + {HistType::kTHnSparseF, {speciesBitAxis, ptAxis, axisrapidity, axisCent, particleAntiAxis, decayTypeAxis}}); + + // Material secondary histograms + histomc.add("histSecondaryMaterialHe3", "He3 from material", HistType::kTH2F, {ptAxis, axisCent}); + histomc.add("histSecondaryMaterialAntiHe3", "Anti-He3 from material", HistType::kTH2F, {ptAxis, axisCent}); + histomc.add("histSecondaryMaterialHe4", "He4 from material", HistType::kTH2F, {ptAxis, axisCent}); + histomc.add("histSecondaryMaterialAntiHe4", "Anti-He4 from material", HistType::kTH2F, {ptAxis, axisCent}); + histomc.add("histVtxZgen", "histVtxZgen", kTH1F, {axisVtxZ}); - histomc.add("ImptParameter", "ImptParameter", kTH1F, {axisImpt}); histomc.add("histEtagen", "histEtagen", kTH1F, {axiseta}); histomc.add("histPtgenHe3", "histPtgenHe3", kTH1F, {ptAxis}); histomc.add("histPtgenAntiHe3", "histPtgenAntiHe3", kTH1F, {ptAxis}); histomc.add("histPtgenHe4", "histPtgenHe4", kTH1F, {ptAxis}); histomc.add("histPtgenAntiHe4", "histPtgenAntiHe4", kTH1F, {ptAxis}); - // Reconstrcuted eta histomc.add("histNevReco", "histNevReco", kTH1F, {axisNev}); histomc.add("histVtxZReco", "histVtxZReco", kTH1F, {axisVtxZ}); histomc.add("histCentFT0CReco", "histCentFT0CReco", kTH1F, {axisCent}); histomc.add("histCentFT0MReco", "histCentFT0MReco", kTH1F, {axisCent}); - histomc.add("histPtRecoHe3", "histPtgenHe3", kTH1F, {ptAxis}); - histomc.add("histPtRecoAntiHe3", "histPtgenAntiHe3", kTH1F, {ptAxis}); - histomc.add("histPtRecoHe4", "histPtgenHe4", kTH1F, {ptAxis}); - histomc.add("histPtRecoAntiHe4", "histPtgenAntiHe4", kTH1F, {ptAxis}); histomc.add("histDeltaPtVsPtGen", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10}, {1000, -0.5, 0.5, "p_{T}(reco) - p_{T}(gen);p_{T}(reco)"}}); histomc.add("histDeltaPtVsPtGenanti", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10}, {1000, -0.5, 0.5, "p_{T}(reco) - p_{T}(gen);p_{T}(reco)"}}); - histomc.add("histPIDtrack", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); - histomc.add("histPIDtrackanti", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); histomc.add("histDeltaPtVsPtGenHe4", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10}, {1000, -0.5, 0.5, "p_{T}(reco) - p_{T}(gen);p_{T}(reco)"}}); histomc.add("histDeltaPtVsPtGenHe4anti", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10}, {1000, -0.5, 0.5, "p_{T}(reco) - p_{T}(gen);p_{T}(reco)"}}); + histomc.add("histPIDtrack", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); + histomc.add("histPIDtrackanti", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); + histomc.add("histPIDtrackhe4", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); + histomc.add("histPIDtrackantihe4", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); + + histomc.add("histWeakDecayPtHe3", "Pt distribution of He3 from weak decays", kTH2F, {ptAxis, axisCent}); + histomc.add("histWeakDecayPtAntiHe3", "Pt distribution of Anti-He3 from weak decays", kTH2F, {ptAxis, axisCent}); + histomc.add("histWeakDecayPtHe4", "Pt distribution of He4 from weak decays", kTH2F, {ptAxis, axisCent}); + histomc.add("histWeakDecayPtAntiHe4", "Pt distribution of Anti-He4 from weak decays", kTH2F, {ptAxis, axisCent}); + histomc.add("histRecoWeakDecayPtHe3", "Pt distribution of reconstructed He3 from weak decays", kTH2F, {ptAxis, axisCent}); + histomc.add("histRecoWeakDecayPtHe4", "Pt distribution of reconstructed He4 from weak decays", kTH2F, {ptAxis, axisCent}); + histomc.add("histRecoSecondaryMaterialHe3", "Reco He3 from material", HistType::kTH2F, {ptAxis, axisCent}); + histomc.add("histRecoSecondaryMaterialAntiHe3", "Reco Anti-He3 from material", HistType::kTH2F, {ptAxis, axisCent}); + histomc.add("histRecoSecondaryMaterialHe4", "Reco He4 from material", HistType::kTH2F, {ptAxis, axisCent}); + histomc.add("histRecoSecondaryMaterialAntiHe4", "Reco Anti-He4 from material", HistType::kTH2F, {ptAxis, axisCent}); + histomc.add("histProcessCodeHe", "Process codes for He isotopes", kTH1F, {{100, 0, 100, "process code"}}); + histomc.add("histProcess23Details", "Process 23 details", kTH2F, {{4, 0.5, 4.5, "particle type"}, {100, 0, 10, "p_{T}"}}); + histomc.add("histAllMaterialSecondariesGen", "All material secondaries (gen)", kTH3F, {{100, 0, 10, "p_{T}"}, {20, -1, 1, "y"}, {5, -0.5, 4.5, "type"}}); + histomc.add("histAllMaterialSecondariesReco", "All material secondaries (reco)", kTH3F, {{100, 0, 10, "p_{T}"}, {20, -1, 1, "y"}, {5, -0.5, 4.5, "type"}}); + + // Add axis labels for type: 0=unknown, 1=He3, 2=anti-He3, 3=He4, 4=anti-He4 } } //---------------------------------------------------------------------------------------------------------------- @@ -286,6 +326,7 @@ struct NucleitpcPbPb { if (removeNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) continue; histos.fill(HIST("histCentFTOC_cut"), collision.centFT0C()); + histos.fill(HIST("histNev"), 2.5); const uint64_t collIdx = collision.globalIndex(); auto tracksByColl = tracksColl.sliceBy(perCollision, collIdx); ///////////////////////////////////////////////////////////////////////////////// @@ -353,12 +394,12 @@ struct NucleitpcPbPb { histos.fill(HIST("Tpcsignal"), getRigidity(track) * track.sign(), track.tpcSignal()); if (cfgFillhspectra && cfgTrackPIDsettings2->get(i, "fillsparsh") == 1) { - histos.fill(HIST("hSpectra"), i, ptMomn, tpcNsigma, sign, collision.centFT0C(), track.dcaZ(), track.dcaXY(), getRapidity(track, i), track.eta()); + histos.fill(HIST("hSpectra"), i, ptMomn, tpcNsigma, sign, collision.centFT0C(), track.dcaZ(), track.dcaXY()); } fillhmassnsigma(track, i, tpcNsigma); if ((std::abs(tpcNsigma) > cfgTrackPIDsettings2->get(i, "maxTPCnsigmaTOF")) && cfgTrackPIDsettings2->get(i, "useTPCnsigmaTOF") < 1) continue; - fillhmass(track, i); + fillhmass(track, i, collision.centFT0C()); if (cfgRequirebetaplot) { histos.fill(HIST("Tofsignal"), getRigidity(track) * track.sign(), o2::pid::tof::Beta::GetBeta(track)); @@ -370,58 +411,275 @@ struct NucleitpcPbPb { } } PROCESS_SWITCH(NucleitpcPbPb, processData, "data analysis", false); + //---------------------------------------------------------------------------------------------------------------- - // MC particles - //-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= + // MC particles - Efficiency x Acceptance and Signal Loss calculations + //-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- struct McCollInfo { bool passedEvSel = false; + float centrality = -1.0f; }; std::vector mcCollInfos; - void processMC(CollisionsFullMC const& collisions, aod::McCollisions const& mcCollisions, TracksFull const& tracks, aod::BCsWithTimestamps const& bcs, aod::McParticles const& particlesMC, aod::McTrackLabels const& trackLabelsMC, aod::McCollisionLabels const& collLabels, aod::TrackAssoc const& tracksColl, CollisionsFull const& colls) + + void processMC(CollisionsFullMC const& collisions, + aod::McCollisions const& mcCollisions, + TracksFull const& tracks, + aod::BCsWithTimestamps const& bcs, + aod::McParticles const& particlesMC, + aod::McTrackLabels const& trackLabelsMC, + aod::McCollisionLabels const& collLabels, + aod::TrackAssoc const& tracksColl, + CollisionsFull const& colls) { - (void)colls; - (void)collLabels; + (void)bcs; + (void)collLabels; + (void)colls; + mcCollInfos.clear(); mcCollInfos.resize(mcCollisions.size()); - // ----------------------------- Generated particles loop ----------------------------- - for (auto const& mcColl : mcCollisions) { - if (std::abs(mcColl.posZ()) > cfgZvertex && cfgZvertexRequireMC) + for (auto const& collision : collisions) { + int mcCollIdx = collision.mcCollisionId(); + if (mcCollIdx < 0 || mcCollIdx >= static_cast(mcCollisions.size())) { + continue; + } + if (std::abs(collision.posZ()) > cfgZvertex) + continue; + if (!collision.sel8()) + continue; + if (collision.centFT0C() > centcut) + continue; + + mcCollInfos[mcCollIdx].passedEvSel = true; + mcCollInfos[mcCollIdx].centrality = collision.centFT0C(); + } + for (auto const& mcCollision : mcCollisions) { + size_t idx = mcCollision.globalIndex(); + + if (idx >= mcCollInfos.size()) continue; - histomc.fill(HIST("histVtxZgen"), mcColl.posZ()); - histomc.fill(HIST("ImptParameter"), mcColl.impactParameter()); + for (auto const& mcParticle : particlesMC) { - if (mcParticle.mcCollisionId() != mcColl.globalIndex()) + + if (mcParticle.mcCollisionId() != mcCollision.globalIndex()) continue; - if (!mcParticle.isPhysicalPrimary()) + + int pdgCode = mcParticle.pdgCode(); + bool isHe3 = (std::abs(pdgCode) == 1000020030); + bool isHe4 = (std::abs(pdgCode) == 1000020040); + + if (std::abs(mcParticle.eta()) > cfgCutEta) + continue; + if (std::abs(mcParticle.y()) > cfgCutRapidity) + continue; + bool isMaterialSecondary = false; + if (!mcParticle.isPhysicalPrimary() && (isHe3 || isHe4)) { + if (!mcParticle.has_mothers()) { + isMaterialSecondary = true; + } else { + auto mother = mcParticle.mothers_as().front(); + if (std::abs(mother.pdgCode()) < 1000000000) { + isMaterialSecondary = true; + } + } + } + if (isHe3 || isHe4) { + histomc.fill(HIST("histProcessCodeHe"), mcParticle.getProcess()); + if (mcParticle.getProcess() == 23) { + histomc.fill(HIST("histProcess23Details"), std::abs(pdgCode), mcParticle.pt()); + } + } + + if (isMaterialSecondary) { + float centrality = mcCollInfos[idx].passedEvSel ? mcCollInfos[idx].centrality : -1.0f; + + if (pdgCode == 1000020030) { + histomc.fill(HIST("histSecondaryMaterialHe3"), mcParticle.pt(), centrality); + } else if (pdgCode == -1000020030) { + histomc.fill(HIST("histSecondaryMaterialAntiHe3"), mcParticle.pt(), centrality); + } else if (pdgCode == 1000020040) { + histomc.fill(HIST("histSecondaryMaterialHe4"), mcParticle.pt(), centrality); + } else if (pdgCode == -1000020040) { + histomc.fill(HIST("histSecondaryMaterialAntiHe4"), mcParticle.pt(), centrality); + } + int type = 0; + if (std::abs(pdgCode) == 1000020030) + type = (pdgCode > 0) ? 1 : 2; + else if (std::abs(pdgCode) == 1000020040) + type = (pdgCode > 0) ? 3 : 4; + histomc.fill(HIST("histAllMaterialSecondariesGen"), mcParticle.pt(), mcParticle.y(), type); + } + + if (!isHe3 && !isHe4) + continue; + + uint16_t flags = 0; + int motherPdgCode = 0; + float motherDecRadius = -1; + int decayType = 0; + int particleAnti = (pdgCode > 0) ? 0 : 1; + + if (mcParticle.isPhysicalPrimary()) { + flags |= kIsPhysicalPrimary; + decayType = 0; + if (mcParticle.has_mothers()) { + for (auto& motherparticle : mcParticle.mothers_as()) { + if (std::find(hfMothCodes.begin(), hfMothCodes.end(), + std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { + flags |= kIsSecondaryFromWeakDecay; + decayType = 1; + motherPdgCode = motherparticle.pdgCode(); + motherDecRadius = std::hypot(mcParticle.vx() - motherparticle.vx(), + mcParticle.vy() - motherparticle.vy()); + break; + } + } + } + } else if (mcParticle.has_mothers()) { + flags |= kIsSecondaryFromWeakDecay; + decayType = 1; + for (auto& motherparticle : mcParticle.mothers_as()) { + motherPdgCode = motherparticle.pdgCode(); + motherDecRadius = std::hypot(mcParticle.vx() - motherparticle.vx(), + mcParticle.vy() - motherparticle.vy()); + } + } else { + flags |= kIsSecondaryFromMaterial; + decayType = 2; + continue; + } + bool isFromWeakDecay = (flags & kIsSecondaryFromWeakDecay) != 0; + + if (!mcParticle.isPhysicalPrimary() && !isFromWeakDecay) + continue; + + int particleType = -1; + if (std::abs(pdgCode) == 1000020030) + particleType = he3; + else if (std::abs(pdgCode) == 1000020040) + particleType = he4; + + if (particleType >= 0) { + float centrality = mcCollInfos[idx].passedEvSel ? mcCollInfos[idx].centrality : -1.0f; + histomc.fill(HIST("hDenomSignalLoss"), particleType, mcParticle.pt(), mcParticle.y(), centrality, particleAnti, decayType); + + if (mcCollInfos[idx].passedEvSel) { + histomc.fill(HIST("hNumerSignalLoss"), particleType, mcParticle.pt(), mcParticle.y(), centrality, particleAnti, decayType); + } + } + + if (isFromWeakDecay) { + float centrality = mcCollInfos[idx].passedEvSel ? mcCollInfos[idx].centrality : -1.0f; + if (pdgCode == 1000020030) { + histomc.fill(HIST("histWeakDecayPtHe3"), mcParticle.pt(), centrality); + } else if (pdgCode == -1000020030) { + histomc.fill(HIST("histWeakDecayPtAntiHe3"), mcParticle.pt(), centrality); + } else if (pdgCode == 1000020040) { + histomc.fill(HIST("histWeakDecayPtHe4"), mcParticle.pt(), centrality); + } else if (pdgCode == -1000020040) { + histomc.fill(HIST("histWeakDecayPtAntiHe4"), mcParticle.pt(), centrality); + } + } + if (pdgCode == 1000020030) { + histomc.fill(HIST("histPtgenHe3"), mcParticle.pt()); + } else if (pdgCode == -1000020030) { + histomc.fill(HIST("histPtgenAntiHe3"), mcParticle.pt()); + } else if (pdgCode == 1000020040) { + histomc.fill(HIST("histPtgenHe4"), mcParticle.pt()); + } else if (pdgCode == -1000020040) { + histomc.fill(HIST("histPtgenAntiHe4"), mcParticle.pt()); + } + } + } + for (auto const& mcCollision : mcCollisions) { + size_t idx = mcCollision.globalIndex(); + + if (idx >= mcCollInfos.size()) + continue; + if (!mcCollInfos[idx].passedEvSel) + continue; + if (std::abs(mcCollision.posZ()) > cfgZvertex) + continue; + + histomc.fill(HIST("histVtxZgen"), mcCollision.posZ()); + float centrality = mcCollInfos[idx].centrality; + + for (auto const& mcParticle : particlesMC) { + if (mcParticle.mcCollisionId() != mcCollision.globalIndex()) continue; + int pdgCode = mcParticle.pdgCode(); - if (std::abs(mcParticle.y()) > cfgCutRapidity && cfgrapidityRequireMC) + bool isHe3 = (std::abs(pdgCode) == 1000020030); + bool isHe4 = (std::abs(pdgCode) == 1000020040); + + if (!isHe3 && !isHe4) + continue; + if (std::abs(mcParticle.eta()) > cfgCutEta) + continue; + if (std::abs(mcParticle.y()) > cfgCutRapidity) continue; - if (std::abs(mcParticle.eta()) > cfgCutEta && cfgetaRequireMC) + + uint16_t flags = 0; + int motherPdgCode = 0; + float motherDecRadius = -1; + int decayType = 0; + int particleAnti = (pdgCode > 0) ? 0 : 1; + + if (mcParticle.isPhysicalPrimary()) { + flags |= kIsPhysicalPrimary; + decayType = 0; + if (mcParticle.has_mothers()) { + for (auto& motherparticle : mcParticle.mothers_as()) { + if (std::find(hfMothCodes.begin(), hfMothCodes.end(), + std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { + flags |= kIsSecondaryFromWeakDecay; + decayType = 1; + motherPdgCode = motherparticle.pdgCode(); + motherDecRadius = std::hypot(mcParticle.vx() - motherparticle.vx(), + mcParticle.vy() - motherparticle.vy()); + break; + } + } + } + } else if (mcParticle.has_mothers()) { + flags |= kIsSecondaryFromWeakDecay; + decayType = 1; + for (auto& motherparticle : mcParticle.mothers_as()) { + motherPdgCode = motherparticle.pdgCode(); + motherDecRadius = std::hypot(mcParticle.vx() - motherparticle.vx(), + mcParticle.vy() - motherparticle.vy()); + } + } else { + flags |= kIsSecondaryFromMaterial; + decayType = 2; continue; - histomc.fill(HIST("histEtagen"), mcParticle.eta()); - float ptScaled = mcParticle.pt(); - if (pdgCode == particlePdgCodes.at(4)) { - histomc.fill(HIST("histPtgenHe3"), ptScaled); - } else if (pdgCode == -particlePdgCodes.at(4)) { - histomc.fill(HIST("histPtgenAntiHe3"), ptScaled); - } else if (pdgCode == particlePdgCodes.at(5)) { - histomc.fill(HIST("histPtgenHe4"), ptScaled); - } else if (pdgCode == -particlePdgCodes.at(5)) { - histomc.fill(HIST("histPtgenAntiHe4"), ptScaled); } - } // mc track loop generated - } // mc collision loop generated - // ----------------------------- Reconstructed track loop ----------------------------- + bool isFromWeakDecay = (flags & kIsSecondaryFromWeakDecay) != 0; + + if (!mcParticle.isPhysicalPrimary() && !isFromWeakDecay) + continue; + + int particleType = -1; + if (std::abs(pdgCode) == 1000020030) + particleType = he3; + else if (std::abs(pdgCode) == 1000020040) + particleType = he4; + + if (particleType >= 0) { + histomc.fill(HIST("hDenomEffAcc"), particleType, mcParticle.pt(), mcParticle.y(), centrality, particleAnti, decayType); + } + } + } for (auto const& collision : collisions) { auto mcCollIdx = collision.mcCollisionId(); - if (mcCollIdx < 0 || mcCollIdx >= mcCollisions.size()) + if (mcCollIdx < 0 || mcCollIdx >= static_cast(mcCollisions.size())) continue; + auto bc = collision.bc_as(); initCCDB(bc); collHasCandidate = false; histomc.fill(HIST("histNevReco"), 0.5); + if (std::abs(collision.posZ()) > cfgZvertex && cfgZvertexRequire) continue; collPassedEvSel = collision.sel8(); @@ -444,17 +702,89 @@ struct NucleitpcPbPb { continue; if (removeNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) continue; - // Loop through TrackAssoc and check if this track belongs to the current collision for (const auto& assoc : tracksColl) { if (assoc.collisionId() != collision.globalIndex()) continue; const auto& track = tracks.rawIteratorAt(assoc.trackId()); auto labelRow = trackLabelsMC.iteratorAt(track.globalIndex()); int label = labelRow.mcParticleId(); - if (label < 0 || label >= particlesMC.size()) + if (label < 0 || label >= static_cast(particlesMC.size())) continue; auto const& matchedMCParticle = particlesMC.iteratorAt(label); + int pdg = matchedMCParticle.pdgCode(); + bool isHe3 = (std::abs(pdg) == 1000020030); + bool isHe4 = (std::abs(pdg) == 1000020040); + + bool isMaterialSecondary = false; + if ((isHe3 || isHe4) && !matchedMCParticle.isPhysicalPrimary()) { + + if (matchedMCParticle.getProcess() == 23) { // kPMaterial + isMaterialSecondary = true; + } else if (!matchedMCParticle.has_mothers()) { + isMaterialSecondary = true; + } + } + + if (isMaterialSecondary) { + if (pdg == 1000020030) { + histomc.fill(HIST("histRecoSecondaryMaterialHe3"), track.pt(), collision.centFT0C()); + } else if (pdg == -1000020030) { + histomc.fill(HIST("histRecoSecondaryMaterialAntiHe3"), track.pt(), collision.centFT0C()); + } else if (pdg == 1000020040) { + histomc.fill(HIST("histRecoSecondaryMaterialHe4"), track.pt(), collision.centFT0C()); + } else if (pdg == -1000020040) { + histomc.fill(HIST("histRecoSecondaryMaterialAntiHe4"), track.pt(), collision.centFT0C()); + } + + int type = 0; + if (std::abs(pdg) == 1000020030) + type = (pdg > 0) ? 1 : 2; + else if (std::abs(pdg) == 1000020040) + type = (pdg > 0) ? 3 : 4; + histomc.fill(HIST("histAllMaterialSecondariesReco"), track.pt(), getRapidity(track, he3), type); + + if (!cfgIncludeMaterialInEfficiency) { + continue; + } + } + uint16_t flags = 0; + int motherPdgCode = 0; + float motherDecRadius = -1; + int decayType = 0; + bool isFromWeakDecay = false; + + if (matchedMCParticle.isPhysicalPrimary()) { + flags |= kIsPhysicalPrimary; + decayType = 0; + if (matchedMCParticle.has_mothers()) { + for (auto& motherparticle : matchedMCParticle.mothers_as()) { + if (std::find(hfMothCodes.begin(), hfMothCodes.end(), + std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { + flags |= kIsSecondaryFromWeakDecay; + isFromWeakDecay = true; + decayType = 1; + motherPdgCode = motherparticle.pdgCode(); + motherDecRadius = std::hypot(matchedMCParticle.vx() - motherparticle.vx(), + matchedMCParticle.vy() - motherparticle.vy()); + break; + } + } + } + } else if (matchedMCParticle.has_mothers()) { + flags |= kIsSecondaryFromWeakDecay; + isFromWeakDecay = true; + decayType = 1; + for (auto& motherparticle : matchedMCParticle.mothers_as()) { + motherPdgCode = motherparticle.pdgCode(); + motherDecRadius = std::hypot(matchedMCParticle.vx() - motherparticle.vx(), + matchedMCParticle.vy() - motherparticle.vy()); + } + } else { + flags |= kIsSecondaryFromMaterial; + decayType = 2; + } + if (!track.isPVContributor() && cfgUsePVcontributors) continue; if (!track.hasITS() && cfgITSrequire) @@ -467,12 +797,12 @@ struct NucleitpcPbPb { continue; if (std::abs(track.eta()) > cfgCutEta && cfgetaRequire) continue; - if (!matchedMCParticle.isPhysicalPrimary()) + if (!matchedMCParticle.isPhysicalPrimary() && !isFromWeakDecay && !isMaterialSecondary) continue; + for (size_t i = 0; i < primaryParticles.size(); i++) { if (std::abs(pdg) != std::abs(particlePdgCodes.at(i))) continue; - float ptMomn; setTrackParCov(track, mTrackParCov); mTrackParCov.setPID(track.pidForTracking()); @@ -488,7 +818,7 @@ struct NucleitpcPbPb { continue; if (track.tpcNClsFound() < cfgTrackPIDsettings->get(i, "minTPCnCls") && cfgTPCNClsfoundRequire) continue; - if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable()) && cfgTPCNClsCrossedRowsRequire) // o2-linter: disable=magic-number (To be checked) + if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable()) && cfgTPCNClsCrossedRowsRequire) continue; if (track.tpcChi2NCl() > cfgTrackPIDsettings->get(i, "maxTPCchi2") && cfgmaxTPCchi2Require) continue; @@ -506,7 +836,7 @@ struct NucleitpcPbPb { if (getMeanItsClsSize(track) < cfgTrackPIDsettings->get(i, "minITSclsSize") && cfgminGetMeanItsClsSizeRequire) continue; - bool insideDCAxy = (std::abs(track.dcaXY()) <= (cfgTrackPIDsettings->get(i, "maxDcaXY") * (0.0105f + 0.0350f / std::pow(ptMomn, 1.1f)))); // o2-linter: disable=magic-number (To be checked) + bool insideDCAxy = (std::abs(track.dcaXY()) <= (cfgTrackPIDsettings->get(i, "maxDcaXY") * (0.0105f + 0.0350f / std::pow(ptMomn, 1.1f)))); if ((!(insideDCAxy) || std::abs(track.dcaZ()) > dcazSigma(ptMomn, cfgTrackPIDsettings->get(i, "maxDcaZ"))) && cfgDCAwithptRequire) continue; @@ -518,24 +848,25 @@ struct NucleitpcPbPb { continue; if (itsSigma > cfgTrackPIDsettings2->get(i, "maxITSnsigma") && cfgTrackPIDsettings2->get(i, "useITSnsigma") < 1) continue; - histos.fill(HIST("Tpcsignal"), getRigidity(track) * track.sign(), track.tpcSignal()); + float ptReco; + mTrackParCov.setPID(track.pidForTracking()); + ptReco = (std::abs(pdg) == particlePdgCodes.at(4) || std::abs(pdg) == particlePdgCodes.at(5)) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); - if (cfgFillhspectra && cfgTrackPIDsettings2->get(i, "fillsparsh") == 1) { - histos.fill(HIST("hSpectra"), i, ptMomn, tpcNsigma, sign, collision.centFT0C(), track.dcaZ(), track.dcaXY(), getRapidity(track, i), track.eta()); - } - fillhmassnsigma(track, i, tpcNsigma); - if ((std::abs(tpcNsigma) > cfgTrackPIDsettings2->get(i, "maxTPCnsigmaTOF")) && cfgTrackPIDsettings2->get(i, "useTPCnsigmaTOF") < 1) - continue; - fillhmass(track, i); + int particleAnti = (pdg > 0) ? 0 : 1; - if (cfgRequirebetaplot) { - histos.fill(HIST("Tofsignal"), getRigidity(track) * track.sign(), o2::pid::tof::Beta::GetBeta(track)); + if (i == he3 || i == he4) { + histomc.fill(HIST("hNumerEffAcc"), i, matchedMCParticle.pt(), matchedMCParticle.y(), collision.centFT0C(), particleAnti, decayType); + } + if (isFromWeakDecay) { + if (std::abs(pdg) == 1000020030) { + histomc.fill(HIST("histRecoWeakDecayPtHe3"), ptReco, collision.centFT0C()); + } else if (std::abs(pdg) == 1000020040) { + histomc.fill(HIST("histRecoWeakDecayPtHe4"), ptReco, collision.centFT0C()); + } } - /*----------------------------------------------------------------------------------------------------------------*/ - float ptReco; - mTrackParCov.setPID(track.pidForTracking()); - ptReco = (std::abs(pdg) == particlePdgCodes.at(4) || std::abs(pdg) == particlePdgCodes.at(5)) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); + histos.fill(HIST("Tpcsignal"), getRigidity(track) * track.sign(), track.tpcSignal()); + if (pdg == -particlePdgCodes.at(5) && cfgmccorrectionhe4Require) { ptReco = ptReco + 0.00765 + 0.503791 * std::exp(-1.10517 * ptReco); } @@ -545,10 +876,10 @@ struct NucleitpcPbPb { int antitriton = 6; if (pidGuess == antitriton) { ptReco = ptReco - 0.464215 + 0.195771 * ptReco - 0.0183111 * ptReco * ptReco; - // LOG(info) << "we have he3" << pidGuess; } } float ptGen = matchedMCParticle.pt(); + float deltaPt = ptReco - ptGen; if (pdg == -particlePdgCodes.at(4)) { @@ -559,31 +890,37 @@ struct NucleitpcPbPb { histomc.fill(HIST("histDeltaPtVsPtGen"), ptReco, deltaPt); histomc.fill(HIST("histPIDtrack"), ptReco, track.pidForTracking()); } - if (pdg == -particlePdgCodes.at(5)) { histomc.fill(HIST("histDeltaPtVsPtGenHe4anti"), ptReco, deltaPt); + histomc.fill(HIST("histPIDtrackantihe4"), ptReco, track.pidForTracking()); } if (pdg == particlePdgCodes.at(5)) { histomc.fill(HIST("histDeltaPtVsPtGenHe4"), ptReco, deltaPt); + histomc.fill(HIST("histPIDtrackhe4"), ptReco, track.pidForTracking()); } - if (pdg == particlePdgCodes.at(4)) { - histomc.fill(HIST("histPtRecoHe3"), ptReco); - } else if (pdg == -particlePdgCodes.at(4)) { - histomc.fill(HIST("histPtRecoAntiHe3"), ptReco); - } else if (pdg == particlePdgCodes.at(5)) { - histomc.fill(HIST("histPtRecoHe4"), ptReco); - } else if (pdg == -particlePdgCodes.at(5)) { - histomc.fill(HIST("histPtRecoAntiHe4"), ptReco); + // For TOF matching efficiency + float ptTOF = -1.0; // Default: no TOF + if (track.hasTOF()) { + ptTOF = ptReco; } - } + + if (cfgTrackPIDsettings2->get(i, "fillsparsh") == 1) { + histomc.fill(HIST("hSpectramc"), i, sign, collision.centFT0C(), + ptReco, ptTOF); + } + + if (cfgRequirebetaplot) { + histos.fill(HIST("Tofsignal"), getRigidity(track) * track.sign(), o2::pid::tof::Beta::GetBeta(track)); + } + } // loop over primaryParticles types histos.fill(HIST("histeta"), track.eta()); - /*----------------------------------------------------------------------------------------------------------------*/ - } // Track loop + } // TrackAssoc loop } // Collision loop } - PROCESS_SWITCH(NucleitpcPbPb, processMC, "MC reco+gen analysis", false); - //-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= + + PROCESS_SWITCH(NucleitpcPbPb, processMC, "MC reco+gen analysis with efficiency corrections", false); + //=-=-=-==-=-=-==-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= void initCCDB(aod::BCsWithTimestamps::iterator const& bc) { if (mRunNumber == bc.runNumber()) { @@ -639,7 +976,7 @@ struct NucleitpcPbPb { } //---------------------------------------------------------------------------------------------------------------- template - void fillhmass(T const& track, int species) + void fillhmass(T const& track, int species, float cent) { if (!track.hasTOF() || !cfgFillmass) return; @@ -652,15 +989,15 @@ struct NucleitpcPbPb { float massTOF = p * charge * std::sqrt(1.f / (beta * beta) - 1.f); // get PDG mass float pdgMass = particleMasses[species]; - float massDiff = massTOF * massTOF - pdgMass * pdgMass; + float massDiff = massTOF - pdgMass; float ptMomn; setTrackParCov(track, mTrackParCov); mTrackParCov.setPID(track.pidForTracking()); ptMomn = (species == he3 || species == he4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); if (track.sign() > 0) { - hmass[2 * species]->Fill(ptMomn, massDiff); + hmass[2 * species]->Fill(ptMomn, massDiff, cent); } else if (track.sign() < 0) { - hmass[2 * species + 1]->Fill(ptMomn, massDiff); + hmass[2 * species + 1]->Fill(ptMomn, massDiff, cent); } } //---------------------------------------------------------------------------------------------------------------- @@ -674,16 +1011,14 @@ struct NucleitpcPbPb { if (beta < eps || beta > 1.0f - eps) return; float charge = (species == he3 || species == he4) ? 2.f : 1.f; - float p = getRigidity(track); // assuming this is the momentum from inner TPC + float p = getRigidity(track); float massTOF = p * charge * std::sqrt(1.f / (beta * beta) - 1.f); // get PDG mass - float pdgMass = particleMasses[species]; - float massDiff = massTOF * massTOF - pdgMass * pdgMass; - + float masssquare = massTOF * massTOF; if (track.sign() > 0) { - hmassnsigma[2 * species]->Fill(sigma, massDiff); + hmassnsigma[2 * species]->Fill(sigma, masssquare); } else if (track.sign() < 0) { - hmassnsigma[2 * species + 1]->Fill(sigma, massDiff); + hmassnsigma[2 * species + 1]->Fill(sigma, masssquare); } } //---------------------------------------------------------------------------------------------------------------- From 5ac67c1876b403c1a1befcdd64cecc620e380c72 Mon Sep 17 00:00:00 2001 From: Jaideep Tanwar <141036812+jtanwar2212@users.noreply.github.com> Date: Fri, 5 Sep 2025 21:02:21 +0530 Subject: [PATCH 2/8] Update nucleitpcpbpb.cxx --- PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx | 96 ++++++++++------------------ 1 file changed, 33 insertions(+), 63 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx index afcc0a596ee..27b985b8654 100644 --- a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx +++ b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx @@ -150,7 +150,7 @@ struct NucleitpcPbPb { Configurable cfgmaxGetMeanItsClsSizeRequire{"cfgmaxGetMeanItsClsSizeRequire", true, "Require maxGetMeanItsClsSize Cut"}; Configurable cfgDCAwithptRequire{"cfgDCAwithptRequire", true, "Require DCA cuts with pt dependance"}; Configurable cfgRequirebetaplot{"cfgRequirebetaplot", true, "Require beta plot"}; - Configurable cfgIncludeMaterialInEfficiency{"cfgIncludeMaterialInEfficiency", false, "Require from material in efficiency"}; + Configurable cfgIncludeMaterialInEfficiency{"cfgIncludeMaterialInEfficiency", true, "Require from material in efficiency"}; Configurable> cfgBetheBlochParams{"cfgBetheBlochParams", {kBetheBlochDefault[0], nParticles, nBetheParams, particleNames, betheBlochParNames}, "TPC Bethe-Bloch parameterisation for light nuclei"}; Configurable> cfgTrackPIDsettings{"cfgTrackPIDsettings", {kTrackPIDSettings[0], nParticles, nTrkSettings, particleNames, trackPIDsettingsNames}, "track selection and PID criteria"}; @@ -465,8 +465,8 @@ struct NucleitpcPbPb { continue; int pdgCode = mcParticle.pdgCode(); - bool isHe3 = (std::abs(pdgCode) == 1000020030); - bool isHe4 = (std::abs(pdgCode) == 1000020040); + bool isHe3 = (std::abs(pdgCode) == particlePdgCodes.at(4)); + bool isHe4 = (std::abs(pdgCode) == particlePdgCodes.at(5)); if (std::abs(mcParticle.eta()) > cfgCutEta) continue; @@ -493,19 +493,19 @@ struct NucleitpcPbPb { if (isMaterialSecondary) { float centrality = mcCollInfos[idx].passedEvSel ? mcCollInfos[idx].centrality : -1.0f; - if (pdgCode == 1000020030) { + if (pdgCode == particlePdgCodes.at(4)) { histomc.fill(HIST("histSecondaryMaterialHe3"), mcParticle.pt(), centrality); - } else if (pdgCode == -1000020030) { + } else if (pdgCode == -particlePdgCodes.at(4)) { histomc.fill(HIST("histSecondaryMaterialAntiHe3"), mcParticle.pt(), centrality); - } else if (pdgCode == 1000020040) { + } else if (pdgCode == particlePdgCodes.at(5)) { histomc.fill(HIST("histSecondaryMaterialHe4"), mcParticle.pt(), centrality); - } else if (pdgCode == -1000020040) { + } else if (pdgCode == -particlePdgCodes.at(5)) { histomc.fill(HIST("histSecondaryMaterialAntiHe4"), mcParticle.pt(), centrality); } int type = 0; - if (std::abs(pdgCode) == 1000020030) + if (std::abs(pdgCode) == particlePdgCodes.at(4)) type = (pdgCode > 0) ? 1 : 2; - else if (std::abs(pdgCode) == 1000020040) + else if (std::abs(pdgCode) == particlePdgCodes.at(5)) type = (pdgCode > 0) ? 3 : 4; histomc.fill(HIST("histAllMaterialSecondariesGen"), mcParticle.pt(), mcParticle.y(), type); } @@ -514,8 +514,6 @@ struct NucleitpcPbPb { continue; uint16_t flags = 0; - int motherPdgCode = 0; - float motherDecRadius = -1; int decayType = 0; int particleAnti = (pdgCode > 0) ? 0 : 1; @@ -528,9 +526,6 @@ struct NucleitpcPbPb { std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { flags |= kIsSecondaryFromWeakDecay; decayType = 1; - motherPdgCode = motherparticle.pdgCode(); - motherDecRadius = std::hypot(mcParticle.vx() - motherparticle.vx(), - mcParticle.vy() - motherparticle.vy()); break; } } @@ -538,11 +533,6 @@ struct NucleitpcPbPb { } else if (mcParticle.has_mothers()) { flags |= kIsSecondaryFromWeakDecay; decayType = 1; - for (auto& motherparticle : mcParticle.mothers_as()) { - motherPdgCode = motherparticle.pdgCode(); - motherDecRadius = std::hypot(mcParticle.vx() - motherparticle.vx(), - mcParticle.vy() - motherparticle.vy()); - } } else { flags |= kIsSecondaryFromMaterial; decayType = 2; @@ -554,9 +544,9 @@ struct NucleitpcPbPb { continue; int particleType = -1; - if (std::abs(pdgCode) == 1000020030) + if (std::abs(pdgCode) == particlePdgCodes.at(4)) particleType = he3; - else if (std::abs(pdgCode) == 1000020040) + else if (std::abs(pdgCode) == particlePdgCodes.at(5)) particleType = he4; if (particleType >= 0) { @@ -570,23 +560,23 @@ struct NucleitpcPbPb { if (isFromWeakDecay) { float centrality = mcCollInfos[idx].passedEvSel ? mcCollInfos[idx].centrality : -1.0f; - if (pdgCode == 1000020030) { + if (pdgCode == particlePdgCodes.at(4)) { histomc.fill(HIST("histWeakDecayPtHe3"), mcParticle.pt(), centrality); - } else if (pdgCode == -1000020030) { + } else if (pdgCode == -particlePdgCodes.at(4)) { histomc.fill(HIST("histWeakDecayPtAntiHe3"), mcParticle.pt(), centrality); - } else if (pdgCode == 1000020040) { + } else if (pdgCode == particlePdgCodes.at(5)) { histomc.fill(HIST("histWeakDecayPtHe4"), mcParticle.pt(), centrality); - } else if (pdgCode == -1000020040) { + } else if (pdgCode == -particlePdgCodes.at(5)) { histomc.fill(HIST("histWeakDecayPtAntiHe4"), mcParticle.pt(), centrality); } } - if (pdgCode == 1000020030) { + if (pdgCode == particlePdgCodes.at(4)) { histomc.fill(HIST("histPtgenHe3"), mcParticle.pt()); - } else if (pdgCode == -1000020030) { + } else if (pdgCode == -particlePdgCodes.at(4)) { histomc.fill(HIST("histPtgenAntiHe3"), mcParticle.pt()); - } else if (pdgCode == 1000020040) { + } else if (pdgCode == particlePdgCodes.at(5)) { histomc.fill(HIST("histPtgenHe4"), mcParticle.pt()); - } else if (pdgCode == -1000020040) { + } else if (pdgCode == -particlePdgCodes.at(5)) { histomc.fill(HIST("histPtgenAntiHe4"), mcParticle.pt()); } } @@ -609,8 +599,8 @@ struct NucleitpcPbPb { continue; int pdgCode = mcParticle.pdgCode(); - bool isHe3 = (std::abs(pdgCode) == 1000020030); - bool isHe4 = (std::abs(pdgCode) == 1000020040); + bool isHe3 = (std::abs(pdgCode) == particlePdgCodes.at(4)); + bool isHe4 = (std::abs(pdgCode) == particlePdgCodes.at(5)); if (!isHe3 && !isHe4) continue; @@ -620,8 +610,6 @@ struct NucleitpcPbPb { continue; uint16_t flags = 0; - int motherPdgCode = 0; - float motherDecRadius = -1; int decayType = 0; int particleAnti = (pdgCode > 0) ? 0 : 1; @@ -634,9 +622,6 @@ struct NucleitpcPbPb { std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { flags |= kIsSecondaryFromWeakDecay; decayType = 1; - motherPdgCode = motherparticle.pdgCode(); - motherDecRadius = std::hypot(mcParticle.vx() - motherparticle.vx(), - mcParticle.vy() - motherparticle.vy()); break; } } @@ -644,11 +629,6 @@ struct NucleitpcPbPb { } else if (mcParticle.has_mothers()) { flags |= kIsSecondaryFromWeakDecay; decayType = 1; - for (auto& motherparticle : mcParticle.mothers_as()) { - motherPdgCode = motherparticle.pdgCode(); - motherDecRadius = std::hypot(mcParticle.vx() - motherparticle.vx(), - mcParticle.vy() - motherparticle.vy()); - } } else { flags |= kIsSecondaryFromMaterial; decayType = 2; @@ -660,9 +640,9 @@ struct NucleitpcPbPb { continue; int particleType = -1; - if (std::abs(pdgCode) == 1000020030) + if (std::abs(pdgCode) == particlePdgCodes.at(4)) particleType = he3; - else if (std::abs(pdgCode) == 1000020040) + else if (std::abs(pdgCode) == particlePdgCodes.at(5)) particleType = he4; if (particleType >= 0) { @@ -713,8 +693,8 @@ struct NucleitpcPbPb { auto const& matchedMCParticle = particlesMC.iteratorAt(label); int pdg = matchedMCParticle.pdgCode(); - bool isHe3 = (std::abs(pdg) == 1000020030); - bool isHe4 = (std::abs(pdg) == 1000020040); + bool isHe3 = (std::abs(pdg) == particlePdgCodes.at(4)); + bool isHe4 = (std::abs(pdg) == particlePdgCodes.at(5)); bool isMaterialSecondary = false; if ((isHe3 || isHe4) && !matchedMCParticle.isPhysicalPrimary()) { @@ -727,20 +707,20 @@ struct NucleitpcPbPb { } if (isMaterialSecondary) { - if (pdg == 1000020030) { + if (pdg == particlePdgCodes.at(4)) { histomc.fill(HIST("histRecoSecondaryMaterialHe3"), track.pt(), collision.centFT0C()); - } else if (pdg == -1000020030) { + } else if (pdg == -particlePdgCodes.at(4)) { histomc.fill(HIST("histRecoSecondaryMaterialAntiHe3"), track.pt(), collision.centFT0C()); - } else if (pdg == 1000020040) { + } else if (pdg == particlePdgCodes.at(5)) { histomc.fill(HIST("histRecoSecondaryMaterialHe4"), track.pt(), collision.centFT0C()); - } else if (pdg == -1000020040) { + } else if (pdg == -particlePdgCodes.at(5)) { histomc.fill(HIST("histRecoSecondaryMaterialAntiHe4"), track.pt(), collision.centFT0C()); } int type = 0; - if (std::abs(pdg) == 1000020030) + if (std::abs(pdg) == particlePdgCodes.at(4)) type = (pdg > 0) ? 1 : 2; - else if (std::abs(pdg) == 1000020040) + else if (std::abs(pdg) == particlePdgCodes.at(5)) type = (pdg > 0) ? 3 : 4; histomc.fill(HIST("histAllMaterialSecondariesReco"), track.pt(), getRapidity(track, he3), type); @@ -749,8 +729,6 @@ struct NucleitpcPbPb { } } uint16_t flags = 0; - int motherPdgCode = 0; - float motherDecRadius = -1; int decayType = 0; bool isFromWeakDecay = false; @@ -764,9 +742,6 @@ struct NucleitpcPbPb { flags |= kIsSecondaryFromWeakDecay; isFromWeakDecay = true; decayType = 1; - motherPdgCode = motherparticle.pdgCode(); - motherDecRadius = std::hypot(matchedMCParticle.vx() - motherparticle.vx(), - matchedMCParticle.vy() - motherparticle.vy()); break; } } @@ -775,11 +750,6 @@ struct NucleitpcPbPb { flags |= kIsSecondaryFromWeakDecay; isFromWeakDecay = true; decayType = 1; - for (auto& motherparticle : matchedMCParticle.mothers_as()) { - motherPdgCode = motherparticle.pdgCode(); - motherDecRadius = std::hypot(matchedMCParticle.vx() - motherparticle.vx(), - matchedMCParticle.vy() - motherparticle.vy()); - } } else { flags |= kIsSecondaryFromMaterial; decayType = 2; @@ -858,9 +828,9 @@ struct NucleitpcPbPb { histomc.fill(HIST("hNumerEffAcc"), i, matchedMCParticle.pt(), matchedMCParticle.y(), collision.centFT0C(), particleAnti, decayType); } if (isFromWeakDecay) { - if (std::abs(pdg) == 1000020030) { + if (std::abs(pdg) == particlePdgCodes.at(4)) { histomc.fill(HIST("histRecoWeakDecayPtHe3"), ptReco, collision.centFT0C()); - } else if (std::abs(pdg) == 1000020040) { + } else if (std::abs(pdg) == particlePdgCodes.at(5)) { histomc.fill(HIST("histRecoWeakDecayPtHe4"), ptReco, collision.centFT0C()); } } From a60d351cc18b4177f95a146f19d25947a0a74ff3 Mon Sep 17 00:00:00 2001 From: Jaideep Tanwar <141036812+jtanwar2212@users.noreply.github.com> Date: Fri, 5 Sep 2025 21:15:00 +0530 Subject: [PATCH 3/8] Update nucleitpcpbpb.cxx --- PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx index 27b985b8654..446b7688486 100644 --- a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx +++ b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx @@ -521,7 +521,7 @@ struct NucleitpcPbPb { flags |= kIsPhysicalPrimary; decayType = 0; if (mcParticle.has_mothers()) { - for (auto& motherparticle : mcParticle.mothers_as()) { + for (const auto& motherparticle : mcParticle.mothers_as()) { if (std::find(hfMothCodes.begin(), hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { flags |= kIsSecondaryFromWeakDecay; @@ -617,7 +617,7 @@ struct NucleitpcPbPb { flags |= kIsPhysicalPrimary; decayType = 0; if (mcParticle.has_mothers()) { - for (auto& motherparticle : mcParticle.mothers_as()) { + for (const auto& motherparticle : mcParticle.mothers_as()) { if (std::find(hfMothCodes.begin(), hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { flags |= kIsSecondaryFromWeakDecay; @@ -736,7 +736,7 @@ struct NucleitpcPbPb { flags |= kIsPhysicalPrimary; decayType = 0; if (matchedMCParticle.has_mothers()) { - for (auto& motherparticle : matchedMCParticle.mothers_as()) { + for (const auto& motherparticle : matchedMCParticle.mothers_as()) { if (std::find(hfMothCodes.begin(), hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { flags |= kIsSecondaryFromWeakDecay; From ee9f2646f69114b8359f226cd98596bcb0c82e00 Mon Sep 17 00:00:00 2001 From: Jaideep Tanwar <141036812+jtanwar2212@users.noreply.github.com> Date: Fri, 5 Sep 2025 23:09:43 +0530 Subject: [PATCH 4/8] Update nucleitpcpbpb.cxx --- PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx | 37 ++++++++-------------------- 1 file changed, 10 insertions(+), 27 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx index 446b7688486..1eff338cec8 100644 --- a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx +++ b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx @@ -60,11 +60,6 @@ static const std::vector particleMasses{o2::constants::physics::MassPion static const std::vector particleCharge{1, 1, 1, 1, 2, 2}; const int nBetheParams = 6; std::vector hfMothCodes = {511, 521, 531, 541, 5122}; // b-mesons + Lambda_b -enum NucleiFlags { - kIsPhysicalPrimary = BIT(0), - kIsSecondaryFromMaterial = BIT(1), - kIsSecondaryFromWeakDecay = BIT(2) -}; static const std::vector betheBlochParNames{"p0", "p1", "p2", "p3", "p4", "resolution"}; constexpr double kBetheBlochDefault[nParticles][nBetheParams]{ {13.611469, 3.598765, -0.021138, 2.039562, 0.651040, 0.09}, // pion @@ -160,6 +155,7 @@ struct NucleitpcPbPb { Configurable cfgFillmassnsigma{"cfgFillmassnsigma", true, "Fill mass vs nsigma histograms"}; Configurable centcut{"centcut", 80.0f, "centrality cut"}; Configurable cfgCutRapidity{"cfgCutRapidity", 0.5f, "Rapidity range"}; + Configurable cfgtpcNClsFindable{"cfgtpcNClsFindable", 0.8f, "tpcNClsFindable over crossedRows"}; ///////////// Configurable cfgZvertex{"cfgZvertex", 10, "Min Z Vertex"}; Configurable cfgZvertexRequire{"cfgZvertexRequire", true, "Pos Z cut require"}; Configurable cfgZvertexRequireMC{"cfgZvertexRequireMC", true, "Pos Z cut require for generated particles"}; @@ -203,6 +199,7 @@ struct NucleitpcPbPb { TRandom3 rand; float he3 = 4; float he4 = 5; + float processmaterial = 23; //---------------------------------------------------------------------------------------------------------------- void init(InitContext const&) { @@ -361,7 +358,7 @@ struct NucleitpcPbPb { continue; if (track.tpcNClsFound() < cfgTrackPIDsettings->get(i, "minTPCnCls") && cfgTPCNClsfoundRequire) continue; - if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable()) && cfgTPCNClsCrossedRowsRequire) // o2-linter: disable=magic-number (To be checked) + if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || track.tpcNClsCrossedRows() < cfgtpcNClsFindable * track.tpcNClsFindable()) && cfgTPCNClsCrossedRowsRequire) continue; if (track.tpcChi2NCl() > cfgTrackPIDsettings->get(i, "maxTPCchi2") && cfgmaxTPCchi2Require) continue; @@ -478,14 +475,15 @@ struct NucleitpcPbPb { isMaterialSecondary = true; } else { auto mother = mcParticle.mothers_as().front(); - if (std::abs(mother.pdgCode()) < 1000000000) { + float notmother = 1000000000; + if (std::abs(mother.pdgCode()) < notmother) { isMaterialSecondary = true; } } } if (isHe3 || isHe4) { histomc.fill(HIST("histProcessCodeHe"), mcParticle.getProcess()); - if (mcParticle.getProcess() == 23) { + if (mcParticle.getProcess() == processmaterial) { histomc.fill(HIST("histProcess23Details"), std::abs(pdgCode), mcParticle.pt()); } } @@ -513,32 +511,27 @@ struct NucleitpcPbPb { if (!isHe3 && !isHe4) continue; - uint16_t flags = 0; int decayType = 0; int particleAnti = (pdgCode > 0) ? 0 : 1; if (mcParticle.isPhysicalPrimary()) { - flags |= kIsPhysicalPrimary; decayType = 0; if (mcParticle.has_mothers()) { for (const auto& motherparticle : mcParticle.mothers_as()) { if (std::find(hfMothCodes.begin(), hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { - flags |= kIsSecondaryFromWeakDecay; decayType = 1; break; } } } } else if (mcParticle.has_mothers()) { - flags |= kIsSecondaryFromWeakDecay; decayType = 1; } else { - flags |= kIsSecondaryFromMaterial; decayType = 2; continue; } - bool isFromWeakDecay = (flags & kIsSecondaryFromWeakDecay) != 0; + bool isFromWeakDecay = (decayType == 1); if (!mcParticle.isPhysicalPrimary() && !isFromWeakDecay) continue; @@ -609,32 +602,27 @@ struct NucleitpcPbPb { if (std::abs(mcParticle.y()) > cfgCutRapidity) continue; - uint16_t flags = 0; int decayType = 0; int particleAnti = (pdgCode > 0) ? 0 : 1; if (mcParticle.isPhysicalPrimary()) { - flags |= kIsPhysicalPrimary; decayType = 0; if (mcParticle.has_mothers()) { for (const auto& motherparticle : mcParticle.mothers_as()) { if (std::find(hfMothCodes.begin(), hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { - flags |= kIsSecondaryFromWeakDecay; decayType = 1; break; } } } } else if (mcParticle.has_mothers()) { - flags |= kIsSecondaryFromWeakDecay; decayType = 1; } else { - flags |= kIsSecondaryFromMaterial; decayType = 2; continue; } - bool isFromWeakDecay = (flags & kIsSecondaryFromWeakDecay) != 0; + bool isFromWeakDecay = (decayType == 1); if (!mcParticle.isPhysicalPrimary() && !isFromWeakDecay) continue; @@ -699,7 +687,7 @@ struct NucleitpcPbPb { bool isMaterialSecondary = false; if ((isHe3 || isHe4) && !matchedMCParticle.isPhysicalPrimary()) { - if (matchedMCParticle.getProcess() == 23) { // kPMaterial + if (matchedMCParticle.getProcess() == processmaterial) { // kPMaterial isMaterialSecondary = true; } else if (!matchedMCParticle.has_mothers()) { isMaterialSecondary = true; @@ -728,18 +716,15 @@ struct NucleitpcPbPb { continue; } } - uint16_t flags = 0; int decayType = 0; bool isFromWeakDecay = false; if (matchedMCParticle.isPhysicalPrimary()) { - flags |= kIsPhysicalPrimary; decayType = 0; if (matchedMCParticle.has_mothers()) { for (const auto& motherparticle : matchedMCParticle.mothers_as()) { if (std::find(hfMothCodes.begin(), hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { - flags |= kIsSecondaryFromWeakDecay; isFromWeakDecay = true; decayType = 1; break; @@ -747,11 +732,9 @@ struct NucleitpcPbPb { } } } else if (matchedMCParticle.has_mothers()) { - flags |= kIsSecondaryFromWeakDecay; isFromWeakDecay = true; decayType = 1; } else { - flags |= kIsSecondaryFromMaterial; decayType = 2; } @@ -788,7 +771,7 @@ struct NucleitpcPbPb { continue; if (track.tpcNClsFound() < cfgTrackPIDsettings->get(i, "minTPCnCls") && cfgTPCNClsfoundRequire) continue; - if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable()) && cfgTPCNClsCrossedRowsRequire) + if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || track.tpcNClsCrossedRows() < cfgtpcNClsFindable * track.tpcNClsFindable()) && cfgTPCNClsCrossedRowsRequire) continue; if (track.tpcChi2NCl() > cfgTrackPIDsettings->get(i, "maxTPCchi2") && cfgmaxTPCchi2Require) continue; From 869f88f7f2897d082c3ee8fd2507516dfb39a1a0 Mon Sep 17 00:00:00 2001 From: Jaideep Tanwar <141036812+jtanwar2212@users.noreply.github.com> Date: Tue, 16 Sep 2025 13:36:49 +0530 Subject: [PATCH 5/8] Update nucleitpcpbpb.cxx --- PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx | 86 +++++++++++++++++++++------- 1 file changed, 64 insertions(+), 22 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx index 1eff338cec8..37e6608875e 100644 --- a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx +++ b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx @@ -124,7 +124,7 @@ struct NucleitpcPbPb { Configurable cfgCutEta{"cfgCutEta", 0.9f, "Eta range for tracks"}; Configurable cfgetaRequire{"cfgetaRequire", true, "eta cut require"}; Configurable cfgetaRequireMC{"cfgetaRequireMC", true, "eta cut require for generated particles"}; - Configurable cfgrapidityRequireMC{"cfgrapidityRequireMC", true, "rapidity cut require for generated particles"}; + Configurable cfgRapidityRequireMC{"cfgRapidityRequireMC", true, "rapidity cut require for generated particles"}; Configurable cfgUsePVcontributors{"cfgUsePVcontributors", true, "use tracks that are PV contibutors"}; Configurable cfgITSrequire{"cfgITSrequire", true, "Additional cut on ITS require"}; Configurable cfgTPCrequire{"cfgTPCrequire", true, "Additional cut on TPC require"}; @@ -146,6 +146,7 @@ struct NucleitpcPbPb { Configurable cfgDCAwithptRequire{"cfgDCAwithptRequire", true, "Require DCA cuts with pt dependance"}; Configurable cfgRequirebetaplot{"cfgRequirebetaplot", true, "Require beta plot"}; Configurable cfgIncludeMaterialInEfficiency{"cfgIncludeMaterialInEfficiency", true, "Require from material in efficiency"}; + Configurable cfgMasscut{"cfgMasscut", true, "Require mass cut on He4 particles"}; Configurable> cfgBetheBlochParams{"cfgBetheBlochParams", {kBetheBlochDefault[0], nParticles, nBetheParams, particleNames, betheBlochParNames}, "TPC Bethe-Bloch parameterisation for light nuclei"}; Configurable> cfgTrackPIDsettings{"cfgTrackPIDsettings", {kTrackPIDSettings[0], nParticles, nTrkSettings, particleNames, trackPIDsettingsNames}, "track selection and PID criteria"}; @@ -155,7 +156,7 @@ struct NucleitpcPbPb { Configurable cfgFillmassnsigma{"cfgFillmassnsigma", true, "Fill mass vs nsigma histograms"}; Configurable centcut{"centcut", 80.0f, "centrality cut"}; Configurable cfgCutRapidity{"cfgCutRapidity", 0.5f, "Rapidity range"}; - Configurable cfgtpcNClsFindable{"cfgtpcNClsFindable", 0.8f, "tpcNClsFindable over crossedRows"}; ///////////// + Configurable cfgtpcNClsFindable{"cfgtpcNClsFindable", 0.8f, "tpcNClsFindable over crossedRows"}; Configurable cfgZvertex{"cfgZvertex", 10, "Min Z Vertex"}; Configurable cfgZvertexRequire{"cfgZvertexRequire", true, "Pos Z cut require"}; Configurable cfgZvertexRequireMC{"cfgZvertexRequireMC", true, "Pos Z cut require for generated particles"}; @@ -179,7 +180,7 @@ struct NucleitpcPbPb { ConfigurableAxis speciesTrackingAxis{"speciesTrackingAxis", {11, -0.5, 10.5}, "particle type 0: pion, 1: proton, 2: deuteron, 3: triton, 4:He3, 5:He4"}; ConfigurableAxis axisDCA{"axisDCA", {400, -10., 10.}, "DCA axis"}; ConfigurableAxis particleAntiAxis{"particleAntiAxis", {2, 0, 2}, "Particle/Anti-particle"}; // 0 = particle, 1 = anti-particle - ConfigurableAxis decayTypeAxis{"decayTypeAxis", {3, 0, 3}, "Decay type"}; // 0 = primary, 1 = from decay + ConfigurableAxis decayTypeAxis{"decayTypeAxis", {3, -0.5, 2.5}, "Decay type"}; // 0 = primary, 1 = from decay, 2 = material // CCDB Service ccdb; @@ -281,6 +282,8 @@ struct NucleitpcPbPb { histomc.add("histPIDtrackanti", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); histomc.add("histPIDtrackhe4", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); histomc.add("histPIDtrackantihe4", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); + histomc.add("hEventLossDenom", "Event loss denominator", kTH1F, {axisCent}); + histomc.add("hEventLossNumer", "Event loss numerator", kTH1F, {axisCent}); histomc.add("histWeakDecayPtHe3", "Pt distribution of He3 from weak decays", kTH2F, {ptAxis, axisCent}); histomc.add("histWeakDecayPtAntiHe3", "Pt distribution of Anti-He3 from weak decays", kTH2F, {ptAxis, axisCent}); @@ -296,8 +299,6 @@ struct NucleitpcPbPb { histomc.add("histProcess23Details", "Process 23 details", kTH2F, {{4, 0.5, 4.5, "particle type"}, {100, 0, 10, "p_{T}"}}); histomc.add("histAllMaterialSecondariesGen", "All material secondaries (gen)", kTH3F, {{100, 0, 10, "p_{T}"}, {20, -1, 1, "y"}, {5, -0.5, 4.5, "type"}}); histomc.add("histAllMaterialSecondariesReco", "All material secondaries (reco)", kTH3F, {{100, 0, 10, "p_{T}"}, {20, -1, 1, "y"}, {5, -0.5, 4.5, "type"}}); - - // Add axis labels for type: 0=unknown, 1=He3, 2=anti-He3, 3=He4, 4=anti-He4 } } //---------------------------------------------------------------------------------------------------------------- @@ -435,20 +436,50 @@ struct NucleitpcPbPb { mcCollInfos.clear(); mcCollInfos.resize(mcCollisions.size()); + + // Store centrality regardless of cuts FIRST for (auto const& collision : collisions) { int mcCollIdx = collision.mcCollisionId(); if (mcCollIdx < 0 || mcCollIdx >= static_cast(mcCollisions.size())) { continue; } - if (std::abs(collision.posZ()) > cfgZvertex) + // STORE CENTRALITY WITHOUT ANY CUTS + mcCollInfos[mcCollIdx].centrality = collision.centFT0C(); + } + + // FILL DENOMINATOR: ONCE per MC collision + for (size_t i = 0; i < mcCollInfos.size(); i++) { + if (mcCollInfos[i].centrality >= 0) { // Only if we found a matching collision + histomc.fill(HIST("hEventLossDenom"), mcCollInfos[i].centrality); + } + } + + for (auto const& collision : collisions) { + int mcCollIdx = collision.mcCollisionId(); + if (mcCollIdx < 0 || mcCollIdx >= static_cast(mcCollisions.size())) { continue; - if (!collision.sel8()) + } + if (std::abs(collision.posZ()) > cfgZvertex && cfgZvertexRequireMC) + continue; + if (!collision.sel8() && cfgsel8Require) continue; if (collision.centFT0C() > centcut) continue; + // Additional cuts + if (removeITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) + continue; + if (removeNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) + continue; + if (requireIsGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) + continue; + if (requireIsVertexITSTPC && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) + continue; + if (removeNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) + continue; + mcCollInfos[mcCollIdx].passedEvSel = true; - mcCollInfos[mcCollIdx].centrality = collision.centFT0C(); + histomc.fill(HIST("hEventLossNumer"), mcCollInfos[mcCollIdx].centrality); } for (auto const& mcCollision : mcCollisions) { size_t idx = mcCollision.globalIndex(); @@ -464,10 +495,20 @@ struct NucleitpcPbPb { int pdgCode = mcParticle.pdgCode(); bool isHe3 = (std::abs(pdgCode) == particlePdgCodes.at(4)); bool isHe4 = (std::abs(pdgCode) == particlePdgCodes.at(5)); - - if (std::abs(mcParticle.eta()) > cfgCutEta) + if (mcParticle.isPhysicalPrimary()) { + if (pdgCode == particlePdgCodes.at(4)) { + histomc.fill(HIST("histPtgenHe3"), mcParticle.pt()); + } else if (pdgCode == -particlePdgCodes.at(4)) { + histomc.fill(HIST("histPtgenAntiHe3"), mcParticle.pt()); + } else if (pdgCode == particlePdgCodes.at(5)) { + histomc.fill(HIST("histPtgenHe4"), mcParticle.pt()); + } else if (pdgCode == -particlePdgCodes.at(5)) { + histomc.fill(HIST("histPtgenAntiHe4"), mcParticle.pt()); + } + } + if (std::abs(mcParticle.eta()) > cfgCutEta && cfgetaRequireMC) continue; - if (std::abs(mcParticle.y()) > cfgCutRapidity) + if (std::abs(mcParticle.y()) > cfgCutRapidity && cfgRapidityRequireMC) continue; bool isMaterialSecondary = false; if (!mcParticle.isPhysicalPrimary() && (isHe3 || isHe4)) { @@ -543,7 +584,7 @@ struct NucleitpcPbPb { particleType = he4; if (particleType >= 0) { - float centrality = mcCollInfos[idx].passedEvSel ? mcCollInfos[idx].centrality : -1.0f; + float centrality = mcCollInfos[idx].centrality; // Always use actual centrality histomc.fill(HIST("hDenomSignalLoss"), particleType, mcParticle.pt(), mcParticle.y(), centrality, particleAnti, decayType); if (mcCollInfos[idx].passedEvSel) { @@ -563,15 +604,6 @@ struct NucleitpcPbPb { histomc.fill(HIST("histWeakDecayPtAntiHe4"), mcParticle.pt(), centrality); } } - if (pdgCode == particlePdgCodes.at(4)) { - histomc.fill(HIST("histPtgenHe3"), mcParticle.pt()); - } else if (pdgCode == -particlePdgCodes.at(4)) { - histomc.fill(HIST("histPtgenAntiHe3"), mcParticle.pt()); - } else if (pdgCode == particlePdgCodes.at(5)) { - histomc.fill(HIST("histPtgenHe4"), mcParticle.pt()); - } else if (pdgCode == -particlePdgCodes.at(5)) { - histomc.fill(HIST("histPtgenAntiHe4"), mcParticle.pt()); - } } } for (auto const& mcCollision : mcCollisions) { @@ -933,6 +965,7 @@ struct NucleitpcPbPb { { if (!track.hasTOF() || !cfgFillmass) return; + float beta{o2::pid::tof::Beta::GetBeta(track)}; const float eps = 1e-6f; if (beta < eps || beta > 1.0f - eps) @@ -942,7 +975,16 @@ struct NucleitpcPbPb { float massTOF = p * charge * std::sqrt(1.f / (beta * beta) - 1.f); // get PDG mass float pdgMass = particleMasses[species]; - float massDiff = massTOF - pdgMass; + float massDiff = 0.0; + if (species != he4) { + massDiff = massTOF - pdgMass; + } + if (species == he4) { + if ((massTOF * massTOF < 6.5 && massTOF * massTOF < 9.138) && cfgMasscut) + return; + massDiff = massTOF - pdgMass; + } + float ptMomn; setTrackParCov(track, mTrackParCov); mTrackParCov.setPID(track.pidForTracking()); From 4d32f3299bf512701252df1baf6a4f529d088101 Mon Sep 17 00:00:00 2001 From: Jaideep Tanwar <141036812+jtanwar2212@users.noreply.github.com> Date: Tue, 16 Sep 2025 15:12:25 +0530 Subject: [PATCH 6/8] Update nucleitpcpbpb.cxx --- PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx index 37e6608875e..06bfc07ec27 100644 --- a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx +++ b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx @@ -980,7 +980,7 @@ struct NucleitpcPbPb { massDiff = massTOF - pdgMass; } if (species == he4) { - if ((massTOF * massTOF < 6.5 && massTOF * massTOF < 9.138) && cfgMasscut) + if (cfgMasscut && (massTOF * massTOF > 6.5 && massTOF * massTOF < 9.138)) return; massDiff = massTOF - pdgMass; } @@ -1008,8 +1008,19 @@ struct NucleitpcPbPb { float charge = (species == he3 || species == he4) ? 2.f : 1.f; float p = getRigidity(track); float massTOF = p * charge * std::sqrt(1.f / (beta * beta) - 1.f); + // get PDG mass float masssquare = massTOF * massTOF; + + if (species != he4) { + masssquare = massTOF * massTOF; + } + if (species == he4) { + if (cfgMasscut && (massTOF * massTOF > 6.5 && massTOF * massTOF < 9.138)) + return; + masssquare = massTOF * massTOF; + } + if (track.sign() > 0) { hmassnsigma[2 * species]->Fill(sigma, masssquare); } else if (track.sign() < 0) { From 708f5305ea4093f85a974f96bfb92ef62320ab3e Mon Sep 17 00:00:00 2001 From: Jaideep Tanwar <141036812+jtanwar2212@users.noreply.github.com> Date: Wed, 24 Sep 2025 13:56:56 +0530 Subject: [PATCH 7/8] Update nucleitpcpbpb.cxx --- PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx | 725 +++++++++++---------------- 1 file changed, 296 insertions(+), 429 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx index 06bfc07ec27..cadac3bf9e1 100644 --- a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx +++ b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx @@ -109,7 +109,9 @@ std::vector> hmassnsigma; } // namespace //---------------------------------------------------------------------------------------------------------------- struct NucleitpcPbPb { - Preslice perCollision = aod::track_association::collisionId; + + Preslice tracksPerCollision = aod::track::collisionId; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry histomc{"histomc", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; Configurable cfgDebug{"cfgDebug", 1, "debug level"}; @@ -121,7 +123,6 @@ struct NucleitpcPbPb { Configurable removeNoTimeFrameBorder{"removeNoTimeFrameBorder", false, "Remove TF border"}; Configurable cfgRigidityCorrection{"cfgRigidityCorrection", false, "apply rigidity correction"}; // Track Selection Cuts - Configurable cfgCutEta{"cfgCutEta", 0.9f, "Eta range for tracks"}; Configurable cfgetaRequire{"cfgetaRequire", true, "eta cut require"}; Configurable cfgetaRequireMC{"cfgetaRequireMC", true, "eta cut require for generated particles"}; Configurable cfgRapidityRequireMC{"cfgRapidityRequireMC", true, "rapidity cut require for generated particles"}; @@ -143,10 +144,10 @@ struct NucleitpcPbPb { Configurable cfgmaxTPCnSigmaRequire{"cfgmaxTPCnSigmaRequire", true, "Require maxTPCnSigma Cut"}; Configurable cfgminGetMeanItsClsSizeRequire{"cfgminGetMeanItsClsSizeRequire", true, "Require minGetMeanItsClsSize Cut"}; Configurable cfgmaxGetMeanItsClsSizeRequire{"cfgmaxGetMeanItsClsSizeRequire", true, "Require maxGetMeanItsClsSize Cut"}; - Configurable cfgDCAwithptRequire{"cfgDCAwithptRequire", true, "Require DCA cuts with pt dependance"}; Configurable cfgRequirebetaplot{"cfgRequirebetaplot", true, "Require beta plot"}; - Configurable cfgIncludeMaterialInEfficiency{"cfgIncludeMaterialInEfficiency", true, "Require from material in efficiency"}; Configurable cfgMasscut{"cfgMasscut", true, "Require mass cut on He4 particles"}; + Configurable cfgdcaxynopt{"cfgdcaxynopt", true, "DCA xy cut without pT dependent"}; + Configurable cfgdcaznopt{"cfgdcaznopt", false, "DCA xy cut without pT dependent"}; Configurable> cfgBetheBlochParams{"cfgBetheBlochParams", {kBetheBlochDefault[0], nParticles, nBetheParams, particleNames, betheBlochParNames}, "TPC Bethe-Bloch parameterisation for light nuclei"}; Configurable> cfgTrackPIDsettings{"cfgTrackPIDsettings", {kTrackPIDSettings[0], nParticles, nTrkSettings, particleNames, trackPIDsettingsNames}, "track selection and PID criteria"}; @@ -155,16 +156,16 @@ struct NucleitpcPbPb { Configurable cfgFillmass{"cfgFillmass", false, "Fill mass histograms"}; Configurable cfgFillmassnsigma{"cfgFillmassnsigma", true, "Fill mass vs nsigma histograms"}; Configurable centcut{"centcut", 80.0f, "centrality cut"}; + Configurable cfgCutEta{"cfgCutEta", 0.9f, "Eta range for tracks"}; Configurable cfgCutRapidity{"cfgCutRapidity", 0.5f, "Rapidity range"}; Configurable cfgtpcNClsFindable{"cfgtpcNClsFindable", 0.8f, "tpcNClsFindable over crossedRows"}; Configurable cfgZvertex{"cfgZvertex", 10, "Min Z Vertex"}; - Configurable cfgZvertexRequire{"cfgZvertexRequire", true, "Pos Z cut require"}; - Configurable cfgZvertexRequireMC{"cfgZvertexRequireMC", true, "Pos Z cut require for generated particles"}; + Configurable cfgZvertexRequireMC{"cfgZvertexRequireMC", true, "Pos Z cut in MC"}; Configurable cfgsel8Require{"cfgsel8Require", true, "sel8 cut require"}; o2::track::TrackParametrizationWithError mTrackParCov; // Binning configuration ConfigurableAxis axisMagField{"axisMagField", {10, -10., 10.}, "magnetic field"}; - ConfigurableAxis axisNev{"axisNev", {5, 0., 5.}, "Number of events"}; + ConfigurableAxis axisNev{"axisNev", {10, 0., 10.}, "Number of events"}; ConfigurableAxis axisRigidity{"axisRigidity", {4000, -10., 10.}, "#it{p}^{TPC}/#it{z}"}; ConfigurableAxis axisdEdx{"axisdEdx", {4000, 0, 4000}, "d#it{E}/d#it{x}"}; ConfigurableAxis axisCent{"axisCent", {100, 0, 100}, "centrality"}; @@ -174,13 +175,13 @@ struct NucleitpcPbPb { ConfigurableAxis axiseta{"axiseta", {100, -1, 1}, "eta"}; ConfigurableAxis axisrapidity{"axisrapidity", {100, -2, 2}, "rapidity"}; ConfigurableAxis axismass{"axismass", {100, -10, 10}, "mass"}; - ConfigurableAxis axismassnsigma{"axismassnsigma", {100, 0, 20}, "mass"}; + ConfigurableAxis axismassnsigma{"axismassnsigma", {100, 0, 20}, "nsigma mass"}; ConfigurableAxis nsigmaAxis{"nsigmaAxis", {160, -10, 10}, "n#sigma_{#pi^{+}}"}; ConfigurableAxis speciesBitAxis{"speciesBitAxis", {8, -0.5, 7.5}, "particle type 0: pion, 1: proton, 2: deuteron, 3: triton, 4:He3, 5:He4"}; ConfigurableAxis speciesTrackingAxis{"speciesTrackingAxis", {11, -0.5, 10.5}, "particle type 0: pion, 1: proton, 2: deuteron, 3: triton, 4:He3, 5:He4"}; ConfigurableAxis axisDCA{"axisDCA", {400, -10., 10.}, "DCA axis"}; - ConfigurableAxis particleAntiAxis{"particleAntiAxis", {2, 0, 2}, "Particle/Anti-particle"}; // 0 = particle, 1 = anti-particle - ConfigurableAxis decayTypeAxis{"decayTypeAxis", {3, -0.5, 2.5}, "Decay type"}; // 0 = primary, 1 = from decay, 2 = material + ConfigurableAxis particleAntiAxis{"particleAntiAxis", {2, -0.5, 1.5}, "Particle/Anti-particle"}; // 0 = particle, 1 = anti-particle + ConfigurableAxis decayTypeAxis{"decayTypeAxis", {3, -0.5, 2.5}, "Decay type"}; // 0 = primary, 1 = from decay, 2 = material // CCDB Service ccdb; @@ -194,13 +195,12 @@ struct NucleitpcPbPb { std::vector primaryParticles; std::vector primVtx, cents; - bool collHasCandidate, collPassedEvSel, collPassedEvSelMc; + bool collHasCandidate, collPassedEvSel; int mRunNumber, occupancy; - float dBz, momn; + float dBz; TRandom3 rand; float he3 = 4; float he4 = 5; - float processmaterial = 23; //---------------------------------------------------------------------------------------------------------------- void init(InitContext const&) { @@ -249,27 +249,28 @@ struct NucleitpcPbPb { if (doprocessMC) { histomc.add("hSpectramc", " ", HistType::kTHnSparseF, {speciesBitAxis, {5, -2.5, 2.5}, axisCent, ptAxis, ptAxis}); + // Efficiency x Acceptance histomc.add("hDenomEffAcc", "Denominator for Efficiency x Acceptance", {HistType::kTHnSparseF, {speciesBitAxis, ptAxis, axisrapidity, axisCent, particleAntiAxis, decayTypeAxis}}); histomc.add("hNumerEffAcc", "Numerator for Efficiency x Acceptance", {HistType::kTHnSparseF, {speciesBitAxis, ptAxis, axisrapidity, axisCent, particleAntiAxis, decayTypeAxis}}); - histomc.add("hDenomSignalLoss", "Denominator for Signal Loss", - {HistType::kTHnSparseF, {speciesBitAxis, ptAxis, axisrapidity, axisCent, particleAntiAxis, decayTypeAxis}}); - histomc.add("hNumerSignalLoss", "Numerator for Signal Loss", - {HistType::kTHnSparseF, {speciesBitAxis, ptAxis, axisrapidity, axisCent, particleAntiAxis, decayTypeAxis}}); - // Material secondary histograms - histomc.add("histSecondaryMaterialHe3", "He3 from material", HistType::kTH2F, {ptAxis, axisCent}); - histomc.add("histSecondaryMaterialAntiHe3", "Anti-He3 from material", HistType::kTH2F, {ptAxis, axisCent}); - histomc.add("histSecondaryMaterialHe4", "He4 from material", HistType::kTH2F, {ptAxis, axisCent}); - histomc.add("histSecondaryMaterialAntiHe4", "Anti-He4 from material", HistType::kTH2F, {ptAxis, axisCent}); + // The Signal loss correction + histomc.add("hHe3SignalLossDenom", "He3 Signal Loss Denominator", kTH1F, {axisCent}); + histomc.add("hHe3SignalLossNumer", "He3 Signal Loss Numerator", kTH1F, {axisCent}); + histomc.add("hHe4SignalLossDenom", "He4 Signal Loss Denominator", kTH1F, {axisCent}); + histomc.add("hHe4SignalLossNumer", "He4 Signal Loss Numerator", kTH1F, {axisCent}); + + histomc.add("haHe3SignalLossDenom", "He3 Signal Loss Denominator", kTH1F, {axisCent}); + histomc.add("haHe3SignalLossNumer", "He3 Signal Loss Numerator", kTH1F, {axisCent}); + histomc.add("haHe4SignalLossDenom", "He4 Signal Loss Denominator", kTH1F, {axisCent}); + histomc.add("haHe4SignalLossNumer", "He4 Signal Loss Numerator", kTH1F, {axisCent}); + + // The event loss correction + histomc.add("hEventLossDenom", "Event loss denominator", kTH1F, {axisCent}); + histomc.add("hEventLossNumer", "Event loss numerator", kTH1F, {axisCent}); histomc.add("histVtxZgen", "histVtxZgen", kTH1F, {axisVtxZ}); - histomc.add("histEtagen", "histEtagen", kTH1F, {axiseta}); - histomc.add("histPtgenHe3", "histPtgenHe3", kTH1F, {ptAxis}); - histomc.add("histPtgenAntiHe3", "histPtgenAntiHe3", kTH1F, {ptAxis}); - histomc.add("histPtgenHe4", "histPtgenHe4", kTH1F, {ptAxis}); - histomc.add("histPtgenAntiHe4", "histPtgenAntiHe4", kTH1F, {ptAxis}); histomc.add("histNevReco", "histNevReco", kTH1F, {axisNev}); histomc.add("histVtxZReco", "histVtxZReco", kTH1F, {axisVtxZ}); histomc.add("histCentFT0CReco", "histCentFT0CReco", kTH1F, {axisCent}); @@ -282,54 +283,47 @@ struct NucleitpcPbPb { histomc.add("histPIDtrackanti", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); histomc.add("histPIDtrackhe4", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); histomc.add("histPIDtrackantihe4", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); - histomc.add("hEventLossDenom", "Event loss denominator", kTH1F, {axisCent}); - histomc.add("hEventLossNumer", "Event loss numerator", kTH1F, {axisCent}); - - histomc.add("histWeakDecayPtHe3", "Pt distribution of He3 from weak decays", kTH2F, {ptAxis, axisCent}); - histomc.add("histWeakDecayPtAntiHe3", "Pt distribution of Anti-He3 from weak decays", kTH2F, {ptAxis, axisCent}); - histomc.add("histWeakDecayPtHe4", "Pt distribution of He4 from weak decays", kTH2F, {ptAxis, axisCent}); - histomc.add("histWeakDecayPtAntiHe4", "Pt distribution of Anti-He4 from weak decays", kTH2F, {ptAxis, axisCent}); - histomc.add("histRecoWeakDecayPtHe3", "Pt distribution of reconstructed He3 from weak decays", kTH2F, {ptAxis, axisCent}); - histomc.add("histRecoWeakDecayPtHe4", "Pt distribution of reconstructed He4 from weak decays", kTH2F, {ptAxis, axisCent}); - histomc.add("histRecoSecondaryMaterialHe3", "Reco He3 from material", HistType::kTH2F, {ptAxis, axisCent}); - histomc.add("histRecoSecondaryMaterialAntiHe3", "Reco Anti-He3 from material", HistType::kTH2F, {ptAxis, axisCent}); - histomc.add("histRecoSecondaryMaterialHe4", "Reco He4 from material", HistType::kTH2F, {ptAxis, axisCent}); - histomc.add("histRecoSecondaryMaterialAntiHe4", "Reco Anti-He4 from material", HistType::kTH2F, {ptAxis, axisCent}); - histomc.add("histProcessCodeHe", "Process codes for He isotopes", kTH1F, {{100, 0, 100, "process code"}}); - histomc.add("histProcess23Details", "Process 23 details", kTH2F, {{4, 0.5, 4.5, "particle type"}, {100, 0, 10, "p_{T}"}}); - histomc.add("histAllMaterialSecondariesGen", "All material secondaries (gen)", kTH3F, {{100, 0, 10, "p_{T}"}, {20, -1, 1, "y"}, {5, -0.5, 4.5, "type"}}); - histomc.add("histAllMaterialSecondariesReco", "All material secondaries (reco)", kTH3F, {{100, 0, 10, "p_{T}"}, {20, -1, 1, "y"}, {5, -0.5, 4.5, "type"}}); } } //---------------------------------------------------------------------------------------------------------------- //---------------------------------------------------------------------------------------------------------------- - void processData(CollisionsFull const& collisions, TracksFull const& tracks, aod::BCsWithTimestamps const&, aod::TrackAssoc const& tracksColl) + void processData(CollisionsFull const& collisions, + TracksFull const& tracks, + aod::BCsWithTimestamps const&) { for (const auto& collision : collisions) { auto bc = collision.bc_as(); initCCDB(bc); initCollision(collision); + if (!collPassedEvSel) continue; if (collision.centFT0C() > centcut) continue; if (removeITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) continue; + histos.fill(HIST("histNev"), 2.5); if (removeNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) continue; + histos.fill(HIST("histNev"), 3.5); + if (requireIsGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) continue; + histos.fill(HIST("histNev"), 4.5); if (requireIsVertexITSTPC && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) continue; + + histos.fill(HIST("histNev"), 5.5); if (removeNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) continue; + histos.fill(HIST("histNev"), 6.5); histos.fill(HIST("histCentFTOC_cut"), collision.centFT0C()); - histos.fill(HIST("histNev"), 2.5); - const uint64_t collIdx = collision.globalIndex(); - auto tracksByColl = tracksColl.sliceBy(perCollision, collIdx); - ///////////////////////////////////////////////////////////////////////////////// - for (const auto& trackId : tracksByColl) { - const auto& track = tracks.rawIteratorAt(trackId.trackId()); + + // new slicing + auto tracksInColl = tracks.sliceBy(tracksPerCollision, collision.globalIndex()); + + // loop over sliced tracks + for (const auto& track : tracksInColl) { if (!track.isPVContributor() && cfgUsePVcontributors) continue; if (!track.hasITS() && cfgITSrequire) @@ -342,24 +336,22 @@ struct NucleitpcPbPb { continue; if (std::abs(track.eta()) > cfgCutEta && cfgetaRequire) continue; - for (size_t i = 0; i < primaryParticles.size(); i++) { + for (size_t i = 0; i < primaryParticles.size(); i++) { float ptMomn; setTrackParCov(track, mTrackParCov); mTrackParCov.setPID(track.pidForTracking()); ptMomn = (i == he3 || i == he4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); - int sign = 0; - if (track.sign() > 0) { - sign = 1; - } - if (track.sign() < 0) { - sign = -1; - } + + int sign = (track.sign() > 0) ? 1 : ((track.sign() < 0) ? -1 : 0); + if (std::abs(getRapidity(track, i)) > cfgCutRapidity && cfgRapidityRequire) continue; if (track.tpcNClsFound() < cfgTrackPIDsettings->get(i, "minTPCnCls") && cfgTPCNClsfoundRequire) continue; - if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || track.tpcNClsCrossedRows() < cfgtpcNClsFindable * track.tpcNClsFindable()) && cfgTPCNClsCrossedRowsRequire) + if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || + track.tpcNClsCrossedRows() < cfgtpcNClsFindable * track.tpcNClsFindable()) && + cfgTPCNClsCrossedRowsRequire) continue; if (track.tpcChi2NCl() > cfgTrackPIDsettings->get(i, "maxTPCchi2") && cfgmaxTPCchi2Require) continue; @@ -367,6 +359,7 @@ struct NucleitpcPbPb { continue; if (track.itsNCls() < cfgTrackPIDsettings->get(i, "minITSnCls") && cfgminITSnClsRequire) continue; + double cosheta = std::cosh(track.eta()); if ((track.itsNCls() / cosheta) < cfgTrackPIDsettings->get(i, "minITSnClscos") && cfgminITSnClscosRequire) continue; @@ -377,9 +370,15 @@ struct NucleitpcPbPb { if (getMeanItsClsSize(track) < cfgTrackPIDsettings->get(i, "minITSclsSize") && cfgminGetMeanItsClsSizeRequire) continue; - bool insideDCAxy = (std::abs(track.dcaXY()) <= (cfgTrackPIDsettings->get(i, "maxDcaXY") * (0.0105f + 0.0350f / std::pow(ptMomn, 1.1f)))); // o2-linter: disable=magic-number (To be checked) - if ((!(insideDCAxy) || std::abs(track.dcaZ()) > dcazSigma(ptMomn, cfgTrackPIDsettings->get(i, "maxDcaZ"))) && cfgDCAwithptRequire) + // DCA XY cut + bool insideDCAxy = cfgdcaxynopt ? (std::abs(track.dcaXY()) <= cfgTrackPIDsettings->get(i, "maxDcaXY")) : (std::abs(track.dcaXY()) <= (cfgTrackPIDsettings->get(i, "maxDcaXY") * (0.0105f + 0.0350f / std::pow(ptMomn, 1.1f)))); + + // DCA Z cut + bool insideDCAz = cfgdcaznopt ? (std::abs(track.dcaZ()) <= cfgTrackPIDsettings->get(i, "maxDcaZ")) : (std::abs(track.dcaZ()) <= dcazSigma(ptMomn, cfgTrackPIDsettings->get(i, "maxDcaZ"))); + + if ((!insideDCAxy || !insideDCAz)) { continue; + } float tpcNsigma = getTPCnSigma(track, primaryParticles.at(i)); if ((std::abs(tpcNsigma) > cfgTrackPIDsettings->get(i, "maxTPCnSigma")) && cfgmaxTPCnSigmaRequire) @@ -389,10 +388,36 @@ struct NucleitpcPbPb { continue; if (itsSigma > cfgTrackPIDsettings2->get(i, "maxITSnsigma") && cfgTrackPIDsettings2->get(i, "useITSnsigma") < 1) continue; + histos.fill(HIST("Tpcsignal"), getRigidity(track) * track.sign(), track.tpcSignal()); if (cfgFillhspectra && cfgTrackPIDsettings2->get(i, "fillsparsh") == 1) { - histos.fill(HIST("hSpectra"), i, ptMomn, tpcNsigma, sign, collision.centFT0C(), track.dcaZ(), track.dcaXY()); + + if (i != he4) { + histos.fill(HIST("hSpectra"), i, ptMomn, tpcNsigma, sign, collision.centFT0C(), track.dcaZ(), track.dcaXY()); + } else { + if (!track.hasTOF()) { + // Fill without TOF + histos.fill(HIST("hSpectra"), i, ptMomn, tpcNsigma, sign, collision.centFT0C(), track.dcaZ(), track.dcaXY()); + } else { + // Has TOF - apply mass cut + float beta = o2::pid::tof::Beta::GetBeta(track); + const float eps = 1e-6f; + if (beta < eps || beta > 1.0f - eps) + continue; + + float charge = 2.f; // he4 has charge 2 + float p = getRigidity(track); + float massTOF = p * charge * std::sqrt(1.f / (beta * beta) - 1.f); + + // Apply mass cut for he4 (mass^2 around 3.73^2 = 13.9) + if (cfgMasscut && (massTOF * massTOF > 6.5 && massTOF * massTOF < 9.138)) { + continue; // Skip if mass cut fails + } + + histos.fill(HIST("hSpectra"), i, ptMomn, tpcNsigma, sign, collision.centFT0C(), track.dcaZ(), track.dcaXY()); + } + } } fillhmassnsigma(track, i, tpcNsigma); if ((std::abs(tpcNsigma) > cfgTrackPIDsettings2->get(i, "maxTPCnsigmaTOF")) && cfgTrackPIDsettings2->get(i, "useTPCnsigmaTOF") < 1) @@ -402,11 +427,11 @@ struct NucleitpcPbPb { if (cfgRequirebetaplot) { histos.fill(HIST("Tofsignal"), getRigidity(track) * track.sign(), o2::pid::tof::Beta::GetBeta(track)); } - } + } // loop primaryParticles + histos.fill(HIST("histeta"), track.eta()); - } // track loop - /////////////////////////////////////////////// - } + } // loop sliced tracks + } // collision loop } PROCESS_SWITCH(NucleitpcPbPb, processData, "data analysis", false); @@ -416,51 +441,31 @@ struct NucleitpcPbPb { struct McCollInfo { bool passedEvSel = false; float centrality = -1.0f; + bool passedEvSelVtZ = false; }; std::vector mcCollInfos; void processMC(CollisionsFullMC const& collisions, aod::McCollisions const& mcCollisions, - TracksFull const& tracks, - aod::BCsWithTimestamps const& bcs, + soa::Join const& tracks, aod::McParticles const& particlesMC, - aod::McTrackLabels const& trackLabelsMC, - aod::McCollisionLabels const& collLabels, - aod::TrackAssoc const& tracksColl, - CollisionsFull const& colls) - { + aod::BCsWithTimestamps const&) - (void)bcs; - (void)collLabels; - (void)colls; + { mcCollInfos.clear(); mcCollInfos.resize(mcCollisions.size()); - // Store centrality regardless of cuts FIRST + // First pass: Store centrality and apply event selection for (auto const& collision : collisions) { int mcCollIdx = collision.mcCollisionId(); if (mcCollIdx < 0 || mcCollIdx >= static_cast(mcCollisions.size())) { continue; } - // STORE CENTRALITY WITHOUT ANY CUTS - mcCollInfos[mcCollIdx].centrality = collision.centFT0C(); - } - // FILL DENOMINATOR: ONCE per MC collision - for (size_t i = 0; i < mcCollInfos.size(); i++) { - if (mcCollInfos[i].centrality >= 0) { // Only if we found a matching collision - histomc.fill(HIST("hEventLossDenom"), mcCollInfos[i].centrality); - } - } + // STORE CENTRALITY WITHOUt CUTS + mcCollInfos[mcCollIdx].centrality = collision.centFT0C(); - for (auto const& collision : collisions) { - int mcCollIdx = collision.mcCollisionId(); - if (mcCollIdx < 0 || mcCollIdx >= static_cast(mcCollisions.size())) { - continue; - } - if (std::abs(collision.posZ()) > cfgZvertex && cfgZvertexRequireMC) - continue; if (!collision.sel8() && cfgsel8Require) continue; if (collision.centFT0C() > centcut) @@ -478,146 +483,75 @@ struct NucleitpcPbPb { if (removeNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) continue; + // Mark this MC collision as passing event selection mcCollInfos[mcCollIdx].passedEvSel = true; - histomc.fill(HIST("hEventLossNumer"), mcCollInfos[mcCollIdx].centrality); - } - for (auto const& mcCollision : mcCollisions) { - size_t idx = mcCollision.globalIndex(); - if (idx >= mcCollInfos.size()) + // Apply event selection cuts + if (std::abs(collision.posZ()) > cfgZvertex && cfgZvertexRequireMC) continue; - for (auto const& mcParticle : particlesMC) { + mcCollInfos[mcCollIdx].passedEvSelVtZ = true; + } - if (mcParticle.mcCollisionId() != mcCollision.globalIndex()) - continue; + // FILL EVENT LOSS AND SIGNAL LOSS: Combined loop per MC collision + for (size_t i = 0; i < mcCollInfos.size(); i++) { + if (mcCollInfos[i].centrality >= 0) { // Only if we found a matching collision + // Event loss denominator + histomc.fill(HIST("hEventLossDenom"), mcCollInfos[i].centrality); - int pdgCode = mcParticle.pdgCode(); - bool isHe3 = (std::abs(pdgCode) == particlePdgCodes.at(4)); - bool isHe4 = (std::abs(pdgCode) == particlePdgCodes.at(5)); - if (mcParticle.isPhysicalPrimary()) { - if (pdgCode == particlePdgCodes.at(4)) { - histomc.fill(HIST("histPtgenHe3"), mcParticle.pt()); - } else if (pdgCode == -particlePdgCodes.at(4)) { - histomc.fill(HIST("histPtgenAntiHe3"), mcParticle.pt()); - } else if (pdgCode == particlePdgCodes.at(5)) { - histomc.fill(HIST("histPtgenHe4"), mcParticle.pt()); - } else if (pdgCode == -particlePdgCodes.at(5)) { - histomc.fill(HIST("histPtgenAntiHe4"), mcParticle.pt()); - } - } - if (std::abs(mcParticle.eta()) > cfgCutEta && cfgetaRequireMC) - continue; - if (std::abs(mcParticle.y()) > cfgCutRapidity && cfgRapidityRequireMC) - continue; - bool isMaterialSecondary = false; - if (!mcParticle.isPhysicalPrimary() && (isHe3 || isHe4)) { - if (!mcParticle.has_mothers()) { - isMaterialSecondary = true; - } else { - auto mother = mcParticle.mothers_as().front(); - float notmother = 1000000000; - if (std::abs(mother.pdgCode()) < notmother) { - isMaterialSecondary = true; - } - } - } - if (isHe3 || isHe4) { - histomc.fill(HIST("histProcessCodeHe"), mcParticle.getProcess()); - if (mcParticle.getProcess() == processmaterial) { - histomc.fill(HIST("histProcess23Details"), std::abs(pdgCode), mcParticle.pt()); - } + // Event loss numerator (if passed selection) + if (mcCollInfos[i].passedEvSel) { + histomc.fill(HIST("hEventLossNumer"), mcCollInfos[i].centrality); } - if (isMaterialSecondary) { - float centrality = mcCollInfos[idx].passedEvSel ? mcCollInfos[idx].centrality : -1.0f; - - if (pdgCode == particlePdgCodes.at(4)) { - histomc.fill(HIST("histSecondaryMaterialHe3"), mcParticle.pt(), centrality); - } else if (pdgCode == -particlePdgCodes.at(4)) { - histomc.fill(HIST("histSecondaryMaterialAntiHe3"), mcParticle.pt(), centrality); - } else if (pdgCode == particlePdgCodes.at(5)) { - histomc.fill(HIST("histSecondaryMaterialHe4"), mcParticle.pt(), centrality); - } else if (pdgCode == -particlePdgCodes.at(5)) { - histomc.fill(HIST("histSecondaryMaterialAntiHe4"), mcParticle.pt(), centrality); + // Fill signal loss for all primary particles in this MC collision + for (auto const& mcParticle : particlesMC) { + if (mcParticle.mcCollisionId() != static_cast(i)) { + continue; } - int type = 0; - if (std::abs(pdgCode) == particlePdgCodes.at(4)) - type = (pdgCode > 0) ? 1 : 2; - else if (std::abs(pdgCode) == particlePdgCodes.at(5)) - type = (pdgCode > 0) ? 3 : 4; - histomc.fill(HIST("histAllMaterialSecondariesGen"), mcParticle.pt(), mcParticle.y(), type); - } - - if (!isHe3 && !isHe4) - continue; - - int decayType = 0; - int particleAnti = (pdgCode > 0) ? 0 : 1; - - if (mcParticle.isPhysicalPrimary()) { - decayType = 0; - if (mcParticle.has_mothers()) { - for (const auto& motherparticle : mcParticle.mothers_as()) { - if (std::find(hfMothCodes.begin(), hfMothCodes.end(), - std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { - decayType = 1; - break; - } - } + if (!mcParticle.isPhysicalPrimary()) { + continue; } - } else if (mcParticle.has_mothers()) { - decayType = 1; - } else { - decayType = 2; - continue; - } - bool isFromWeakDecay = (decayType == 1); - - if (!mcParticle.isPhysicalPrimary() && !isFromWeakDecay) - continue; - - int particleType = -1; - if (std::abs(pdgCode) == particlePdgCodes.at(4)) - particleType = he3; - else if (std::abs(pdgCode) == particlePdgCodes.at(5)) - particleType = he4; - - if (particleType >= 0) { - float centrality = mcCollInfos[idx].centrality; // Always use actual centrality - histomc.fill(HIST("hDenomSignalLoss"), particleType, mcParticle.pt(), mcParticle.y(), centrality, particleAnti, decayType); - if (mcCollInfos[idx].passedEvSel) { - histomc.fill(HIST("hNumerSignalLoss"), particleType, mcParticle.pt(), mcParticle.y(), centrality, particleAnti, decayType); + // Signal loss denominator + if (mcParticle.pdgCode() == particlePdgCodes.at(4)) { // He3 + histomc.fill(HIST("hHe3SignalLossDenom"), mcCollInfos[i].centrality); + } else if (mcParticle.pdgCode() == particlePdgCodes.at(5)) { // He4 + histomc.fill(HIST("hHe4SignalLossDenom"), mcCollInfos[i].centrality); + } else if (mcParticle.pdgCode() == -particlePdgCodes.at(4)) { // anti-He3 + histomc.fill(HIST("haHe3SignalLossDenom"), mcCollInfos[i].centrality); + } else if (mcParticle.pdgCode() == -particlePdgCodes.at(5)) { // He4 + histomc.fill(HIST("haHe4SignalLossDenom"), mcCollInfos[i].centrality); } - } - if (isFromWeakDecay) { - float centrality = mcCollInfos[idx].passedEvSel ? mcCollInfos[idx].centrality : -1.0f; - if (pdgCode == particlePdgCodes.at(4)) { - histomc.fill(HIST("histWeakDecayPtHe3"), mcParticle.pt(), centrality); - } else if (pdgCode == -particlePdgCodes.at(4)) { - histomc.fill(HIST("histWeakDecayPtAntiHe3"), mcParticle.pt(), centrality); - } else if (pdgCode == particlePdgCodes.at(5)) { - histomc.fill(HIST("histWeakDecayPtHe4"), mcParticle.pt(), centrality); - } else if (pdgCode == -particlePdgCodes.at(5)) { - histomc.fill(HIST("histWeakDecayPtAntiHe4"), mcParticle.pt(), centrality); + // Signal loss numerator (if event passed selection) + if (mcCollInfos[i].passedEvSel) { + if (mcParticle.pdgCode() == particlePdgCodes.at(4)) { // He3 + histomc.fill(HIST("hHe3SignalLossNumer"), mcCollInfos[i].centrality); + } else if (mcParticle.pdgCode() == particlePdgCodes.at(5)) { // He4 + histomc.fill(HIST("hHe4SignalLossNumer"), mcCollInfos[i].centrality); + } else if (mcParticle.pdgCode() == -particlePdgCodes.at(4)) { // anti-He3 + histomc.fill(HIST("haHe3SignalLossNumer"), mcCollInfos[i].centrality); + } else if (mcParticle.pdgCode() == -particlePdgCodes.at(5)) { // anti-He4 + histomc.fill(HIST("haHe4SignalLossNumer"), mcCollInfos[i].centrality); + } } } } } + + // Process MC collisions for efficiency and reconstructed collisions for (auto const& mcCollision : mcCollisions) { size_t idx = mcCollision.globalIndex(); - if (idx >= mcCollInfos.size()) continue; - if (!mcCollInfos[idx].passedEvSel) - continue; - if (std::abs(mcCollision.posZ()) > cfgZvertex) - continue; - histomc.fill(HIST("histVtxZgen"), mcCollision.posZ()); float centrality = mcCollInfos[idx].centrality; + // bool passedEvSel = mcCollInfos[idx].passedEvSel; + bool passedEvSelVtZ = mcCollInfos[idx].passedEvSelVtZ; + + // Process generated particles for efficiency denominators + histomc.fill(HIST("histVtxZgen"), mcCollision.posZ()); for (auto const& mcParticle : particlesMC) { if (mcParticle.mcCollisionId() != mcCollision.globalIndex()) @@ -629,9 +563,10 @@ struct NucleitpcPbPb { if (!isHe3 && !isHe4) continue; - if (std::abs(mcParticle.eta()) > cfgCutEta) + + if (std::abs(mcParticle.eta()) > cfgCutEta && cfgetaRequireMC) continue; - if (std::abs(mcParticle.y()) > cfgCutRapidity) + if (std::abs(mcParticle.y()) > cfgCutRapidity && cfgRapidityRequireMC) continue; int decayType = 0; @@ -656,7 +591,8 @@ struct NucleitpcPbPb { } bool isFromWeakDecay = (decayType == 1); - if (!mcParticle.isPhysicalPrimary() && !isFromWeakDecay) + // if (!mcParticle.isPhysicalPrimary() && !isFromWeakDecay) + if (!mcParticle.isPhysicalPrimary() && isFromWeakDecay) continue; int particleType = -1; @@ -666,244 +602,175 @@ struct NucleitpcPbPb { particleType = he4; if (particleType >= 0) { - histomc.fill(HIST("hDenomEffAcc"), particleType, mcParticle.pt(), mcParticle.y(), centrality, particleAnti, decayType); - } - } - } - for (auto const& collision : collisions) { - auto mcCollIdx = collision.mcCollisionId(); - if (mcCollIdx < 0 || mcCollIdx >= static_cast(mcCollisions.size())) - continue; - - auto bc = collision.bc_as(); - initCCDB(bc); - collHasCandidate = false; - histomc.fill(HIST("histNevReco"), 0.5); - - if (std::abs(collision.posZ()) > cfgZvertex && cfgZvertexRequire) - continue; - collPassedEvSel = collision.sel8(); - if (!collPassedEvSel && cfgsel8Require) - continue; - histomc.fill(HIST("histNevReco"), 1.5); - histomc.fill(HIST("histVtxZReco"), collision.posZ()); - histomc.fill(HIST("histCentFT0CReco"), collision.centFT0C()); - histomc.fill(HIST("histCentFT0MReco"), collision.centFT0M()); - if (collision.centFT0C() > centcut) - continue; - histomc.fill(HIST("histNevReco"), 2.5); - if (removeITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) - continue; - if (removeNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) - continue; - if (requireIsGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) - continue; - if (requireIsVertexITSTPC && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) - continue; - if (removeNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) - continue; - for (const auto& assoc : tracksColl) { - if (assoc.collisionId() != collision.globalIndex()) - continue; - const auto& track = tracks.rawIteratorAt(assoc.trackId()); - auto labelRow = trackLabelsMC.iteratorAt(track.globalIndex()); - int label = labelRow.mcParticleId(); - if (label < 0 || label >= static_cast(particlesMC.size())) - continue; - auto const& matchedMCParticle = particlesMC.iteratorAt(label); - - int pdg = matchedMCParticle.pdgCode(); - bool isHe3 = (std::abs(pdg) == particlePdgCodes.at(4)); - bool isHe4 = (std::abs(pdg) == particlePdgCodes.at(5)); - bool isMaterialSecondary = false; - if ((isHe3 || isHe4) && !matchedMCParticle.isPhysicalPrimary()) { - - if (matchedMCParticle.getProcess() == processmaterial) { // kPMaterial - isMaterialSecondary = true; - } else if (!matchedMCParticle.has_mothers()) { - isMaterialSecondary = true; + // Efficiency x Acceptance histograms + if (passedEvSelVtZ) { + histomc.fill(HIST("hDenomEffAcc"), particleType, mcParticle.pt(), mcParticle.y(), centrality, particleAnti, decayType); } } + } - if (isMaterialSecondary) { - if (pdg == particlePdgCodes.at(4)) { - histomc.fill(HIST("histRecoSecondaryMaterialHe3"), track.pt(), collision.centFT0C()); - } else if (pdg == -particlePdgCodes.at(4)) { - histomc.fill(HIST("histRecoSecondaryMaterialAntiHe3"), track.pt(), collision.centFT0C()); - } else if (pdg == particlePdgCodes.at(5)) { - histomc.fill(HIST("histRecoSecondaryMaterialHe4"), track.pt(), collision.centFT0C()); - } else if (pdg == -particlePdgCodes.at(5)) { - histomc.fill(HIST("histRecoSecondaryMaterialAntiHe4"), track.pt(), collision.centFT0C()); - } - - int type = 0; - if (std::abs(pdg) == particlePdgCodes.at(4)) - type = (pdg > 0) ? 1 : 2; - else if (std::abs(pdg) == particlePdgCodes.at(5)) - type = (pdg > 0) ? 3 : 4; - histomc.fill(HIST("histAllMaterialSecondariesReco"), track.pt(), getRapidity(track, he3), type); - - if (!cfgIncludeMaterialInEfficiency) { + // Process reconstructed collisions for this MC collision + if (passedEvSelVtZ) { + // Find the corresponding reconstructed collision + for (auto const& collision : collisions) { + if (collision.mcCollisionId() != static_cast(idx)) continue; - } - } - int decayType = 0; - bool isFromWeakDecay = false; - if (matchedMCParticle.isPhysicalPrimary()) { - decayType = 0; - if (matchedMCParticle.has_mothers()) { - for (const auto& motherparticle : matchedMCParticle.mothers_as()) { - if (std::find(hfMothCodes.begin(), hfMothCodes.end(), - std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { - isFromWeakDecay = true; - decayType = 1; - break; - } - } - } - } else if (matchedMCParticle.has_mothers()) { - isFromWeakDecay = true; - decayType = 1; - } else { - decayType = 2; - } + auto bc = collision.bc_as(); + initCCDB(bc); - if (!track.isPVContributor() && cfgUsePVcontributors) - continue; - if (!track.hasITS() && cfgITSrequire) - continue; - if (!track.hasTPC() && cfgTPCrequire) - continue; - if (!track.passedITSRefit() && cfgPassedITSRefit) - continue; - if (!track.passedTPCRefit() && cfgPassedTPCRefit) - continue; - if (std::abs(track.eta()) > cfgCutEta && cfgetaRequire) - continue; - if (!matchedMCParticle.isPhysicalPrimary() && !isFromWeakDecay && !isMaterialSecondary) - continue; + histomc.fill(HIST("histNevReco"), 0.5); + histomc.fill(HIST("histVtxZReco"), collision.posZ()); + histomc.fill(HIST("histCentFT0CReco"), collision.centFT0C()); + histomc.fill(HIST("histCentFT0MReco"), collision.centFT0M()); - for (size_t i = 0; i < primaryParticles.size(); i++) { - if (std::abs(pdg) != std::abs(particlePdgCodes.at(i))) - continue; - float ptMomn; - setTrackParCov(track, mTrackParCov); - mTrackParCov.setPID(track.pidForTracking()); - ptMomn = (i == he3 || i == he4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); - int sign = 0; - if (track.sign() > 0) { - sign = 1; - } - if (track.sign() < 0) { - sign = -1; - } - if (std::abs(getRapidity(track, i)) > cfgCutRapidity && cfgRapidityRequire) - continue; - if (track.tpcNClsFound() < cfgTrackPIDsettings->get(i, "minTPCnCls") && cfgTPCNClsfoundRequire) - continue; - if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || track.tpcNClsCrossedRows() < cfgtpcNClsFindable * track.tpcNClsFindable()) && cfgTPCNClsCrossedRowsRequire) - continue; - if (track.tpcChi2NCl() > cfgTrackPIDsettings->get(i, "maxTPCchi2") && cfgmaxTPCchi2Require) - continue; - if (track.tpcChi2NCl() < cfgTrackPIDsettings->get(i, "minTPCchi2") && cfgminTPCchi2Require) - continue; - if (track.itsNCls() < cfgTrackPIDsettings->get(i, "minITSnCls") && cfgminITSnClsRequire) - continue; - double cosheta = std::cosh(track.eta()); - if ((track.itsNCls() / cosheta) < cfgTrackPIDsettings->get(i, "minITSnClscos") && cfgminITSnClscosRequire) - continue; - if ((track.itsNClsInnerBarrel() < cfgTrackPIDsettings->get(i, "minReqClusterITSib")) && cfgminReqClusterITSibRequire) - continue; - if (track.itsChi2NCl() > cfgTrackPIDsettings->get(i, "maxITSchi2") && cfgmaxITSchi2Require) - continue; - if (getMeanItsClsSize(track) < cfgTrackPIDsettings->get(i, "minITSclsSize") && cfgminGetMeanItsClsSizeRequire) - continue; + auto tracksInColl = tracks.sliceBy(tracksPerCollision, collision.globalIndex()); - bool insideDCAxy = (std::abs(track.dcaXY()) <= (cfgTrackPIDsettings->get(i, "maxDcaXY") * (0.0105f + 0.0350f / std::pow(ptMomn, 1.1f)))); - if ((!(insideDCAxy) || std::abs(track.dcaZ()) > dcazSigma(ptMomn, cfgTrackPIDsettings->get(i, "maxDcaZ"))) && cfgDCAwithptRequire) - continue; + for (auto const& track : tracksInColl) { + if (!track.has_mcParticle()) + continue; // skip un-matched reco tracks - float tpcNsigma = getTPCnSigma(track, primaryParticles.at(i)); - if ((std::abs(tpcNsigma) > cfgTrackPIDsettings->get(i, "maxTPCnSigma")) && cfgmaxTPCnSigmaRequire) - continue; - float itsSigma = getITSnSigma(track, primaryParticles.at(i)); - if (itsSigma < cfgTrackPIDsettings2->get(i, "minITSnsigma") && cfgTrackPIDsettings2->get(i, "useITSnsigma") < 1) - continue; - if (itsSigma > cfgTrackPIDsettings2->get(i, "maxITSnsigma") && cfgTrackPIDsettings2->get(i, "useITSnsigma") < 1) - continue; - float ptReco; - mTrackParCov.setPID(track.pidForTracking()); - ptReco = (std::abs(pdg) == particlePdgCodes.at(4) || std::abs(pdg) == particlePdgCodes.at(5)) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); + auto const& matchedMCParticle = track.mcParticle_as(); - int particleAnti = (pdg > 0) ? 0 : 1; + // Only process particles from this MC collision + if (matchedMCParticle.mcCollisionId() != mcCollision.globalIndex()) + continue; - if (i == he3 || i == he4) { - histomc.fill(HIST("hNumerEffAcc"), i, matchedMCParticle.pt(), matchedMCParticle.y(), collision.centFT0C(), particleAnti, decayType); - } - if (isFromWeakDecay) { - if (std::abs(pdg) == particlePdgCodes.at(4)) { - histomc.fill(HIST("histRecoWeakDecayPtHe3"), ptReco, collision.centFT0C()); - } else if (std::abs(pdg) == particlePdgCodes.at(5)) { - histomc.fill(HIST("histRecoWeakDecayPtHe4"), ptReco, collision.centFT0C()); - } - } + int pdg = matchedMCParticle.pdgCode(); + bool isHe3 = (std::abs(pdg) == particlePdgCodes.at(4)); + bool isHe4 = (std::abs(pdg) == particlePdgCodes.at(5)); - histos.fill(HIST("Tpcsignal"), getRigidity(track) * track.sign(), track.tpcSignal()); + if (!isHe3 && !isHe4) + continue; - if (pdg == -particlePdgCodes.at(5) && cfgmccorrectionhe4Require) { - ptReco = ptReco + 0.00765 + 0.503791 * std::exp(-1.10517 * ptReco); - } + int decayType = 0; + bool isFromWeakDecay = false; - if (pdg == -particlePdgCodes.at(4) && cfgmccorrectionhe4Require) { - int pidGuess = track.pidForTracking(); - int antitriton = 6; - if (pidGuess == antitriton) { - ptReco = ptReco - 0.464215 + 0.195771 * ptReco - 0.0183111 * ptReco * ptReco; + if (matchedMCParticle.isPhysicalPrimary()) { + decayType = 0; + if (matchedMCParticle.has_mothers()) { + for (const auto& motherparticle : matchedMCParticle.mothers_as()) { + if (std::find(hfMothCodes.begin(), hfMothCodes.end(), + std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { + isFromWeakDecay = true; + decayType = 1; + break; + } + } + } + } else if (matchedMCParticle.has_mothers()) { + isFromWeakDecay = true; + decayType = 1; + } else { + decayType = 2; } - } - float ptGen = matchedMCParticle.pt(); - float deltaPt = ptReco - ptGen; + if (!track.isPVContributor() && cfgUsePVcontributors) + continue; + if (!track.hasITS() && cfgITSrequire) + continue; + if (!track.hasTPC() && cfgTPCrequire) + continue; + if (!track.passedITSRefit() && cfgPassedITSRefit) + continue; + if (!track.passedTPCRefit() && cfgPassedTPCRefit) + continue; + if (std::abs(track.eta()) > cfgCutEta && cfgetaRequire) + continue; + if (!matchedMCParticle.isPhysicalPrimary() && isFromWeakDecay) + continue; + + for (size_t i = 0; i < primaryParticles.size(); i++) { + if (std::abs(pdg) != std::abs(particlePdgCodes.at(i))) + continue; + + if (std::abs(getRapidity(track, i)) > cfgCutRapidity && cfgRapidityRequire) + continue; + + float ptReco; + setTrackParCov(track, mTrackParCov); + mTrackParCov.setPID(track.pidForTracking()); + + ptReco = (std::abs(pdg) == particlePdgCodes.at(4) || std::abs(pdg) == particlePdgCodes.at(5)) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); + + int particleAnti = (pdg > 0) ? 0 : 1; + + if (pdg == -particlePdgCodes.at(5) && cfgmccorrectionhe4Require) { + ptReco = ptReco + 0.00765 + 0.503791 * std::exp(-1.10517 * ptReco); + } - if (pdg == -particlePdgCodes.at(4)) { - histomc.fill(HIST("histDeltaPtVsPtGenanti"), ptReco, deltaPt); - histomc.fill(HIST("histPIDtrackanti"), ptReco, track.pidForTracking()); - } - if (pdg == particlePdgCodes.at(4)) { - histomc.fill(HIST("histDeltaPtVsPtGen"), ptReco, deltaPt); - histomc.fill(HIST("histPIDtrack"), ptReco, track.pidForTracking()); - } - if (pdg == -particlePdgCodes.at(5)) { - histomc.fill(HIST("histDeltaPtVsPtGenHe4anti"), ptReco, deltaPt); - histomc.fill(HIST("histPIDtrackantihe4"), ptReco, track.pidForTracking()); - } - if (pdg == particlePdgCodes.at(5)) { - histomc.fill(HIST("histDeltaPtVsPtGenHe4"), ptReco, deltaPt); - histomc.fill(HIST("histPIDtrackhe4"), ptReco, track.pidForTracking()); - } + if (pdg == -particlePdgCodes.at(4) && cfgmccorrectionhe4Require) { + int pidGuess = track.pidForTracking(); + int antitriton = 6; + if (pidGuess == antitriton) { + ptReco = ptReco - 0.464215 + 0.195771 * ptReco - 0.0183111 * ptReco * ptReco; + } + } - // For TOF matching efficiency - float ptTOF = -1.0; // Default: no TOF - if (track.hasTOF()) { - ptTOF = ptReco; - } + if (track.tpcNClsFound() < cfgTrackPIDsettings->get(i, "minTPCnCls") && cfgTPCNClsfoundRequire) + continue; + if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || track.tpcNClsCrossedRows() < cfgtpcNClsFindable * track.tpcNClsFindable()) && cfgTPCNClsCrossedRowsRequire) + continue; + if (track.tpcChi2NCl() > cfgTrackPIDsettings->get(i, "maxTPCchi2") && cfgmaxTPCchi2Require) + continue; + if (track.tpcChi2NCl() < cfgTrackPIDsettings->get(i, "minTPCchi2") && cfgminTPCchi2Require) + continue; + if (track.itsNCls() < cfgTrackPIDsettings->get(i, "minITSnCls") && cfgminITSnClsRequire) + continue; + double cosheta = std::cosh(track.eta()); + if ((track.itsNCls() / cosheta) < cfgTrackPIDsettings->get(i, "minITSnClscos") && cfgminITSnClscosRequire) + continue; + if ((track.itsNClsInnerBarrel() < cfgTrackPIDsettings->get(i, "minReqClusterITSib")) && cfgminReqClusterITSibRequire) + continue; + if (track.itsChi2NCl() > cfgTrackPIDsettings->get(i, "maxITSchi2") && cfgmaxITSchi2Require) + continue; + if (getMeanItsClsSize(track) < cfgTrackPIDsettings->get(i, "minITSclsSize") && cfgminGetMeanItsClsSizeRequire) + continue; + + // DCA XY cut + bool insideDCAxy = cfgdcaxynopt ? (std::abs(track.dcaXY()) <= cfgTrackPIDsettings->get(i, "maxDcaXY")) : (std::abs(track.dcaXY()) <= (cfgTrackPIDsettings->get(i, "maxDcaXY") * (0.0105f + 0.0350f / std::pow(ptReco, 1.1f)))); + + // DCA Z cut + bool insideDCAz = cfgdcaznopt ? (std::abs(track.dcaZ()) <= cfgTrackPIDsettings->get(i, "maxDcaZ")) : (std::abs(track.dcaZ()) <= dcazSigma(ptReco, cfgTrackPIDsettings->get(i, "maxDcaZ"))); + + if ((!insideDCAxy || !insideDCAz)) { + continue; + } - if (cfgTrackPIDsettings2->get(i, "fillsparsh") == 1) { - histomc.fill(HIST("hSpectramc"), i, sign, collision.centFT0C(), - ptReco, ptTOF); - } + if (i == he3 || i == he4) { + histomc.fill(HIST("hNumerEffAcc"), i, ptReco, matchedMCParticle.y(), collision.centFT0C(), particleAnti, decayType); + } - if (cfgRequirebetaplot) { - histos.fill(HIST("Tofsignal"), getRigidity(track) * track.sign(), o2::pid::tof::Beta::GetBeta(track)); + histos.fill(HIST("Tpcsignal"), getRigidity(track) * track.sign(), track.tpcSignal()); + + // Fill the requested histograms + float ptGen = matchedMCParticle.pt(); + float deltaPt = ptReco - ptGen; + + if (pdg == -particlePdgCodes.at(4)) { + histomc.fill(HIST("histDeltaPtVsPtGenanti"), ptReco, deltaPt); + histomc.fill(HIST("histPIDtrackanti"), ptReco, track.pidForTracking()); + } + if (pdg == particlePdgCodes.at(4)) { + histomc.fill(HIST("histDeltaPtVsPtGen"), ptReco, deltaPt); + histomc.fill(HIST("histPIDtrack"), ptReco, track.pidForTracking()); + } + if (pdg == -particlePdgCodes.at(5)) { + histomc.fill(HIST("histDeltaPtVsPtGenHe4anti"), ptReco, deltaPt); + histomc.fill(HIST("histPIDtrackantihe4"), ptReco, track.pidForTracking()); + } + if (pdg == particlePdgCodes.at(5)) { + histomc.fill(HIST("histDeltaPtVsPtGenHe4"), ptReco, deltaPt); + histomc.fill(HIST("histPIDtrackhe4"), ptReco, track.pidForTracking()); + } + } } - } // loop over primaryParticles types - histos.fill(HIST("histeta"), track.eta()); - } // TrackAssoc loop - } // Collision loop + break; // Found the matching collision, break out of collision loop + } + } + } } - PROCESS_SWITCH(NucleitpcPbPb, processMC, "MC reco+gen analysis with efficiency corrections", false); //=-=-=-==-=-=-==-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= void initCCDB(aod::BCsWithTimestamps::iterator const& bc) From 2e2a07e59e0c6b4009a1a244f6d5ecb5860a9202 Mon Sep 17 00:00:00 2001 From: Jaideep Tanwar <141036812+jtanwar2212@users.noreply.github.com> Date: Wed, 24 Sep 2025 14:25:50 +0530 Subject: [PATCH 8/8] Update nucleitpcpbpb.cxx --- PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx index cadac3bf9e1..d1634c4f294 100644 --- a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx +++ b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx @@ -742,6 +742,16 @@ struct NucleitpcPbPb { histomc.fill(HIST("hNumerEffAcc"), i, ptReco, matchedMCParticle.y(), collision.centFT0C(), particleAnti, decayType); } + float ptTOF = -1.0; // Default: no TOF + if (track.hasTOF()) { + ptTOF = ptReco; + } + + if (cfgTrackPIDsettings2->get(i, "fillsparsh") == 1) { + histomc.fill(HIST("hSpectramc"), i, particleAnti, collision.centFT0C(), + ptReco, ptTOF); + } + histos.fill(HIST("Tpcsignal"), getRigidity(track) * track.sign(), track.tpcSignal()); // Fill the requested histograms