From 52b4f4090c3479fb7312a4fb397f748dbeeda3f9 Mon Sep 17 00:00:00 2001 From: Samuele Cattaruzzi <124249902+scattaru@users.noreply.github.com> Date: Wed, 3 Apr 2024 17:29:45 +0200 Subject: [PATCH] PWGHF: Ds-h correlations secondary particles (#5457) * Added secondary particles histograms and multiplicity MC in eff * clang-format * Post processing macros minor error adjustment --- PWGHF/HFC/DataModel/CorrelationTables.h | 20 +- PWGHF/HFC/Macros/DhCorrelationExtraction.cxx | 13 +- .../HFC/TableProducer/correlatorDsHadrons.cxx | 244 ++++++++++-------- PWGHF/HFC/Tasks/taskCorrelationDsHadrons.cxx | 5 + 4 files changed, 167 insertions(+), 115 deletions(-) diff --git a/PWGHF/HFC/DataModel/CorrelationTables.h b/PWGHF/HFC/DataModel/CorrelationTables.h index 44a095bb6c5..876c70ac5ac 100644 --- a/PWGHF/HFC/DataModel/CorrelationTables.h +++ b/PWGHF/HFC/DataModel/CorrelationTables.h @@ -118,14 +118,15 @@ DECLARE_SOA_TABLE(LcHadronRecoInfo, "AOD", "LCHRECOINFO", //! Lc-Hadrons pairs R // definition of columns and tables for Ds-Hadron correlation pairs namespace hf_correlation_ds_hadron { -DECLARE_SOA_COLUMN(DeltaPhi, deltaPhi, float); //! DeltaPhi between Ds and Hadrons -DECLARE_SOA_COLUMN(DeltaEta, deltaEta, float); //! DeltaEta between Ds and Hadrons -DECLARE_SOA_COLUMN(PtD, ptD, float); //! Transverse momentum of Ds -DECLARE_SOA_COLUMN(PtHadron, ptHadron, float); //! Transverse momentum of Hadron -DECLARE_SOA_COLUMN(MD, mD, float); //! Invariant mass of Ds -DECLARE_SOA_COLUMN(PoolBin, poolBin, int); //! Pool Bin for the MixedEvent -DECLARE_SOA_COLUMN(IsSignal, isSignal, bool); //! Used in MC-Rec, Ds Signal -DECLARE_SOA_COLUMN(IsPrompt, isPrompt, bool); //! Used in MC-Rec, Ds Prompt or Non-Prompt +DECLARE_SOA_COLUMN(DeltaPhi, deltaPhi, float); //! DeltaPhi between Ds and Hadrons +DECLARE_SOA_COLUMN(DeltaEta, deltaEta, float); //! DeltaEta between Ds and Hadrons +DECLARE_SOA_COLUMN(PtD, ptD, float); //! Transverse momentum of Ds +DECLARE_SOA_COLUMN(PtHadron, ptHadron, float); //! Transverse momentum of Hadron +DECLARE_SOA_COLUMN(MD, mD, float); //! Invariant mass of Ds +DECLARE_SOA_COLUMN(PoolBin, poolBin, int); //! Pool Bin for the MixedEvent +DECLARE_SOA_COLUMN(IsSignal, isSignal, bool); //! Used in MC-Rec, Ds Signal +DECLARE_SOA_COLUMN(IsPrompt, isPrompt, bool); //! Used in MC-Rec, Ds Prompt or Non-Prompt +DECLARE_SOA_COLUMN(IsPhysicalPrimary, isPhysicalPrimary, bool); //! Used in MC-Rec, primary associated particles } // namespace hf_correlation_ds_hadron DECLARE_SOA_TABLE(DsHadronPair, "AOD", "DSHPAIR", //! Ds-Hadrons pairs Informations @@ -140,7 +141,8 @@ DECLARE_SOA_TABLE(DsHadronRecoInfo, "AOD", "DSHRECOINFO", //! Ds-Hadrons pairs R aod::hf_correlation_ds_hadron::IsSignal); DECLARE_SOA_TABLE(DsHadronGenInfo, "AOD", "DSHGENINFO", //! Ds-Hadrons pairs Generated Informations - aod::hf_correlation_ds_hadron::IsPrompt); + aod::hf_correlation_ds_hadron::IsPrompt, + aod::hf_correlation_ds_hadron::IsPhysicalPrimary); // definition of columns and tables for Dplus properties namespace hf_dplus_meson diff --git a/PWGHF/HFC/Macros/DhCorrelationExtraction.cxx b/PWGHF/HFC/Macros/DhCorrelationExtraction.cxx index 59f157f9a42..2c222fe5982 100644 --- a/PWGHF/HFC/Macros/DhCorrelationExtraction.cxx +++ b/PWGHF/HFC/Macros/DhCorrelationExtraction.cxx @@ -485,14 +485,17 @@ void DhCorrelationExtraction::NormalizeMEplot(TH2D*& histoME, TH2D*& histoMEsoft // evaluate the normalization (from ALL tracks, including possible fake softpions) -> **histoME indeed includes bin1+bin2 of THnSparse, i.e. all the tracks** Double_t factorNorm = 0; - for (int in = -1; in <= 0; in++) - factorNorm += histoME->GetBinContent(bin0phi + in, bin0eta); - for (int in = -1; in <= 0; in++) - factorNorm += histoME->GetBinContent(bin0phi + in, bin0eta - 1); + for (int in = -1; in <= 0; in++) { + factorNorm += histoME->GetBinContent(bin0eta, bin0phi + in); + } + for (int in = -1; in <= 0; in++) { + factorNorm += histoME->GetBinContent(bin0eta - 1, bin0phi + in); + } factorNorm /= 4.; - if (fSubtractSoftPiME) + if (fSubtractSoftPiME) { histoME->Add(histoMEsoftPi, -1); // remove the tracks compatible with soft pion (if requested) + } // apply the normalization histoME->Scale(1. / factorNorm); diff --git a/PWGHF/HFC/TableProducer/correlatorDsHadrons.cxx b/PWGHF/HFC/TableProducer/correlatorDsHadrons.cxx index 1fc31333d65..921dd29e8f6 100644 --- a/PWGHF/HFC/TableProducer/correlatorDsHadrons.cxx +++ b/PWGHF/HFC/TableProducer/correlatorDsHadrons.cxx @@ -143,6 +143,7 @@ struct HfCorrelatorDsHadrons { Configurable numberEventsMixed{"numberEventsMixed", 5, "Number of events mixed in ME process"}; Configurable applyEfficiency{"applyEfficiency", true, "Flag for applying D-meson efficiency weights"}; Configurable yCandMax{"yCandMax", 0.8, "max. cand. rapidity"}; + Configurable yCandGenMax{"yCandGenMax", 0.5, "max. gen. cand. rapidity"}; Configurable etaTrackMax{"etaTrackMax", 0.8, "max. eta of tracks"}; Configurable dcaXYTrackMax{"dcaXYTrackMax", 1., "max. DCA_xy of tracks"}; Configurable dcaZTrackMax{"dcaZTrackMax", 1., "max. DCA_z of tracks"}; @@ -179,6 +180,7 @@ struct HfCorrelatorDsHadrons { enum AssocTrackStep { kAssocTrackStepMcGen = 0, kAssocTrackStepMcGenInAcceptance, kAssocTrackStepRecoAll, + kAssocTrackStepRecoMcMatch, kAssocTrackStepRecoPrimaries, kAssocTrackStepFake, kAssocTrackNSteps }; @@ -235,6 +237,8 @@ struct HfCorrelatorDsHadrons { registry.add("hPtCandMcRecSigNonPrompt", "Ds,Hadron candidates Non Prompt - MC Reco", {HistType::kTH1F, {axisPtD}}); registry.add("hPtCandMcRecBkg", "Ds,Hadron candidates - MC Reco", {HistType::kTH1F, {axisPtD}}); registry.add("hPtParticleAssocMcRec", "Associated Particle - MC Rec", {HistType::kTH1F, {axisPtHadron}}); + registry.add("hPtParticleAssocVsCandMcRec", "Associated Particle - MC Rec", {HistType::kTH2F, {{axisPtHadron}, {axisPtD}}}); + registry.add("hPtPrimaryParticleAssocVsCandMcRec", "Associated Particle - MC Rec", {HistType::kTH2F, {{axisPtHadron}, {axisPtD}}}); registry.add("hEtaMcRecSig", "Ds,Hadron candidates - MC Reco", {HistType::kTH1F, {axisEta}}); registry.add("hPhiMcRecSig", "Ds,Hadron candidates - MC Reco", {HistType::kTH1F, {axisPhi}}); registry.add("hEtaMcRecBkg", "Ds,Hadron candidates - MC Reco", {HistType::kTH1F, {axisEta}}); @@ -245,6 +249,7 @@ struct HfCorrelatorDsHadrons { registry.add("hMassDsVsPtMcRec", "Ds signal candidates - MC Reco", {HistType::kTH2F, {{axisMassD}, {axisPtD}}}); registry.add("hMassDsMcRecSig", "Ds signal candidates - MC Reco", {HistType::kTH2F, {{axisMassD}, {axisPtD}}}); registry.add("hMassDsMcRecBkg", "Ds background candidates - MC Reco", {HistType::kTH2F, {{axisMassD}, {axisPtD}}}); + registry.add("hFakeTracksMcRec", "Fake tracks - MC Rec", {HistType::kTH1F, {axisPtHadron}}); // Histograms for MC Gen analysis registry.add("hPtCandMcGen", "Ds,Hadron particles - MC Gen", {HistType::kTH1F, {axisPtD}}); registry.add("hPtCandMcGenPrompt", "Ds,Hadron particles - MC Gen Prompt", {HistType::kTH1F, {axisPtD}}); @@ -393,7 +398,7 @@ struct HfCorrelatorDsHadrons { track.pt(), poolBin); entryDsHadronRecoInfo(hfHelper.invMassDsToKKPi(candidate), false); - entryDsHadronGenInfo(false); + entryDsHadronGenInfo(false, false); } else if (candidate.isSelDsToPiKK() >= selectionFlagDs) { entryDsHadronPair(getDeltaPhi(track.phi(), candidate.phi()), track.eta() - candidate.eta(), @@ -401,111 +406,141 @@ struct HfCorrelatorDsHadrons { track.pt(), poolBin); entryDsHadronRecoInfo(hfHelper.invMassDsToPiKK(candidate), false); - entryDsHadronGenInfo(false); + entryDsHadronGenInfo(false, false); } - } - } + } // end track loop + } // end candidate loop } PROCESS_SWITCH(HfCorrelatorDsHadrons, processData, "Process data", true); /// Ds-Hadron correlation pair builder - for MC reco-level analysis (candidates matched to true signal only, but also the various bkg sources are studied) void processMcRec(SelCollisionsWithDs::iterator const& collision, CandDsMcReco const& candidates, - MyTracksData const& tracks) + TracksWithMc const& tracks, + aod::McParticles const&) { BinningType corrBinning{{zPoolBins, multPoolBins}, true}; - if (candidates.size() > 0) { - registry.fill(HIST("hZVtx"), collision.posZ()); - registry.fill(HIST("hMultFT0M"), collision.multFT0M()); - int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M())); - registry.fill(HIST("hCollisionPoolBin"), poolBin); - - // MC reco level - bool isDsPrompt = false; - bool isDsSignal = false; - float multiplicityFT0M = collision.multFT0M(); - for (const auto& candidate : candidates) { - // prompt and non-prompt division - isDsPrompt = candidate.originMcRec() == RecoDecay::OriginType::Prompt; - // Ds Signal - isDsSignal = std::abs(candidate.flagMcMatchRec()) == 1 << aod::hf_cand_3prong::DecayType::DsToKKPi; + registry.fill(HIST("hZVtx"), collision.posZ()); + registry.fill(HIST("hMultFT0M"), collision.multFT0M()); + int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M())); + registry.fill(HIST("hCollisionPoolBin"), poolBin); - if (std::abs(hfHelper.yDs(candidate)) > yCandMax || candidate.pt() < ptCandMin || candidate.pt() > ptCandMax) { - continue; + // MC reco level + bool isDsPrompt = false; + bool isDsSignal = false; + bool isAlreadyFilledEvent = false; + float multiplicityFT0M = collision.multFT0M(); + for (const auto& candidate : candidates) { + // prompt and non-prompt division + isDsPrompt = candidate.originMcRec() == RecoDecay::OriginType::Prompt; + // Ds Signal + isDsSignal = std::abs(candidate.flagMcMatchRec()) == 1 << aod::hf_cand_3prong::DecayType::DsToKKPi; + + if (std::abs(hfHelper.yDs(candidate)) > yCandMax || candidate.pt() < ptCandMin || candidate.pt() > ptCandMax) { + continue; + } + + double efficiencyWeightD = 1.; + if (applyEfficiency) { + efficiencyWeightD = 1. / efficiencyD->at(o2::analysis::findBin(binsPtEfficiencyD, candidate.pt())); + } + if (isDsSignal) { + fillHistoMcRecSig(candidate, multiplicityFT0M); + // DsToKKPi and DsToPiKK division + if (candidate.isSelDsToKKPi() >= selectionFlagDs) { + registry.fill(HIST("hMassDsMCRec"), hfHelper.invMassDsToKKPi(candidate), efficiencyWeightD); + registry.fill(HIST("hMassDsMCRecSig"), hfHelper.invMassDsToKKPi(candidate), candidate.pt(), efficiencyWeightD); + registry.fill(HIST("hMassDsVsPtMCRec"), hfHelper.invMassDsToKKPi(candidate), candidate.pt(), efficiencyWeightD); + registry.fill(HIST("hSelectionStatusDsToKKPi"), candidate.isSelDsToKKPi()); + } + if (candidate.isSelDsToPiKK() >= selectionFlagDs) { + registry.fill(HIST("hMassDsMCRec"), hfHelper.invMassDsToPiKK(candidate), efficiencyWeightD); + registry.fill(HIST("hMassDsMCRecSig"), hfHelper.invMassDsToPiKK(candidate), candidate.pt(), efficiencyWeightD); + registry.fill(HIST("hMassDsVsPtMCRec"), hfHelper.invMassDsToPiKK(candidate), candidate.pt(), efficiencyWeightD); + registry.fill(HIST("hSelectionStatusDsToPiKK"), candidate.isSelDsToPiKK()); + } + if (candidate.isSelDsToKKPi() >= selectionFlagDs && candidate.isSelDsToPiKK() >= selectionFlagDs) { + registry.fill(HIST("hCountSelectionStatusDsToKKPiAndToPiKK"), 0.); + } + } else { + fillHistoMcRecBkg(candidate); + // DsToKKPi and DsToPiKK division + if (candidate.isSelDsToKKPi() >= selectionFlagDs) { + registry.fill(HIST("hMassDsMCRec"), hfHelper.invMassDsToKKPi(candidate), efficiencyWeightD); + registry.fill(HIST("hMassDsMCRecBkg"), hfHelper.invMassDsToKKPi(candidate), candidate.pt(), efficiencyWeightD); + registry.fill(HIST("hMassDsVsPtMCRec"), hfHelper.invMassDsToKKPi(candidate), candidate.pt(), efficiencyWeightD); + registry.fill(HIST("hSelectionStatusDsToKKPiMcRec"), candidate.isSelDsToKKPi()); + } else if (candidate.isSelDsToPiKK() >= selectionFlagDs) { + registry.fill(HIST("hMassDsMCRec"), hfHelper.invMassDsToPiKK(candidate), efficiencyWeightD); + registry.fill(HIST("hMassDsMCRecBkg"), hfHelper.invMassDsToPiKK(candidate), candidate.pt(), efficiencyWeightD); + registry.fill(HIST("hMassDsVsPtMCRec"), hfHelper.invMassDsToPiKK(candidate), candidate.pt(), efficiencyWeightD); + registry.fill(HIST("hSelectionStatusDsToPiKKMcRec"), candidate.isSelDsToPiKK()); } + } - double efficiencyWeightD = 1.; - if (applyEfficiency) { - efficiencyWeightD = 1. / efficiencyD->at(o2::analysis::findBin(binsPtEfficiencyD, candidate.pt())); + // Ds-Hadron correlation dedicated section + // if the candidate is selected as Ds, search for Hadron and evaluate correlations + for (const auto& track : tracks) { + // Removing Ds daughters by checking track indices + if ((candidate.prong0Id() == track.globalIndex()) || (candidate.prong1Id() == track.globalIndex()) || (candidate.prong2Id() == track.globalIndex())) { + continue; } - if (isDsSignal) { - fillHistoMcRecSig(candidate, multiplicityFT0M); - // DsToKKPi and DsToPiKK division - if (candidate.isSelDsToKKPi() >= selectionFlagDs) { - registry.fill(HIST("hMassDsMCRec"), hfHelper.invMassDsToKKPi(candidate), efficiencyWeightD); - registry.fill(HIST("hMassDsMCRecSig"), hfHelper.invMassDsToKKPi(candidate), candidate.pt(), efficiencyWeightD); - registry.fill(HIST("hMassDsVsPtMCRec"), hfHelper.invMassDsToKKPi(candidate), candidate.pt(), efficiencyWeightD); - registry.fill(HIST("hSelectionStatusDsToKKPi"), candidate.isSelDsToKKPi()); - } - if (candidate.isSelDsToPiKK() >= selectionFlagDs) { - registry.fill(HIST("hMassDsMCRec"), hfHelper.invMassDsToPiKK(candidate), efficiencyWeightD); - registry.fill(HIST("hMassDsMCRecSig"), hfHelper.invMassDsToPiKK(candidate), candidate.pt(), efficiencyWeightD); - registry.fill(HIST("hMassDsVsPtMCRec"), hfHelper.invMassDsToPiKK(candidate), candidate.pt(), efficiencyWeightD); - registry.fill(HIST("hSelectionStatusDsToPiKK"), candidate.isSelDsToPiKK()); - } - if (candidate.isSelDsToKKPi() >= selectionFlagDs && candidate.isSelDsToPiKK() >= selectionFlagDs) { - registry.fill(HIST("hCountSelectionStatusDsToKKPiAndToPiKK"), 0.); + bool isPhysicalPrimary = false; + // DsToKKPi and DsToPiKK division + if (candidate.isSelDsToKKPi() >= selectionFlagDs) { + entryDsHadronPair(getDeltaPhi(track.phi(), candidate.phi()), + track.eta() - candidate.eta(), + candidate.pt(), + track.pt(), + poolBin); + entryDsHadronRecoInfo(hfHelper.invMassDsToKKPi(candidate), isDsSignal); + if (track.has_mcParticle()) { + auto mcParticle = track.template mcParticle_as(); + isPhysicalPrimary = mcParticle.isPhysicalPrimary(); + entryDsHadronGenInfo(isDsPrompt, isPhysicalPrimary); + } else { + entryDsHadronGenInfo(isDsPrompt, isPhysicalPrimary); + registry.fill(HIST("hFakeTracksMcRec"), track.pt()); } - } else { - fillHistoMcRecBkg(candidate); - // DsToKKPi and DsToPiKK division - if (candidate.isSelDsToKKPi() >= selectionFlagDs) { - registry.fill(HIST("hMassDsMCRec"), hfHelper.invMassDsToKKPi(candidate), efficiencyWeightD); - registry.fill(HIST("hMassDsMCRecBkg"), hfHelper.invMassDsToKKPi(candidate), candidate.pt(), efficiencyWeightD); - registry.fill(HIST("hMassDsVsPtMCRec"), hfHelper.invMassDsToKKPi(candidate), candidate.pt(), efficiencyWeightD); - registry.fill(HIST("hSelectionStatusDsToKKPiMcRec"), candidate.isSelDsToKKPi()); - } else if (candidate.isSelDsToPiKK() >= selectionFlagDs) { - registry.fill(HIST("hMassDsMCRec"), hfHelper.invMassDsToPiKK(candidate), efficiencyWeightD); - registry.fill(HIST("hMassDsMCRecBkg"), hfHelper.invMassDsToPiKK(candidate), candidate.pt(), efficiencyWeightD); - registry.fill(HIST("hMassDsVsPtMCRec"), hfHelper.invMassDsToPiKK(candidate), candidate.pt(), efficiencyWeightD); - registry.fill(HIST("hSelectionStatusDsToPiKKMcRec"), candidate.isSelDsToPiKK()); + // for secondary particle fraction estimation + if (!isAlreadyFilledEvent) { + registry.fill(HIST("hPtParticleAssocVsCandMcRec"), track.pt(), candidate.pt()); + if (isPhysicalPrimary) { + registry.fill(HIST("hPtPrimaryParticleAssocVsCandMcRec"), track.pt(), candidate.pt()); + } } - } - - // Ds-Hadron correlation dedicated section - // if the candidate is selected as Ds, search for Hadron and evaluate correlations - for (const auto& track : tracks) { - // Removing Ds daughters by checking track indices - if ((candidate.prong0Id() == track.globalIndex()) || (candidate.prong1Id() == track.globalIndex()) || (candidate.prong2Id() == track.globalIndex())) { - continue; + } else if (candidate.isSelDsToPiKK() >= selectionFlagDs) { + entryDsHadronPair(getDeltaPhi(track.phi(), candidate.phi()), + track.eta() - candidate.eta(), + candidate.pt(), + track.pt(), + poolBin); + entryDsHadronRecoInfo(hfHelper.invMassDsToPiKK(candidate), isDsSignal); + if (track.has_mcParticle()) { + auto mcParticle = track.template mcParticle_as(); + isPhysicalPrimary = mcParticle.isPhysicalPrimary(); + entryDsHadronGenInfo(isDsPrompt, isPhysicalPrimary); + } else { + entryDsHadronGenInfo(isDsPrompt, false); + registry.fill(HIST("hFakeTracksMcRec"), track.pt()); } - registry.fill(HIST("hPtParticleAssocMcRec"), track.pt()); - // DsToKKPi and DsToPiKK division - if (candidate.isSelDsToKKPi() >= selectionFlagDs) { - entryDsHadronPair(getDeltaPhi(track.phi(), candidate.phi()), - track.eta() - candidate.eta(), - candidate.pt(), - track.pt(), - poolBin); - entryDsHadronRecoInfo(hfHelper.invMassDsToKKPi(candidate), isDsSignal); - entryDsHadronGenInfo(isDsPrompt); - } else if (candidate.isSelDsToPiKK() >= selectionFlagDs) { - entryDsHadronPair(getDeltaPhi(track.phi(), candidate.phi()), - track.eta() - candidate.eta(), - candidate.pt(), - track.pt(), - poolBin); - entryDsHadronRecoInfo(hfHelper.invMassDsToPiKK(candidate), isDsSignal); - entryDsHadronGenInfo(isDsPrompt); + // for secondary particle fraction estimation + if (!isAlreadyFilledEvent) { + registry.fill(HIST("hPtParticleAssocVsCandMcRec"), track.pt(), candidate.pt()); + if (isPhysicalPrimary) { + registry.fill(HIST("hPtPrimaryParticleAssocVsCandMcRec"), track.pt(), candidate.pt()); + } } } - } - } + } // end track loop + isAlreadyFilledEvent = true; + } // end candidate loop } PROCESS_SWITCH(HfCorrelatorDsHadrons, processMcRec, "Process MC Reco mode", false); /// Ds-Hadron correlation - for calculating candidate reconstruction efficiency using MC reco-level analysis void processMcCandEfficiency(soa::Join const&, + soa::Join const&, CandDsMcGen const& mcParticles, CandDsMcReco const& candidates, aod::TracksWMc const&) @@ -513,25 +548,27 @@ struct HfCorrelatorDsHadrons { auto hCandidates = registry.get(HIST("hCandidates")); /// Gen loop - double multiplicity = 1.; + float multiplicity = -1.; for (auto& mcParticle : mcParticles) { // generated candidates if (std::abs(mcParticle.pdgCode()) == Pdg::kDS) { + auto mcCollision = mcParticle.template mcCollision_as>(); + multiplicity = mcCollision.multMCFT0A() + mcCollision.multMCFT0C(); // multFT0M = multFt0A + multFT0C hCandidates->Fill(kCandidateStepMcGenAll, mcParticle.pt(), multiplicity, mcParticle.originMcGen()); if (std::abs(mcParticle.flagMcMatchGen()) == 1 << aod::hf_cand_3prong::DecayType::DsToKKPi) { hCandidates->Fill(kCandidateStepMcGenDsToKKPi, mcParticle.pt(), multiplicity, mcParticle.originMcGen()); auto yDs = RecoDecay::y(std::array{mcParticle.px(), mcParticle.py(), mcParticle.pz()}, o2::constants::physics::MassDS); - if (std::abs(yDs) <= yCandMax) { + if (std::abs(yDs) <= yCandGenMax) { hCandidates->Fill(kCandidateStepMcCandInAcceptance, mcParticle.pt(), multiplicity, mcParticle.originMcGen()); } - bool inAcceptance = true; + bool isDaughterInAcceptance = true; auto daughters = mcParticle.template daughters_as(); for (const auto& daughter : daughters) { if (daughter.pt() < ptDaughterMin || std::abs(daughter.eta()) > etaTrackMax) { - inAcceptance = false; + isDaughterInAcceptance = false; } } - if (inAcceptance) { + if (isDaughterInAcceptance) { hCandidates->Fill(kCandidateStepMcDaughtersInAcceptance, mcParticle.pt(), multiplicity, mcParticle.originMcGen()); fillHistoMcGen(mcParticle); } @@ -561,17 +598,21 @@ struct HfCorrelatorDsHadrons { /// Ds-Hadron correlation - for calculating associated particle tracking efficiency using MC reco-level analysis void processMcTrackEfficiency(soa::Join const&, + soa::Join const&, CandDsMcGen const& mcParticles, TracksWithMc const& tracksData) { auto hAssocTracks = registry.get(HIST("hAssocTracks")); /// Gen loop - double multiplicity = 1.; - double posZ = 0.; + float multiplicity = -1.; + float posZ = -20.; for (auto& mcParticle : mcParticles) { // generated tracks if (mcParticle.isPhysicalPrimary() && ((std::abs(mcParticle.pdgCode()) == kElectron) || (std::abs(mcParticle.pdgCode()) == kMuonMinus) || (std::abs(mcParticle.pdgCode()) == kPiPlus) || (std::abs(mcParticle.pdgCode()) == kKPlus) || (std::abs(mcParticle.pdgCode()) == kProton))) { + auto mcCollision = mcParticle.template mcCollision_as>(); + multiplicity = mcCollision.multMCFT0A() + mcCollision.multMCFT0C(); // multFT0M = multFt0A + multFT0C + posZ = mcCollision.posZ(); hAssocTracks->Fill(kAssocTrackStepMcGen, mcParticle.eta(), mcParticle.pt(), multiplicity, posZ); if (mcParticle.pt() > ptTrackMin && std::abs(mcParticle.eta()) < etaTrackMax) { hAssocTracks->Fill(kAssocTrackStepMcGenInAcceptance, mcParticle.eta(), mcParticle.pt(), multiplicity, posZ); @@ -588,9 +629,10 @@ struct HfCorrelatorDsHadrons { posZ = collision.posZ(); registry.fill(HIST("hZVtx"), posZ); registry.fill(HIST("hMultFT0M"), multiplicity); + hAssocTracks->Fill(kAssocTrackStepRecoAll, track.eta(), track.pt(), multiplicity, posZ); if (track.has_mcParticle()) { auto mcParticle = track.template mcParticle_as(); - hAssocTracks->Fill(kAssocTrackStepRecoAll, mcParticle.eta(), mcParticle.pt(), multiplicity, posZ); + hAssocTracks->Fill(kAssocTrackStepRecoMcMatch, mcParticle.eta(), mcParticle.pt(), multiplicity, posZ); if (mcParticle.isPhysicalPrimary()) { hAssocTracks->Fill(kAssocTrackStepRecoPrimaries, mcParticle.eta(), mcParticle.pt(), multiplicity, posZ); registry.fill(HIST("hPtParticleAssocMcRec"), track.pt()); @@ -631,7 +673,7 @@ struct HfCorrelatorDsHadrons { } if (std::abs(particle.flagMcMatchGen()) == 1 << aod::hf_cand_3prong::DecayType::DsToKKPi) { double yD = RecoDecay::y(std::array{particle.px(), particle.py(), particle.pz()}, MassDS); // TODO - if (yCandMax >= 0. && std::abs(yD) > yCandMax) { + if (yCandGenMax >= 0. && std::abs(yD) > yCandGenMax) { continue; } if (ptCandMin >= 0. && particle.pt() < ptCandMin) { @@ -660,7 +702,7 @@ struct HfCorrelatorDsHadrons { particleAssoc.pt(), poolBin); entryDsHadronRecoInfo(MassDS, true); - entryDsHadronGenInfo(isDsPrompt); + entryDsHadronGenInfo(isDsPrompt, particleAssoc.isPhysicalPrimary()); } } } @@ -672,15 +714,17 @@ struct HfCorrelatorDsHadrons { CandDsData const& candidates, MyTracksData const& tracks) { - std::cout << "1" << std::endl; BinningType corrBinning{{zPoolBins, multPoolBins}, true}; + for (const auto& collision : collisions) { + registry.fill(HIST("hMultFT0M"), collision.multFT0M()); + registry.fill(HIST("hZVtx"), collision.posZ()); + } for (const auto& candidate : candidates) { if (std::abs(hfHelper.yDs(candidate)) > yCandMax || candidate.pt() < ptCandMin || candidate.pt() > ptCandMax) { continue; } fillHisto(candidate); } - std::cout << "2" << std::endl; auto tracksTuple = std::make_tuple(candidates, tracks); Pair pairData{corrBinning, numberEventsMixed, -1, collisions, tracksTuple, &cache}; @@ -691,8 +735,6 @@ struct HfCorrelatorDsHadrons { // LOGF(info, "Mixed event collisions: Index = (%d, %d), tracks Size: (%d, %d), Z Vertex: (%f, %f), Pool Bin: (%d, %d)", c1.globalIndex(), c2.globalIndex(), tracks1.size(), tracks2.size(), c1.posZ(), c2.posZ(), corrBinning.getBin(std::make_tuple(c1.posZ(), c1.multFT0M())), corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFT0M()))); int poolBin = corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFT0M())); int poolBinDs = corrBinning.getBin(std::make_tuple(c1.posZ(), c1.multFT0M())); - registry.fill(HIST("hMultFT0M"), c1.multFT0M()); - registry.fill(HIST("hZVtx"), c1.posZ()); registry.fill(HIST("hTracksPoolBin"), poolBin); registry.fill(HIST("hDsPoolBin"), poolBinDs); for (const auto& [cand, pAssoc] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { @@ -712,7 +754,7 @@ struct HfCorrelatorDsHadrons { pAssoc.pt(), poolBin); entryDsHadronRecoInfo(hfHelper.invMassDsToKKPi(cand), false); - entryDsHadronGenInfo(false); + entryDsHadronGenInfo(false, false); } else if (cand.isSelDsToPiKK() >= selectionFlagDs) { // LOGF(info, "Mixed event tracks pair: (%d, %d) from events (%d, %d), track event: (%d, %d), PiKK", cand.index(), pAssoc.index(), c1.index(), c2.index(), cand.collision().index(), pAssoc.collision().index()); entryDsHadronPair(getDeltaPhi(pAssoc.phi(), cand.phi()), @@ -721,7 +763,7 @@ struct HfCorrelatorDsHadrons { pAssoc.pt(), poolBin); entryDsHadronRecoInfo(hfHelper.invMassDsToPiKK(cand), false); - entryDsHadronGenInfo(false); + entryDsHadronGenInfo(false, false); } } } @@ -779,7 +821,7 @@ struct HfCorrelatorDsHadrons { pAssoc.pt(), poolBin); entryDsHadronRecoInfo(hfHelper.invMassDsToKKPi(candidate), isDsSignal); - entryDsHadronGenInfo(isDsPrompt); + entryDsHadronGenInfo(isDsPrompt, false); } else if (candidate.isSelDsToPiKK() >= selectionFlagDs) { entryDsHadronPair(getDeltaPhi(pAssoc.phi(), candidate.phi()), pAssoc.eta() - candidate.eta(), @@ -787,7 +829,7 @@ struct HfCorrelatorDsHadrons { pAssoc.pt(), poolBin); entryDsHadronRecoInfo(hfHelper.invMassDsToPiKK(candidate), isDsSignal); - entryDsHadronGenInfo(isDsPrompt); + entryDsHadronGenInfo(isDsPrompt, false); } } } diff --git a/PWGHF/HFC/Tasks/taskCorrelationDsHadrons.cxx b/PWGHF/HFC/Tasks/taskCorrelationDsHadrons.cxx index d80fc17c6bc..57a1171a2fa 100644 --- a/PWGHF/HFC/Tasks/taskCorrelationDsHadrons.cxx +++ b/PWGHF/HFC/Tasks/taskCorrelationDsHadrons.cxx @@ -71,6 +71,7 @@ struct HfTaskCorrelationDsHadrons { registry.add("hDeltaEtaPtIntSignalRegionMcRec", "Ds-h deltaEta signal region MC reco", {HistType::kTH1F, {axisDetlaEta}}); registry.add("hDeltaPhiPtIntSignalRegionMcRec", "Ds-h deltaPhi signal region MC reco", {HistType::kTH1F, {axisDetlaPhi}}); registry.add("hCorrel2DVsPtSignalRegionMcRec", "Ds-h correlations signal region MC reco", {HistType::kTHnSparseD, {{axisDetlaPhi}, {axisDetlaEta}, {axisPtD}, {axisPtHadron}, {axisPoolBin}}}); + registry.add("hCorrel2DVsPtPhysicalPrimaryMcRec", "Ds-h correlations signal region (only true primary particles) MC reco", {HistType::kTHnSparseD, {{axisDetlaPhi}, {axisDetlaEta}, {axisPtD}, {axisPtHadron}}}); registry.add("hCorrel2DVsPtSignalRegionMcRecPromptDivision", "Ds-h correlations signal region Prompt-NonPrompt MC reco", {HistType::kTHnSparseD, {{axisDetlaPhi}, {axisDetlaEta}, {axisPtD}, {axisPtHadron}, {axisDsPrompt}, {axisPoolBin}}}); registry.add("hDeltaEtaPtIntSidebandsMcRec", "Ds-h deltaEta sideband region MC reco", {HistType::kTH1F, {axisDetlaEta}}); registry.add("hDeltaPhiPtIntSidebandsMcRec", "Ds-h deltaPhi sideband region MC reco", {HistType::kTH1F, {axisDetlaPhi}}); @@ -135,6 +136,7 @@ struct HfTaskCorrelationDsHadrons { int poolBin = pairEntry.poolBin(); double massD = pairEntry.mD(); int statusDsPrompt = static_cast(pairEntry.isPrompt()); + bool isPhysicalPrimary = pairEntry.isPhysicalPrimary(); int ptBinD = o2::analysis::findBin(binsPtD, ptD); double efficiencyWeight = 1.; @@ -146,6 +148,9 @@ struct HfTaskCorrelationDsHadrons { registry.fill(HIST("hCorrel2DVsPtSignalRegionMcRec"), deltaPhi, deltaEta, ptD, ptHadron, poolBin, efficiencyWeight); registry.fill(HIST("hDeltaEtaPtIntSignalRegionMcRec"), deltaEta, efficiencyWeight); registry.fill(HIST("hDeltaPhiPtIntSignalRegionMcRec"), deltaPhi, efficiencyWeight); + if (isPhysicalPrimary) { + registry.fill(HIST("hCorrel2DVsPtPhysicalPrimaryMcRec"), deltaPhi, deltaEta, ptD, ptHadron, efficiencyWeight); + } // prompt and non-prompt division if (pairEntry.isSignal()) { registry.fill(HIST("hCorrel2DVsPtSignalRegionMcRecPromptDivision"), deltaPhi, deltaEta, ptD, ptHadron, statusDsPrompt, poolBin, efficiencyWeight);