From 40a8cfd5923f99c465a2cb093373b9899a482ed1 Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Tue, 26 Jul 2022 17:33:41 +0200 Subject: [PATCH 1/3] PWGHF: add information about tracks with ambiguous collision association --- PWGHF/Tasks/HFMCValidation.cxx | 97 +++++++++++++++++++++++++++++----- 1 file changed, 85 insertions(+), 12 deletions(-) diff --git a/PWGHF/Tasks/HFMCValidation.cxx b/PWGHF/Tasks/HFMCValidation.cxx index 59fa5e755e5..d44cd34de57 100644 --- a/PWGHF/Tasks/HFMCValidation.cxx +++ b/PWGHF/Tasks/HFMCValidation.cxx @@ -39,6 +39,63 @@ static const std::array particleNames = {"Dplus", "D static const std::array originNames = {"Prompt", "NonPrompt"}; } // namespace +/// Add relevant information about ambiguous tracks +namespace o2::aod +{ +// Columns to store the DCA to the primary vertex +DECLARE_SOA_COLUMN(IsAmbiguousTrack, isAmbiguousTrack, bool); //! +DECLARE_SOA_SELF_ARRAY_INDEX_COLUMN(AmbiguousCollisionIndices, ambiguousCollisionIndices); //! +DECLARE_SOA_TABLE(TracksWithAmbiguousCollisionInfo, "AOD", "TRACKSWAMBINFO", //! + IsAmbiguousTrack, + AmbiguousCollisionIndicesIds); +} // namespace o2::aod + +struct AddAmbiguousTrackInfo { + + Produces trackWithAmbiguousInfo; + using TracksWithSel = soa::Join; + + void process(TracksWithSel const& tracks, + aod::AmbiguousTracks const& ambitracks, + aod::Collisions const& collisions, + aod::BCs const& bcs) + { + // loop over ambiguous tracks + std::vector trackIndices{}; + std::vector ambTrackIndices{}; + for (auto& ambitrack : ambitracks) { + auto track = ambitrack.track_as(); // Obtain the corresponding track + if (track.isGlobalTrackWoDCA()) { // add info only for global tracks + trackIndices.push_back(track.globalIndex()); + ambTrackIndices.push_back(ambitrack.globalIndex()); + } + } + // loop over tracks + for (auto& track : tracks) { + std::vector collIndices{}; + bool isAmbiguous = false; + if (track.isGlobalTrackWoDCA()) { // add info only for global tracks + auto trackIdx = track.globalIndex(); + auto iter = std::find(trackIndices.begin(), trackIndices.end(), trackIdx); + if (iter != trackIndices.end()) { + isAmbiguous = true; + auto ambitrack = ambitracks.rawIteratorAt(ambTrackIndices[std::distance(trackIndices.begin(), iter)]); + for (auto& collision : collisions) { + uint64_t mostProbableBC = collision.bc().globalBC(); + for (auto& bc : ambitrack.bc()) { + if (bc.globalBC() == mostProbableBC) { + collIndices.push_back(collision.globalIndex()); + break; + } + } + } + } + } + trackWithAmbiguousInfo(isAmbiguous, collIndices); + } + } +}; + /// Gen Level Validation /// /// - Number of HF quarks produced per collision @@ -57,8 +114,8 @@ struct ValidationGenLevel { Configurable xVertexMax{"xVertexMax", 100., "max. x of generated primary vertex [cm]"}; Configurable yVertexMin{"yVertexMin", -100., "min. y of generated primary vertex [cm]"}; Configurable yVertexMax{"yVertexMax", 100., "max. y of generated primary vertex [cm]"}; - Configurable zVertexMin{"zVertexMin", -10., "min. z of generated primary vertex [cm]"}; - Configurable zVertexMax{"zVertexMax", 10., "max. z of generated primary vertex [cm]"}; + Configurable zVertexMin{"zVertexMin", -100., "min. z of generated primary vertex [cm]"}; + Configurable zVertexMax{"zVertexMax", 100., "max. z of generated primary vertex [cm]"}; std::shared_ptr hPromptCharmHadronsPtDistr, hPromptCharmHadronsYDistr, hNonPromptCharmHadronsPtDistr, hNonPromptCharmHadronsYDistr, hPromptCharmHadronsDecLenDistr, hNonPromptCharmHadronsDecLenDistr, hQuarkPerEvent; @@ -274,20 +331,23 @@ struct ValidationRecLevel { std::array, 3>, 2>, nCharmHadrons> histPtDau, histEtaDau, histImpactParameterDau; std::array, 2>, nCharmHadrons> histPtReco; std::array, 4> histOriginTracks; + std::shared_ptr histAmbiguousTracks, histTracks; HistogramRegistry registry{"registry", {}}; void init(o2::framework::InitContext&) { - histOriginTracks[0] = registry.add("histOriginNonAssociatedTracks", ";origin;#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {{4, -1.5, 2.5}, {50, 0., 10.}, {40, -1., 1.}, {200, -1., 1.}}); // tracks not associated to any collision - histOriginTracks[1] = registry.add("histOriginAssociatedTracks", ";origin;#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {{4, -1.5, 2.5}, {50, 0., 10.}, {40, -1., 1.}, {200, -1., 1.}}); // tracks associasted to a collision - histOriginTracks[2] = registry.add("histOriginGoodAssociatedTracks", ";origin;#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {{4, -1.5, 2.5}, {50, 0., 10.}, {40, -1., 1.}, {200, -1., 1.}}); // tracks associated to the correct collision (considering the MC collision index) - histOriginTracks[3] = registry.add("histOriginWrongAssociatedTracks", ";origin;#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {{4, -1.5, 2.5}, {50, 0., 10.}, {40, -1., 1.}, {200, -1., 1.}}); // tracks associated to the wrong collision (considering the MC collision index) + histOriginTracks[0] = registry.add("histOriginNonAssociatedTracks", ";origin;#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm); isContributor", HistType::kTHnSparseF, {{4, -1.5, 2.5}, {50, 0., 10.}, {40, -1., 1.}, {200, -1., 1.}, {2, -0.5, 1.5}}); // tracks not associated to any collision + histOriginTracks[1] = registry.add("histOriginAssociatedTracks", ";origin;#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm); isContributor", HistType::kTHnSparseF, {{4, -1.5, 2.5}, {50, 0., 10.}, {40, -1., 1.}, {200, -1., 1.}, {2, -0.5, 1.5}}); // tracks associasted to a collision + histOriginTracks[2] = registry.add("histOriginGoodAssociatedTracks", ";origin;#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm); isContributor", HistType::kTHnSparseF, {{4, -1.5, 2.5}, {50, 0., 10.}, {40, -1., 1.}, {200, -1., 1.}, {2, -0.5, 1.5}}); // tracks associated to the correct collision considering only first reco collision (based on the MC collision index) + histOriginTracks[3] = registry.add("histOriginGoodAssociatedTracksAmbiguous", ";origin;#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm); isContributor", HistType::kTHnSparseF, {{4, -1.5, 2.5}, {50, 0., 10.}, {40, -1., 1.}, {200, -1., 1.}, {2, -0.5, 1.5}}); // tracks associated to the correct collision considering all ambiguous reco collisions (based on the MC collision index) for (std::size_t iHist{0}; iHist < histOriginTracks.size(); ++iHist) { histOriginTracks[iHist]->GetAxis(0)->SetBinLabel(1, "no MC particle"); histOriginTracks[iHist]->GetAxis(0)->SetBinLabel(2, "no quark"); histOriginTracks[iHist]->GetAxis(0)->SetBinLabel(3, "charm"); histOriginTracks[iHist]->GetAxis(0)->SetBinLabel(4, "beauty"); } + histAmbiguousTracks = registry.add("histAmbiguousTracks", "Tracks that are ambiguous vs. origin;#it{p}_{T}^{reco} (GeV/#it{c});entries", HistType::kTH2F, {{4, -1.5, 2.5}, {50, 0., 10.}}); + histTracks = registry.add("histTracks", "Tracks vs. origin;#it{p}_{T}^{reco} (GeV/#it{c});entries", HistType::kTH2F, {{4, -1.5, 2.5}, {50, 0., 10.}}); for (auto iHad = 0; iHad < nCharmHadrons; ++iHad) { histDeltaPt[iHad] = registry.add(Form("histDeltaPt%s", particleNames[iHad].data()), Form("Pt difference reco - MC %s; #it{p}_{T}^{reco} - #it{p}_{T}^{gen} (GeV/#it{c}); entries", labels[iHad].data()), HistType::kTH1F, {{2000, -1., 1.}}); histDeltaPx[iHad] = registry.add(Form("histDeltaPx%s", particleNames[iHad].data()), Form("Px difference reco - MC %s; #it{p}_{x}^{reco} - #it{p}_{x}^{gen} (GeV/#it{c}); entries", labels[iHad].data()), HistType::kTH1F, {{2000, -1., 1.}}); @@ -311,7 +371,7 @@ struct ValidationRecLevel { using HfCandProng2WithMCRec = soa::Join; using HfCandProng3WithMCRec = soa::Join; using CollisionsWithMCLabels = soa::Join; - using TracksWithSel = soa::Join; + using TracksWithSel = soa::Join; void process(HfCandProng2WithMCRec const& cand2Prongs, HfCandProng3WithMCRec const& cand3Prongs, TracksWithSel const& tracks, aod::McParticles const& particlesMC, aod::McCollisions const& mcCollisions, CollisionsWithMCLabels const& collisions) { @@ -324,20 +384,32 @@ struct ValidationRecLevel { if (track.has_mcParticle()) { auto particle = track.mcParticle(); // get corresponding MC particle to check origin auto origin = RecoDecay::getCharmHadronOrigin(particlesMC, particle, true); + histTracks->Fill(origin, track.pt()); + if (track.isAmbiguousTrack()) { + histAmbiguousTracks->Fill(origin, track.pt()); + } float deltaZ = -999.f; if (index) { auto collision = track.collision_as(); auto mcCollision = particle.mcCollision(); deltaZ = collision.posZ() - mcCollision.posZ(); - uint index2 = 2; if (collision.mcCollisionId() == particle.mcCollisionId()) { - index2 = 1; + histOriginTracks[index + 1]->Fill(origin, track.pt(), track.eta(), deltaZ, track.isPVContributor()); + } else { // if the most probable collision is not the good one, check if the tracks is ambiguous + if (track.isAmbiguousTrack()) { + for (auto& collIdx : track.ambiguousCollisionIndicesIds()) { + auto ambCollision = collisions.rawIteratorAt(collIdx); + if (ambCollision.mcCollisionId() == particle.mcCollisionId()) { + histOriginTracks[index + 2]->Fill(origin, track.pt(), track.eta(), deltaZ, track.isPVContributor()); + break; + } + } + } } - histOriginTracks[index + index2]->Fill(origin, track.pt(), track.eta(), deltaZ); } - histOriginTracks[index]->Fill(origin, track.pt(), track.eta(), deltaZ); + histOriginTracks[index]->Fill(origin, track.pt(), track.eta(), deltaZ, track.isPVContributor()); } else { - histOriginTracks[index]->Fill(-1.f, track.pt(), track.eta(), -999.f); + histOriginTracks[index]->Fill(-1.f, track.pt(), track.eta(), -999.f, track.isPVContributor()); } } @@ -475,6 +547,7 @@ struct ValidationRecLevel { WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { WorkflowSpec workflow{ + adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc, TaskName{"hf-mc-validation-gen"}), adaptAnalysisTask(cfgc, TaskName{"hf-mc-validation-rec"})}; return workflow; From cac151e1da5bf70a861d964c2c0751e17a1c5e48 Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Tue, 26 Jul 2022 17:36:40 +0200 Subject: [PATCH 2/3] Fix comment --- PWGHF/Tasks/HFMCValidation.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/Tasks/HFMCValidation.cxx b/PWGHF/Tasks/HFMCValidation.cxx index d44cd34de57..cfd81af0035 100644 --- a/PWGHF/Tasks/HFMCValidation.cxx +++ b/PWGHF/Tasks/HFMCValidation.cxx @@ -42,7 +42,7 @@ static const std::array originNames = {"Prompt", "NonPrompt"}; /// Add relevant information about ambiguous tracks namespace o2::aod { -// Columns to store the DCA to the primary vertex +// Columns to store the information about ambiguous tracks joinable with the track table DECLARE_SOA_COLUMN(IsAmbiguousTrack, isAmbiguousTrack, bool); //! DECLARE_SOA_SELF_ARRAY_INDEX_COLUMN(AmbiguousCollisionIndices, ambiguousCollisionIndices); //! DECLARE_SOA_TABLE(TracksWithAmbiguousCollisionInfo, "AOD", "TRACKSWAMBINFO", //! From b49f78a2369a4cb60d7c2d2aed3378cd939c64c9 Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Tue, 26 Jul 2022 17:38:31 +0200 Subject: [PATCH 3/3] Add axis labels --- PWGHF/Tasks/HFMCValidation.cxx | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/PWGHF/Tasks/HFMCValidation.cxx b/PWGHF/Tasks/HFMCValidation.cxx index cfd81af0035..59452343eba 100644 --- a/PWGHF/Tasks/HFMCValidation.cxx +++ b/PWGHF/Tasks/HFMCValidation.cxx @@ -348,6 +348,14 @@ struct ValidationRecLevel { } histAmbiguousTracks = registry.add("histAmbiguousTracks", "Tracks that are ambiguous vs. origin;#it{p}_{T}^{reco} (GeV/#it{c});entries", HistType::kTH2F, {{4, -1.5, 2.5}, {50, 0., 10.}}); histTracks = registry.add("histTracks", "Tracks vs. origin;#it{p}_{T}^{reco} (GeV/#it{c});entries", HistType::kTH2F, {{4, -1.5, 2.5}, {50, 0., 10.}}); + histTracks->GetXaxis()->SetBinLabel(1, "no MC particle"); + histTracks->GetXaxis()->SetBinLabel(2, "no quark"); + histTracks->GetXaxis()->SetBinLabel(3, "charm"); + histTracks->GetXaxis()->SetBinLabel(4, "beauty"); + histAmbiguousTracks->GetXaxis()->SetBinLabel(1, "no MC particle"); + histAmbiguousTracks->GetXaxis()->SetBinLabel(2, "no quark"); + histAmbiguousTracks->GetXaxis()->SetBinLabel(3, "charm"); + histAmbiguousTracks->GetXaxis()->SetBinLabel(4, "beauty"); for (auto iHad = 0; iHad < nCharmHadrons; ++iHad) { histDeltaPt[iHad] = registry.add(Form("histDeltaPt%s", particleNames[iHad].data()), Form("Pt difference reco - MC %s; #it{p}_{T}^{reco} - #it{p}_{T}^{gen} (GeV/#it{c}); entries", labels[iHad].data()), HistType::kTH1F, {{2000, -1., 1.}}); histDeltaPx[iHad] = registry.add(Form("histDeltaPx%s", particleNames[iHad].data()), Form("Px difference reco - MC %s; #it{p}_{x}^{reco} - #it{p}_{x}^{gen} (GeV/#it{c}); entries", labels[iHad].data()), HistType::kTH1F, {{2000, -1., 1.}});