diff --git a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx index 237c91ca787..f1ee505a4e5 100644 --- a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx +++ b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx @@ -128,7 +128,6 @@ struct NucleitpcPbPb { HistogramRegistry histomc{"histomc", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; // Event Selection Configurables - Configurable removeITSROFrameBorder{"removeITSROFrameBorder", false, "Remove ITS RO frame border"}; Configurable removeNoSameBunchPileup{"removeNoSameBunchPileup", false, "Remove no same bunch pileup"}; Configurable requireIsGoodZvtxFT0vsPV{"requireIsGoodZvtxFT0vsPV", false, "Require is good Zvtx FT0 vs PV"}; Configurable requireIsVertexITSTPC{"requireIsVertexITSTPC", false, "Require is vertex ITS TPC"}; @@ -183,7 +182,6 @@ struct NucleitpcPbPb { Configurable cfgRigidityCorrection{"cfgRigidityCorrection", false, "apply rigidity correction"}; Configurable cfgRequirebetaplot{"cfgRequirebetaplot", true, "Require beta plot"}; Configurable cfgmass2{"cfgmass2", true, "Fill mass square difference"}; - Configurable cfgFillhspectra{"cfgFillhspectra", true, "fill data sparsh"}; Configurable cfgFillmass{"cfgFillmass", false, "Fill mass histograms"}; Configurable cfgFillmassnsigma{"cfgFillmassnsigma", true, "Fill mass vs nsigma histograms"}; @@ -207,12 +205,9 @@ struct NucleitpcPbPb { ConfigurableAxis axismass{"axismass", {1200, 0, 12}, "mass"}; ConfigurableAxis axismassnsigma{"axismassnsigma", {1200, 0, 12}, "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.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 - ConfigurableAxis axisTPCsig{"axisTPCsig", {1000, 0, 2000}, "TPC signal (a.u.)"}; // CCDB Service ccdb; @@ -230,6 +225,8 @@ struct NucleitpcPbPb { int mRunNumber, occupancy; float dBz; TRandom3 rand; + + float de = 2; float triton = 3; float he3 = 4; float he4 = 5; @@ -254,99 +251,80 @@ 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, axisOccupancy}); - histos.add("DCAxy_vs_pT_He3_data", "DCA_{xy} vs p_{T} for He3 (Data);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH2F, {ptAxis, axisDCA}}); - histos.add("DCAxy_vs_pT_antiHe3_data", "DCA_{xy} vs p_{T} for anti-He3 (Data);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH2F, {ptAxis, axisDCA}}); - - histos.add("TPCsig_triton", "TPC signal for triton;p/z (GeV/c);TPC signal", - {HistType::kTH2F, {axisRigidity, axisTPCsig}}); - histos.add("TPCsig_he3", "TPC signal for he3;p/z (GeV/c);TPC signal", - {HistType::kTH2F, {axisRigidity, axisTPCsig}}); - histos.add("TPCsig_he4", "TPC signal for he4;p/z (GeV/c);TPC signal", - {HistType::kTH2F, {axisRigidity, axisTPCsig}}); + histos.add("hSpectra", " ", HistType::kTHnSparseF, {ptAxis, nsigmaAxis, {5, -2.5, 2.5}, axisCent}); + histos.add("DCAxy_vs_pT_data", "DCA_{xy} vs p_{T} for He3 (Data);p_{T} (GeV/c);DCA_{xy} (cm)", + {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); + histos.add("DCAxy_vs_pT_anti_data", "DCA_{xy} vs p_{T} for anti-He3 (Data);p_{T} (GeV/c);DCA_{xy} (cm)", + {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); + + hmass.resize(2 * nParticles + 2); + hmassnsigma.resize(2 * nParticles + 2); + + 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}; 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()), ";p_T{TPC} (GeV/#it{c}); mass^{2}/z^{2}; Centrality(%)", HistType::kTH3F, {ptAxis, axismassnsigma, axisCent}); + hmassnsigma[2 * i + 1] = histos.add(Form("histmass_nsigmaanti/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}/z^{2}; centrality(%)", HistType::kTH3F, {ptAxis, axismassnsigma, axisCent}); + } + } + + histos.add("histeta", "histeta", kTH1F, {axiseta}); + + histos.add("histEvents", "histEvents", kTH2F, {axisCent, axisOccupancy}); } - histos.add("histeta", "histeta", kTH1F, {axiseta}); - histos.add("histEvents", "histEvents", kTH2F, {axisCent, axisOccupancy}); + histos.add("dcaZ", "dcaZ", kTH2F, {ptAxis, axisDCA}); histos.add("dcaXY", "dcaXY", kTH2F, {ptAxis, axisDCA}); histos.add("Tofsignal", "Tofsignal", kTH2F, {axisRigidity, {4000, 0.2, 1.2, "#beta"}}); histos.add("Tpcsignal", "Tpcsignal", kTH2F, {axisRigidity, axisdEdx}); - hmass.resize(2 * nParticles + 2); - hmassnsigma.resize(2 * nParticles + 2); - - 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}; 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()), ";p_T{TPC} (GeV/#it{c}); mass^{2}/z^{2}; Centrality(%)", HistType::kTH3F, {ptAxis, axismassnsigma, axisCent}); - hmassnsigma[2 * i + 1] = histos.add(Form("histmass_nsigmaanti/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}/z^{2}; centrality(%)", HistType::kTH3F, {ptAxis, axismassnsigma, axisCent}); - } - } - if (doprocessMC) { - histomc.add("hSpectramc", " ", HistType::kTHnSparseF, {speciesBitAxis, {5, -2.5, 2.5}, axisCent, ptAxis, ptAxis}); + histomc.add("hSpectramc", " ", HistType::kTHnSparseF, {{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}}); + {HistType::kTHnSparseF, {ptAxis, axisCent, particleAntiAxis, decayTypeAxis}}); histomc.add("hNumerEffAcc", "Numerator for Efficiency x Acceptance", - {HistType::kTHnSparseF, {speciesBitAxis, ptAxis, axisrapidity, axisCent, particleAntiAxis, decayTypeAxis}}); + {HistType::kTHnSparseF, {ptAxis, axisCent, particleAntiAxis, decayTypeAxis}}); // 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("hSignalLossDenom", "Signal Loss Denominator", kTH2F, {axisCent, ptAxis}); + histomc.add("hSignalLossNumer", "Signal Loss Numerator", kTH2F, {axisCent, ptAxis}); - 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}); + histomc.add("haSignalLossDenom", "anti particle Signal Loss Denominator", kTH2F, {axisCent, ptAxis}); + histomc.add("haSignalLossNumer", "antiparticle Signal Loss Numerator", kTH2F, {axisCent, ptAxis}); // 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("histNevReco", "histNevReco", kTH1F, {axisNev}); histomc.add("histVtxZReco", "histVtxZReco", kTH1F, {axisVtxZ}); - histomc.add("histCentFT0CReco", "histCentFT0CReco", kTH1F, {axisCent}); - histomc.add("histCentFT0MReco", "histCentFT0MReco", kTH1F, {axisCent}); - - histomc.add("histdetapttriton", " delta pt vs pt rec for trition detected", HistType::kTH2F, {{1000, 0, 10}, {1000, -0.5, 0.5, "p_{T}(reco) - p_{T}(gen);p_{T}(reco)"}}); - histomc.add("histdetapttritonanti", " delta pt vs pt rec for trition detected", HistType::kTH2F, {{1000, 0, 10}, {1000, -0.5, 0.5, "p_{T}(reco) - p_{T}(gen);p_{T}(reco)"}}); - 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("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("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("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)"}}); } if (doprocessDCA) { - histomc.add("hSpectraDCA", " ", HistType::kTHnSparseF, {speciesBitAxis, {5, -2.5, 2.5}, axisCent, ptAxis, ptAxis, decayTypeAxis, axisDCA}); // NEW: Add DCAxy vs pT histograms for MC - histomc.add("DCAxy_vs_pT_He3_transport", "DCA_{xy} vs p_{T} for He3 (Transport);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH2F, {ptAxis, axisDCA}}); - histomc.add("DCAxy_vs_pT_He3_weakdecay", "DCA_{xy} vs p_{T} for He3 (Weak Decay);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH2F, {ptAxis, axisDCA}}); - histomc.add("DCAxy_vs_pT_He3_total", "DCA_{xy} vs p_{T} for He3 (Total);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH2F, {ptAxis, axisDCA}}); - histomc.add("DCAxy_vs_pT_antiHe3_total", "DCA_{xy} vs p_{T} for anti-He3 (Total);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH2F, {ptAxis, axisDCA}}); + histomc.add("DCAxy_vs_pT_transport", "DCA_{xy} vs p_{T} (Transport);p_{T} (GeV/c);DCA_{xy} (cm)", + {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); + histomc.add("DCAxy_vs_pT_weakdecay", "DCA_{xy} vs p_{T} (Weak Decay);p_{T} (GeV/c);DCA_{xy} (cm)", + {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); + histomc.add("DCAxy_vs_pT_total", "DCA_{xy} vs p_{T} (Total);p_{T} (GeV/c);DCA_{xy} (cm)", + {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); + histomc.add("DCAxy_vs_pT_anti_total", "DCA_{xy} vs p_{T} for anti (Total);p_{T} (GeV/c);DCA_{xy} (cm)", + {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); } } //---------------------------------------------------------------------------------------------------------------- @@ -362,39 +340,34 @@ struct NucleitpcPbPb { if (!collPassedEvSel) 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); + histos.fill(HIST("histNev"), 2.5); if (requireIsGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) continue; - histos.fill(HIST("histNev"), 4.5); + histos.fill(HIST("histNev"), 3.5); if (requireIsVertexITSTPC && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) continue; - histos.fill(HIST("histNev"), 5.5); + histos.fill(HIST("histNev"), 4.5); if (removeNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) continue; - histos.fill(HIST("histNev"), 6.5); + histos.fill(HIST("histNev"), 5.5); histos.fill(HIST("histEvents"), collision.centFT0C(), occupancy); histos.fill(HIST("histVtxZ"), collision.posZ()); histos.fill(HIST("histCentFT0C"), collision.centFT0C()); histos.fill(HIST("histCentFT0M"), collision.centFT0M()); if (collision.centFT0C() > centcut) continue; - histos.fill(HIST("histNev"), 7.5); + histos.fill(HIST("histNev"), 6.5); histos.fill(HIST("histCentFTOC_cut"), collision.centFT0C()); // 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.hasTPC()) continue; if (!track.passedITSRefit() && cfgPassedITSRefit) @@ -403,13 +376,13 @@ struct NucleitpcPbPb { continue; if (std::abs(track.eta()) > cfgCutEta && cfgetaRequire) continue; - for (size_t i = 0; i < primaryParticles.size(); i++) { + if (cfgTrackPIDsettings2->get(i, "fillsparsh") != 1) + continue; float ptMomn; setTrackParCov(track, mTrackParCov); mTrackParCov.setPID(track.pidForTracking()); ptMomn = (i == he3 || i == he4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); - float tpcNsigma = getTPCnSigma(track, primaryParticles.at(i)); if ((std::abs(tpcNsigma) > cfgTrackPIDsettings->get(i, "maxTPCnSigma")) && cfgmaxTPCnSigmaRequire) continue; @@ -422,23 +395,19 @@ struct NucleitpcPbPb { } else if (i == he4) { param = (track.sign() > 0) ? 2 : 3; } - if (param >= 0) { a = cfgktrackcorrection->get(param, "a"); b = cfgktrackcorrection->get(param, "b"); c = cfgktrackcorrection->get(param, "c"); } - if (i == he4 && ptMomn < pTcorrectHe4) { ptMomn = ptMomn + a + b * std::exp(c * ptMomn); } - if (i == he3 && ptMomn < pTcorrectHe3) { ptMomn = ptMomn - a + b * ptMomn - c * ptMomn * ptMomn; } } 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")) @@ -453,7 +422,6 @@ 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; @@ -463,9 +431,7 @@ struct NucleitpcPbPb { continue; if (getMeanItsClsSize(track) < cfgTrackPIDsettings->get(i, "minITSclsSize") && cfgminGetMeanItsClsSizeRequire) continue; - bool insideDCAxy = false; - if (cfgUseNewDCAxyCut) { double sigmaFactor = cfgTrackPIDsettings->get(i, "maxDcaXY"); double sigma_base = (0.0118 * std::exp(-0.6889 * ptMomn) + 0.0017); @@ -514,51 +480,37 @@ struct NucleitpcPbPb { if (i == he4) { tpcNsigmaDe = track.tpcNSigmaDe(); } - if (i == triton || i == he3 || i == he4) { - float rigidity = getRigidity(track); - if (i == triton) { - histos.fill(HIST("TPCsig_triton"), rigidity * track.sign(), track.tpcSignal()); - } else if (i == he3) { - histos.fill(HIST("TPCsig_he3"), rigidity * track.sign(), track.tpcSignal()); - } else if (i == he4) { - histos.fill(HIST("TPCsig_he4"), rigidity * track.sign(), track.tpcSignal()); - } - } - if (i == he3) { - if (track.sign() > 0) { - histos.fill(HIST("DCAxy_vs_pT_He3_data"), ptMomn, track.dcaXY()); - } else if (track.sign() < 0) { - histos.fill(HIST("DCAxy_vs_pT_antiHe3_data"), ptMomn, track.dcaXY()); - } + if (track.sign() > 0) { + histos.fill(HIST("DCAxy_vs_pT_data"), ptMomn, track.dcaXY(), collision.centFT0C()); + } else if (track.sign() < 0) { + histos.fill(HIST("DCAxy_vs_pT_anti_data"), ptMomn, track.dcaXY(), collision.centFT0C()); } - if (cfgFillhspectra && cfgTrackPIDsettings2->get(i, "fillsparsh") == 1) { - if (i != he4) { - histos.fill(HIST("hSpectra"), i, ptMomn, tpcNsigma, sign, collision.centFT0C(), track.dcaZ(), track.dcaXY(), collision.trackOccupancyInTimeRange()); + if (i != he4) { + histos.fill(HIST("hSpectra"), ptMomn, tpcNsigma, sign, collision.centFT0C()); + } else { + if (!track.hasTOF()) { + if (std::abs(tpcNsigmaDe) > deuteronsigmarejection) { + histos.fill(HIST("hSpectra"), ptMomn, tpcNsigma, sign, collision.centFT0C()); + } } else { - if (!track.hasTOF()) { - if (std::abs(tpcNsigmaDe) > deuteronsigmarejection) { - histos.fill(HIST("hSpectra"), i, ptMomn, tpcNsigma, sign, collision.centFT0C(), track.dcaZ(), track.dcaXY(), collision.trackOccupancyInTimeRange()); - } - } 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 (cfghe3massrejreq && (massTOF * massTOF > cfgminmassrejection && massTOF * massTOF < cfgmaxmassrejection)) { - continue; // Skip if mass cut fails - } - if (std::abs(tpcNsigmaDe) > deuteronsigmarejection) { - histos.fill(HIST("hSpectra"), i, ptMomn, tpcNsigma, sign, collision.centFT0C(), track.dcaZ(), track.dcaXY(), collision.trackOccupancyInTimeRange()); - } + // 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 (cfghe3massrejreq && (massTOF * massTOF > cfgminmassrejection && massTOF * massTOF < cfgmaxmassrejection)) { + continue; // Skip if mass cut fails + } + if (std::abs(tpcNsigmaDe) > deuteronsigmarejection) { + histos.fill(HIST("hSpectra"), ptMomn, tpcNsigma, sign, collision.centFT0C()); } } } @@ -631,9 +583,6 @@ struct NucleitpcPbPb { 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)) @@ -656,7 +605,7 @@ struct NucleitpcPbPb { // 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 + // Event loss denominator (always filled regardless of fillsparsh) histomc.fill(HIST("hEventLossDenom"), mcCollInfos[i].centrality); // Event loss numerator (if passed selection) @@ -673,29 +622,31 @@ struct NucleitpcPbPb { continue; } - // Signal loss denominator - if (mcParticle.pdgCode() == particlePdgCodes.at(4)) { // He3 - histomc.fill(HIST("hHe3SignalLossDenom"), mcCollInfos[i].centrality); + // Check if this particle type has fillsparsh enabled + int particleType = -1; + for (size_t j = 0; j < primaryParticles.size(); j++) { + if (std::abs(mcParticle.pdgCode()) == std::abs(particlePdgCodes.at(j))) { + particleType = j; + break; + } + } - } 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); + // Skip if fillsparsh is not enabled for this particle + if (particleType < 0 || cfgTrackPIDsettings2->get(particleType, "fillsparsh") != 1) + continue; - } else if (mcParticle.pdgCode() == -particlePdgCodes.at(5)) { // He4 - histomc.fill(HIST("haHe4SignalLossDenom"), mcCollInfos[i].centrality); + // Signal loss denominator + if (mcParticle.pdgCode() == particlePdgCodes.at(particleType)) { // particle + histomc.fill(HIST("hSignalLossDenom"), mcCollInfos[i].centrality, mcParticle.pt()); + } else if (mcParticle.pdgCode() == -particlePdgCodes.at(particleType)) { // anti-particle + histomc.fill(HIST("haSignalLossDenom"), mcCollInfos[i].centrality, mcParticle.pt()); } - // 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); + if (mcParticle.pdgCode() == particlePdgCodes.at(particleType)) { // particle + histomc.fill(HIST("hSignalLossNumer"), mcCollInfos[i].centrality, mcParticle.pt()); + } else if (mcParticle.pdgCode() == -particlePdgCodes.at(particleType)) { // anti-particle + histomc.fill(HIST("haSignalLossNumer"), mcCollInfos[i].centrality, mcParticle.pt()); } } } @@ -709,7 +660,6 @@ struct NucleitpcPbPb { continue; float centrality = mcCollInfos[idx].centrality; - // bool passedEvSel = mcCollInfos[idx].passedEvSel; bool passedEvSelVtZ = mcCollInfos[idx].passedEvSelVtZ; // Process generated particles for efficiency denominators @@ -720,10 +670,18 @@ struct NucleitpcPbPb { continue; int pdgCode = mcParticle.pdgCode(); - bool isHe3 = (std::abs(pdgCode) == particlePdgCodes.at(4)); - bool isHe4 = (std::abs(pdgCode) == particlePdgCodes.at(5)); - if (!isHe3 && !isHe4) + // Check which particle type this is + int particleType = -1; + for (size_t j = 0; j < primaryParticles.size(); j++) { + if (std::abs(pdgCode) == std::abs(particlePdgCodes.at(j))) { + particleType = j; + break; + } + } + + // Only process if this particle type has fillsparsh enabled + if (particleType < 0 || cfgTrackPIDsettings2->get(particleType, "fillsparsh") != 1) continue; if (std::abs(mcParticle.eta()) > cfgCutEta && cfgetaRequireMC) @@ -751,24 +709,10 @@ struct NucleitpcPbPb { decayType = 2; continue; } - bool isFromWeakDecay = (decayType == 1); - if (!mcParticle.isPhysicalPrimary() && !isFromWeakDecay) - // if (!mcParticle.isPhysicalPrimary()) - 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) { - - // Efficiency x Acceptance histograms - if (passedEvSelVtZ) { - histomc.fill(HIST("hDenomEffAcc"), particleType, mcParticle.pt(), mcParticle.y(), centrality, particleAnti, decayType); - } + // Efficiency x Acceptance histograms - Denominator + if (passedEvSelVtZ) { + histomc.fill(HIST("hDenomEffAcc"), mcParticle.pt(), centrality, particleAnti, decayType); } } @@ -782,10 +726,7 @@ struct NucleitpcPbPb { auto bc = collision.bc_as(); initCCDB(bc); - histomc.fill(HIST("histNevReco"), 0.5); histomc.fill(HIST("histVtxZReco"), collision.posZ()); - histomc.fill(HIST("histCentFT0CReco"), collision.centFT0C()); - histomc.fill(HIST("histCentFT0MReco"), collision.centFT0M()); auto tracksInColl = tracks.sliceBy(tracksPerCollision, collision.globalIndex()); @@ -800,11 +741,6 @@ struct NucleitpcPbPb { continue; int pdg = matchedMCParticle.pdgCode(); - bool isHe3 = (std::abs(pdg) == particlePdgCodes.at(4)); - bool isHe4 = (std::abs(pdg) == particlePdgCodes.at(5)); - - if (!isHe3 && !isHe4) - continue; int decayType = 0; @@ -836,13 +772,15 @@ struct NucleitpcPbPb { continue; if (std::abs(track.eta()) > cfgCutEta && cfgetaRequire) continue; - if (!matchedMCParticle.isPhysicalPrimary()) - continue; for (size_t i = 0; i < primaryParticles.size(); i++) { if (std::abs(pdg) != std::abs(particlePdgCodes.at(i))) continue; + // CRITICAL: Only process if fillsparsh is enabled for this particle + if (cfgTrackPIDsettings2->get(i, "fillsparsh") != 1) + continue; + float ptReco; setTrackParCov(track, mTrackParCov); mTrackParCov.setPID(track.pidForTracking()); @@ -855,9 +793,9 @@ struct NucleitpcPbPb { int param = -1; if (i == he3) { - param = (-particlePdgCodes.at(4) > 0) ? 0 : 1; + param = (pdg > 0) ? 0 : 1; } else if (i == he4) { - param = (-particlePdgCodes.at(4) > 0) ? 2 : 3; + param = (pdg > 0) ? 2 : 3; } if (param >= 0) { @@ -907,31 +845,22 @@ struct NucleitpcPbPb { double sigmaFactor = cfgTrackPIDsettings->get(i, "maxDcaXY"); double sigma_base = (0.0118 * std::exp(-0.6889 * ptReco) + 0.0017); double sigma_new = sigmaFactor * sigma_base; - insideDCAxy = (std::abs(track.dcaXY()) <= sigma_new); } else { - 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)))); + 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)))); } bool insideDCAz = false; if (cfgUseNewDCAzCut) { - double sigmaFactorZ = cfgTrackPIDsettings->get(i, "maxDcaZ"); - double p0 = 0.1014; double p1 = 1.7512; double p2 = 0.0024; - double sigma_base_z = p0 * std::exp(-p1 * ptReco) + p2; - double sigma_new_z = sigmaFactorZ * sigma_base_z; - insideDCAz = (std::abs(track.dcaZ()) <= sigma_new_z); } else { insideDCAz = (std::abs(track.dcaZ()) <= cfgTrackPIDsettings->get(i, "maxDcaZ")); @@ -945,56 +874,32 @@ struct NucleitpcPbPb { if ((std::abs(tpcNsigma) > cfgTrackPIDsettings->get(i, "maxTPCnSigma")) && cfgmaxTPCnSigmaRequire) continue; - if (i == he3 || i == he4) { - histomc.fill(HIST("hNumerEffAcc"), i, ptReco, getRapidity(track, i), collision.centFT0C(), particleAnti, decayType); - } + // Efficiency x Acceptance - Numerator + histomc.fill(HIST("hNumerEffAcc"), ptReco, collision.centFT0C(), particleAnti, decayType); - float ptTOF = -1.0; // Default: no TOF + // Sparse histogram + float ptTOF = -1.0; if (track.hasTOF()) { ptTOF = ptReco; } + histomc.fill(HIST("hSpectramc"), particleAnti, collision.centFT0C(), ptReco, ptTOF); - if (cfgTrackPIDsettings2->get(i, "fillsparsh") == 1) { - histomc.fill(HIST("hSpectramc"), i, particleAnti, collision.centFT0C(), - ptReco, ptTOF); - } - + // Basic track histograms histos.fill(HIST("dcaXY"), ptReco, track.dcaXY()); histos.fill(HIST("dcaZ"), ptReco, track.dcaZ()); - histos.fill(HIST("Tpcsignal"), getRigidity(track) * track.sign(), track.tpcSignal()); - // Fill the requested histograms + // Delta Pt histograms float ptGen = matchedMCParticle.pt(); float deltaPt = ptReco - ptGen; - if (pdg == -particlePdgCodes.at(4)) { + if (pdg == -particlePdgCodes.at(i)) { // Anti-particle histomc.fill(HIST("histDeltaPtVsPtGenanti"), ptReco, deltaPt); histomc.fill(HIST("histPIDtrackanti"), ptReco, track.pidForTracking()); - - int pidGuess = track.pidForTracking(); - int antitriton = 6; - if (pidGuess == antitriton) { - histomc.fill(HIST("histdetapttritonanti"), ptReco, deltaPt); - } } - if (pdg == particlePdgCodes.at(4)) { + if (pdg == particlePdgCodes.at(i)) { // Particle histomc.fill(HIST("histDeltaPtVsPtGen"), ptReco, deltaPt); histomc.fill(HIST("histPIDtrack"), ptReco, track.pidForTracking()); - - int pidGuess = track.pidForTracking(); - int antitriton = 6; - if (pidGuess == antitriton) { - histomc.fill(HIST("histdetapttriton"), ptReco, deltaPt); - } - } - 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()); } } } @@ -1005,9 +910,6 @@ struct NucleitpcPbPb { } PROCESS_SWITCH(NucleitpcPbPb, processMC, "MC reco+gen analysis with efficiency corrections", false); //=-=-=-==-=-=-==-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= - //---------------------------------------------------------------------------------------------------------------- - // MC particles - DCA secondary fraction - //-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= void processDCA(CollisionsFullMC const& collisions, aod::McCollisions const& mcCollisions, aod::McParticles const& particlesMC, @@ -1034,9 +936,6 @@ struct NucleitpcPbPb { 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)) @@ -1086,11 +985,6 @@ struct NucleitpcPbPb { continue; int pdg = matchedMCParticle.pdgCode(); - bool isHe3 = (std::abs(pdg) == particlePdgCodes.at(4)); - bool isHe4 = (std::abs(pdg) == particlePdgCodes.at(5)); - - if (!isHe3 && !isHe4) - continue; if (!track.isPVContributor() && cfgUsePVcontributors) continue; @@ -1108,14 +1002,15 @@ struct NucleitpcPbPb { if (std::abs(pdg) != std::abs(particlePdgCodes.at(i))) continue; + if (cfgTrackPIDsettings2->get(i, "fillsparsh") != 1) + 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; - double a = 0, b = 0, c = 0; int param = -1; @@ -1199,67 +1094,23 @@ struct NucleitpcPbPb { decayType = 2; } - if (i == he3) { - if (pdg == particlePdgCodes.at(4)) { // He3 - histomc.fill(HIST("DCAxy_vs_pT_He3_total"), ptReco, track.dcaXY()); - if (decayType == 2) { // Transport/Material - histomc.fill(HIST("DCAxy_vs_pT_He3_transport"), ptReco, track.dcaXY()); - } else if (decayType == 1) { // Weak decay (including HF) - histomc.fill(HIST("DCAxy_vs_pT_He3_weakdecay"), ptReco, track.dcaXY()); - } - } else if (pdg == -particlePdgCodes.at(4)) { // anti-He3 - histomc.fill(HIST("DCAxy_vs_pT_antiHe3_total"), ptReco, track.dcaXY()); - } - } - - bool insideDCAxy = false; - - if (cfgUseNewDCAxyCut) { - double sigmaFactor = cfgTrackPIDsettings->get(i, "maxDcaXY"); - double sigma_base = (0.0118 * std::exp(-0.6889 * ptReco) + 0.0017); - double sigma_new = sigmaFactor * sigma_base; - - insideDCAxy = (std::abs(track.dcaXY()) <= sigma_new); - } else { - 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)))); - } - - bool insideDCAz = false; - - if (cfgUseNewDCAzCut) { - - double sigmaFactorZ = cfgTrackPIDsettings->get(i, "maxDcaZ"); - - double p0 = 0.1014; - double p1 = 1.7512; - double p2 = 0.0024; - - double sigma_base_z = p0 * std::exp(-p1 * ptReco) + p2; + float ptDCA; - double sigma_new_z = sigmaFactorZ * sigma_base_z; - - insideDCAz = (std::abs(track.dcaZ()) <= sigma_new_z); + if (i == he3 || i == he4) { + ptDCA = ptReco; } else { - insideDCAz = (std::abs(track.dcaZ()) <= cfgTrackPIDsettings->get(i, "maxDcaZ")); + ptDCA = 2 * ptReco; } - if ((!insideDCAxy || !insideDCAz)) { - continue; - } - - float ptTOF = -1.0; // Default: no TOF - if (track.hasTOF()) { - ptTOF = ptReco; - } - - if (cfgTrackPIDsettings2->get(i, "fillsparsh") == 1) { - histomc.fill(HIST("hSpectraDCA"), i, particleAnti, collision.centFT0C(), - ptReco, ptTOF, decayType, track.dcaXY()); + if (pdg == particlePdgCodes.at(i)) { // He3 + histomc.fill(HIST("DCAxy_vs_pT_total"), ptDCA, track.dcaXY(), collision.centFT0C()); + if (decayType == 2) { // Transport/Material + histomc.fill(HIST("DCAxy_vs_pT_transport"), ptDCA, track.dcaXY(), collision.centFT0C()); + } else if (decayType == 1) { // Weak decay (including HF) + histomc.fill(HIST("DCAxy_vs_pT_weakdecay"), ptDCA, track.dcaXY(), collision.centFT0C()); + } + } else if (pdg == -particlePdgCodes.at(i)) { // anti-He3 + histomc.fill(HIST("DCAxy_vs_pT_anti_total"), ptDCA, track.dcaXY(), collision.centFT0C()); } } } @@ -1392,7 +1243,6 @@ struct NucleitpcPbPb { return; } - // Get pT similar to fillhmass float ptMomn; setTrackParCov(track, mTrackParCov); mTrackParCov.setPID(track.pidForTracking());