diff --git a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx index fc99ebcd30b..1eff338cec8 100644 --- a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx +++ b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx @@ -59,6 +59,7 @@ 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 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 +104,7 @@ struct PrimParticles { } }; // struct PrimParticles //---------------------------------------------------------------------------------------------------------------- -std::vector> hmass; +std::vector> hmass; std::vector> hmassnsigma; } // namespace //---------------------------------------------------------------------------------------------------------------- @@ -144,6 +145,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", 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"}; @@ -153,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"}; @@ -165,14 +168,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; @@ -192,6 +199,7 @@ struct NucleitpcPbPb { TRandom3 rand; float he3 = 4; float he4 = 5; + float processmaterial = 23; //---------------------------------------------------------------------------------------------------------------- void init(InitContext const&) { @@ -213,6 +221,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 +233,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 +323,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); ///////////////////////////////////////////////////////////////////////////////// @@ -320,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; @@ -353,12 +391,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 +408,246 @@ 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()) - continue; - if (!mcParticle.isPhysicalPrimary()) + + if (mcParticle.mcCollisionId() != mcCollision.globalIndex()) continue; + int pdgCode = mcParticle.pdgCode(); - if (std::abs(mcParticle.y()) > cfgCutRapidity && cfgrapidityRequireMC) + bool isHe3 = (std::abs(pdgCode) == particlePdgCodes.at(4)); + bool isHe4 = (std::abs(pdgCode) == particlePdgCodes.at(5)); + + 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(); + 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()); + } + } + + 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); + } + 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; - if (std::abs(mcParticle.eta()) > cfgCutEta && cfgetaRequireMC) + + 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; + } + } + } + } else if (mcParticle.has_mothers()) { + decayType = 1; + } else { + decayType = 2; + continue; + } + bool isFromWeakDecay = (decayType == 1); + + if (!mcParticle.isPhysicalPrimary() && !isFromWeakDecay) continue; - histomc.fill(HIST("histEtagen"), mcParticle.eta()); - float ptScaled = mcParticle.pt(); + + 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].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 == 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); + } + } if (pdgCode == particlePdgCodes.at(4)) { - histomc.fill(HIST("histPtgenHe3"), ptScaled); + histomc.fill(HIST("histPtgenHe3"), mcParticle.pt()); } else if (pdgCode == -particlePdgCodes.at(4)) { - histomc.fill(HIST("histPtgenAntiHe3"), ptScaled); + histomc.fill(HIST("histPtgenAntiHe3"), mcParticle.pt()); } else if (pdgCode == particlePdgCodes.at(5)) { - histomc.fill(HIST("histPtgenHe4"), ptScaled); + histomc.fill(HIST("histPtgenHe4"), mcParticle.pt()); } else if (pdgCode == -particlePdgCodes.at(5)) { - histomc.fill(HIST("histPtgenAntiHe4"), ptScaled); + histomc.fill(HIST("histPtgenAntiHe4"), mcParticle.pt()); } - } // mc track loop generated - } // mc collision loop generated - // ----------------------------- Reconstructed track loop ----------------------------- + } + } + 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(); + bool isHe3 = (std::abs(pdgCode) == particlePdgCodes.at(4)); + bool isHe4 = (std::abs(pdgCode) == particlePdgCodes.at(5)); + + if (!isHe3 && !isHe4) + continue; + if (std::abs(mcParticle.eta()) > cfgCutEta) + continue; + if (std::abs(mcParticle.y()) > cfgCutRapidity) + 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; + } + } + } + } 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) { + 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 +670,74 @@ 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) == 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; + } + } + + 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) { + 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; + } + if (!track.isPVContributor() && cfgUsePVcontributors) continue; if (!track.hasITS() && cfgITSrequire) @@ -467,12 +750,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 +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) // 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; @@ -506,7 +789,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 +801,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) == 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()); + } } - /*----------------------------------------------------------------------------------------------------------------*/ - 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 +829,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 +843,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 +929,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 +942,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 +964,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); } } //----------------------------------------------------------------------------------------------------------------