Skip to content

PWGHF: add information about tracks with ambiguous collision association #1057

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
fgrosa merged 3 commits into from
Jul 26, 2022
Merged

PWGHF: add information about tracks with ambiguous collision association #1057

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
40a8cfd
PWGHF: add information about tracks with ambiguous collision association
fgrosa Jul 26, 2022
cac151e
Fix comment
fgrosa Jul 26, 2022
b49f78a
Add axis labels
fgrosa Jul 26, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
105 changes: 93 additions & 12 deletions PWGHF/Tasks/HFMCValidation.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,63 @@ static const std::array<std::string, nCharmHadrons> particleNames = {"Dplus", "D
static const std::array<std::string, 2> originNames = {"Prompt", "NonPrompt"};
} // namespace

/// Add relevant information about ambiguous tracks
namespace o2::aod
{
// 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", //!
IsAmbiguousTrack,
AmbiguousCollisionIndicesIds);
} // namespace o2::aod

struct AddAmbiguousTrackInfo {

Produces<o2::aod::TracksWithAmbiguousCollisionInfo> trackWithAmbiguousInfo;
using TracksWithSel = soa::Join<aod::Tracks, aod::TrackSelection>;

void process(TracksWithSel const& tracks,
aod::AmbiguousTracks const& ambitracks,
aod::Collisions const& collisions,
aod::BCs const& bcs)
{
// loop over ambiguous tracks
std::vector<int> trackIndices{};
std::vector<int> ambTrackIndices{};
for (auto& ambitrack : ambitracks) {
auto track = ambitrack.track_as<TracksWithSel>(); // 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<int> 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
Expand All @@ -57,8 +114,8 @@ struct ValidationGenLevel {
Configurable<double> xVertexMax{"xVertexMax", 100., "max. x of generated primary vertex [cm]"};
Configurable<double> yVertexMin{"yVertexMin", -100., "min. y of generated primary vertex [cm]"};
Configurable<double> yVertexMax{"yVertexMax", 100., "max. y of generated primary vertex [cm]"};
Configurable<double> zVertexMin{"zVertexMin", -10., "min. z of generated primary vertex [cm]"};
Configurable<double> zVertexMax{"zVertexMax", 10., "max. z of generated primary vertex [cm]"};
Configurable<double> zVertexMin{"zVertexMin", -100., "min. z of generated primary vertex [cm]"};
Configurable<double> zVertexMax{"zVertexMax", 100., "max. z of generated primary vertex [cm]"};

std::shared_ptr<TH1> hPromptCharmHadronsPtDistr, hPromptCharmHadronsYDistr, hNonPromptCharmHadronsPtDistr, hNonPromptCharmHadronsYDistr, hPromptCharmHadronsDecLenDistr, hNonPromptCharmHadronsDecLenDistr, hQuarkPerEvent;

Expand Down Expand Up @@ -274,20 +331,31 @@ struct ValidationRecLevel {
std::array<std::array<std::array<std::shared_ptr<TH1>, 3>, 2>, nCharmHadrons> histPtDau, histEtaDau, histImpactParameterDau;
std::array<std::array<std::shared_ptr<TH1>, 2>, nCharmHadrons> histPtReco;
std::array<std::shared_ptr<THnSparse>, 4> histOriginTracks;
std::shared_ptr<TH2> histAmbiguousTracks, histTracks;

HistogramRegistry registry{"registry", {}};
void init(o2::framework::InitContext&)
{
histOriginTracks[0] = registry.add<THnSparse>("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<THnSparse>("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<THnSparse>("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<THnSparse>("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<THnSparse>("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<THnSparse>("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<THnSparse>("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<THnSparse>("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<TH2>("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<TH2>("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<TH1>(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<TH1>(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.}});
Expand All @@ -311,7 +379,7 @@ struct ValidationRecLevel {
using HfCandProng2WithMCRec = soa::Join<aod::HfCandProng2, aod::HfCandProng2MCRec>;
using HfCandProng3WithMCRec = soa::Join<aod::HfCandProng3, aod::HfCandProng3MCRec>;
using CollisionsWithMCLabels = soa::Join<aod::Collisions, aod::McCollisionLabels>;
using TracksWithSel = soa::Join<aod::BigTracksMC, aod::TrackSelection>;
using TracksWithSel = soa::Join<aod::BigTracksMC, aod::TrackSelection, aod::TracksWithAmbiguousCollisionInfo>;

void process(HfCandProng2WithMCRec const& cand2Prongs, HfCandProng3WithMCRec const& cand3Prongs, TracksWithSel const& tracks, aod::McParticles const& particlesMC, aod::McCollisions const& mcCollisions, CollisionsWithMCLabels const& collisions)
{
Expand All @@ -324,20 +392,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<CollisionsWithMCLabels>();
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());
}
}

Expand Down Expand Up @@ -475,6 +555,7 @@ struct ValidationRecLevel {
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
{
WorkflowSpec workflow{
adaptAnalysisTask<AddAmbiguousTrackInfo>(cfgc),
adaptAnalysisTask<ValidationGenLevel>(cfgc, TaskName{"hf-mc-validation-gen"}),
adaptAnalysisTask<ValidationRecLevel>(cfgc, TaskName{"hf-mc-validation-rec"})};
return workflow;
Expand Down