Skip to content

[PWGCF] FemtoUniverse: update in the MC functions for D0s #12619

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
prchakra merged 2 commits into from
Aug 18, 2025
Merged

[PWGCF] FemtoUniverse: update in the MC functions for D0s #12619

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
61f6e21
FemtoUniverse: updating the MC checks for D0s
kgwizdzi Aug 18, 2025
0de705b
FU: fixing clang
kgwizdzi Aug 18, 2025
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
12 changes: 10 additions & 2 deletions PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerMCTruthTask.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
#include "PWGCF/FemtoUniverse/DataModel/FemtoDerived.h"

#include "Common/CCDB/TriggerAliases.h"
#include "Common/Core/RecoDecay.h"
#include "Common/DataModel/EventSelection.h"
#include "Common/DataModel/Multiplicity.h"

Expand Down Expand Up @@ -171,10 +172,10 @@ struct FemtoUniverseProducerMCTruthTask {
if (confAnalysisWithPID) {
bool pass = false;
std::vector<int> tmpPDGCodes = confPDGCodes; // necessary due to some features of the Configurable
for (const int& pdg : tmpPDGCodes) {
for (auto const& pdg : tmpPDGCodes) {
if (pdgCode == Pdg::kPhi) { // phi meson
pass = true;
} else if (pdgCode == Pdg::kD0) { // D0 meson
} else if (std::abs(pdgCode) == Pdg::kD0) { // D0(bar) meson
pass = true;
} else if (pdgCode == Pdg::kDPlus) { // D+ meson
pass = true;
Expand All @@ -187,6 +188,13 @@ struct FemtoUniverseProducerMCTruthTask {
continue;
}

/// check if we end-up with the correct final state using MC info
int8_t sign = 0;
if (std::abs(pdgCode) == Pdg::kD0 && !RecoDecay::isMatchedMCGen(tracks, particle, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign)) {
/// check if we have D0(bar) → π± K∓
continue;
}

// we cannot use isSelectedMinimal since it takes Ncls
// if (!trackCuts.isSelectedMinimal(track)) {
// continue;
Expand Down
94 changes: 31 additions & 63 deletions PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@
#include <CCDB/BasicCCDBManager.h>

#include "Math/Vector4D.h"
#include "TLorentzVector.h"

Check failure on line 57 in PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
#include "TMath.h"
#include <TPDGCode.h>

Expand Down Expand Up @@ -328,10 +328,7 @@

// D0/D0bar mesons
struct : o2::framework::ConfigurableGroup {
Configurable<float> confD0D0barCandMaxY{"confD0D0barCandMaxY", -1., "max. cand. rapidity"};
Configurable<float> confD0D0barCandEtaCut{"confD0D0barCandEtaCut", 0.8, "max. cand. pseudorapidity"};
Configurable<float> yD0D0barCandRecoMax{"yD0D0barCandRecoMax", 0.8, "MC Reco, max. rapidity of D0/D0bar cand."};
Configurable<float> yD0D0barCandGenMax{"yD0D0barCandGenMax", 0.8, "MC Truth, max. rapidity of D0/D0bar cand."};
Configurable<float> trackD0pTGenMin{"trackD0pTGenMin", 0.0, "MC Truth, min. pT for tracks and D0/D0bar cand."};
Configurable<float> trackD0pTGenMax{"trackD0pTGenMax", 24.0, "MC Truth, max. pT for tracks and D0/D0bar cand."};
Configurable<bool> storeD0D0barDoubleMassHypo{"storeD0D0barDoubleMassHypo", false, "Store D0/D0bar cand. which pass selection criteria for both, D0 and D0bar"};
Expand Down Expand Up @@ -458,11 +455,11 @@
}

template <typename TrackType>
aod::femtouniverseparticle::CutContainerType PIDBitmask(const TrackType& track)

Check failure on line 458 in PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/function-variable] 

Use lowerCamelCase for names of functions and variables.
{
static const o2::track::PID pids[] = {o2::track::PID::Proton, o2::track::PID::Pion, o2::track::PID::Kaon};
aod::femtouniverseparticle::CutContainerType mask = 0u;
for (UInt_t i = 0; i < 3; ++i) {

Check failure on line 462 in PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.

Check failure on line 462 in PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/entity] 

Replace ROOT entities with equivalents from standard C++ or from O2.
if (isNSigmaTPC(trackCuts.getNsigmaTPC(track, pids[i])))
mask |= (1u << i);
if (isNSigmaTOF(track.p(), trackCuts.getNsigmaTOF(track, pids[i]), track.hasTOF()))
Expand Down Expand Up @@ -882,8 +879,8 @@
phiOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kFake;
}

TLorentzVector part1Vec;

Check failure on line 882 in PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
TLorentzVector part2Vec;

Check failure on line 883 in PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.

const auto mMassOne = o2::constants::physics::MassKPlus; // FIXME: Get from the PDG service of the common header
const auto mMassTwo = o2::constants::physics::MassKMinus; // FIXME: Get from the PDG service of the common header
Expand All @@ -891,7 +888,7 @@
part1Vec.SetPtEtaPhiM(kaon1MC.pt(), kaon1MC.eta(), kaon1MC.phi(), mMassOne);
part2Vec.SetPtEtaPhiM(kaon2MC.pt(), kaon2MC.eta(), kaon2MC.phi(), mMassTwo);

TLorentzVector sumVec(part1Vec);

Check failure on line 891 in PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
sumVec += part2Vec;

float phiEta = sumVec.Eta();
Expand Down Expand Up @@ -1357,10 +1354,6 @@
continue;
}

if (ConfD0Selection.confD0D0barCandMaxY >= 0. && std::abs(hfHelper.yD0(hfCand)) > ConfD0Selection.confD0D0barCandMaxY) {
continue;
}

if (std::abs(hfCand.eta()) > ConfD0Selection.confD0D0barCandEtaCut) {
continue;
}
Expand Down Expand Up @@ -1480,10 +1473,6 @@
continue;
}

if (ConfD0Selection.confD0D0barCandMaxY >= 0. && std::abs(hfHelper.yD0(hfCand)) > ConfD0Selection.confD0D0barCandMaxY) {
continue;
}

if (std::abs(hfCand.eta()) > ConfD0Selection.confD0D0barCandEtaCut) {
continue;
}
Expand Down Expand Up @@ -1610,10 +1599,6 @@
continue;
}

if (ConfD0Selection.confD0D0barCandMaxY >= 0. && std::abs(hfHelper.yD0(hfCand)) > ConfD0Selection.confD0D0barCandMaxY) {
continue;
}

if (std::abs(hfCand.eta()) > ConfD0Selection.confD0D0barCandEtaCut) {
continue;
}
Expand Down Expand Up @@ -1774,8 +1759,8 @@
continue;
}

TLorentzVector part1Vec;

Check failure on line 1762 in PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
TLorentzVector part2Vec;

Check failure on line 1763 in PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.

const auto mMassOne = o2::constants::physics::MassKPlus; // FIXME: Get from the PDG service of the common header
const auto mMassTwo = o2::constants::physics::MassKMinus; // FIXME: Get from the PDG service of the common header
Expand All @@ -1783,7 +1768,7 @@
part1Vec.SetPtEtaPhiM(p1.pt(), p1.eta(), p1.phi(), mMassOne);
part2Vec.SetPtEtaPhiM(p2.pt(), p2.eta(), p2.phi(), mMassTwo);

TLorentzVector sumVec(part1Vec);

Check failure on line 1771 in PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
sumVec += part2Vec;

float phiEta = sumVec.Eta();
Expand Down Expand Up @@ -1997,20 +1982,18 @@
}

template <typename TrackType, bool transientLabels = false>
void fillMCTruthParticlesD0(TrackType const& mcParts, std::optional<std::reference_wrapper<const std::set<int>>> recoMcIds = std::nullopt)
void fillMCTruthParticlesD0(TrackType const& mcParts)
{
std::vector<int> childIDs = {0, 0}; // these IDs are necessary to keep track of the children
std::vector<int> tmpIDtrack;
float ptGenB = -1;

for (const auto& particle : mcParts) {

if (particle.eta() < -ConfFilterCuts.confEtaFilterCut || particle.eta() > ConfFilterCuts.confEtaFilterCut)
continue;
if (particle.pt() < ConfD0Selection.trackD0pTGenMin || particle.pt() > ConfD0Selection.trackD0pTGenMax)
continue;

uint32_t pdgCode = static_cast<uint32_t>(particle.pdgCode());
int pdgCode = particle.pdgCode();

if (ConfGeneral.confMCTruthAnalysisWithPID) {
bool pass = false;
Expand All @@ -2019,10 +2002,10 @@
if (static_cast<int>(pdg) == static_cast<int>(pdgCode)) {
if (pdgCode == Pdg::kPhi) {
pass = true;
} else if (pdgCode == Pdg::kD0) {
} else if (std::abs(pdgCode) == Pdg::kD0) {
pass = true;
} else {
if (particle.isPhysicalPrimary() || (ConfGeneral.confActivateSecondaries && recoMcIds && recoMcIds->get().contains(particle.globalIndex())))
if (particle.isPhysicalPrimary())
pass = true;
}
}
Expand All @@ -2031,41 +2014,31 @@
continue;
}

if (pdgCode == Pdg::kD0) {
if (std::abs(particle.flagMcMatchGen()) == o2::hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiK) {
if (ConfD0Selection.yD0D0barCandGenMax >= 0. && std::abs(RecoDecay::y(particle.pVector(), o2::constants::physics::MassD0)) > ConfD0Selection.yD0D0barCandGenMax) {
continue;
}
if (particle.originMcGen() == RecoDecay::OriginType::Prompt) {
ptGenB = -1;
} else {
ptGenB = mcParts.rawIteratorAt(particle.idxBhadMotherPart()).pt();
}
outputParts(outputCollision.lastIndex(),
particle.pt(),
particle.eta(),
particle.phi(),
aod::femtouniverseparticle::ParticleType::kMCTruthTrack,
-999.,
pdgCode,
pdgCode, // getter tempFitVar
childIDs,
particle.flagMcMatchGen(),
ptGenB); // pT of the B hadron (mother particle, only when non-prompt D0)
}
} else {
outputParts(outputCollision.lastIndex(),
particle.pt(),
particle.eta(),
particle.phi(),
aod::femtouniverseparticle::ParticleType::kMCTruthTrack,
-999.,
pdgCode,
pdgCode,
childIDs,
-999.,
-999.);
/// check if we end-up with the correct final state using MC info
int8_t sign = 0;
int8_t origin = -99;
int8_t mcGenFlag = -99;
if (std::abs(particle.pdgCode()) == Pdg::kD0 && !RecoDecay::isMatchedMCGen(mcParts, particle, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign)) {
/// check if we have D0(bar) → π± K∓
continue;
}
if (std::abs(particle.pdgCode()) == Pdg::kD0) {
origin = RecoDecay::getCharmHadronOrigin(mcParts, particle);
mcGenFlag = particle.flagMcMatchGen();
}

outputParts(outputCollision.lastIndex(),
particle.pt(),
particle.eta(),
particle.phi(),
aod::femtouniverseparticle::ParticleType::kMCTruthTrack,
-999.,
pdgCode,
pdgCode, // getter tempFitVar
childIDs,
mcGenFlag,
origin); // D0(bar) origin

if (confIsDebug) {
fillDebugParticle<false, true, false>(particle);
}
Expand All @@ -2074,7 +2047,7 @@
// aligned, so that they can be joined in the task.
if constexpr (transientLabels) {
outputPartsMCLabels(-1);
outputDebugPartsMC(9999);
outputDebugPartsMC(-999);
}
}
}
Expand Down Expand Up @@ -2534,7 +2507,6 @@
aod::McParticles const& mcParts)
{
// MC Reco
std::set<int> recoMcIds;
for (const auto& col : collisions) {
auto groupedTracks = tracks.sliceBy(perCollisionTracks, col.globalIndex());
auto groupedD0s = hfMcRecoCands.sliceBy(perCollisionD0s, col.globalIndex());
Expand All @@ -2545,18 +2517,14 @@
if (colcheck) {
fillTracks<true>(groupedTracks);
fillD0D0barMcMl<true>(col, groupedTracks, groupedD0s, mcParts);
for (const auto& track : groupedTracks) {
if (trackCuts.isSelectedMinimal(track))
recoMcIds.insert(track.mcParticleId());
}
}
}
// MC Truth
for (const auto& mccol : mccols) {
auto groupedMCParticles = hfMcGenCands.sliceBy(mcPartPerMcColl, mccol.globalIndex());
auto groupedCollisions = collisions.sliceBy(recoCollsPerMCColl, mccol.globalIndex());
fillMCTruthCollisions(groupedCollisions, groupedMCParticles); // fills the reco collisions for mc collision
fillMCTruthParticlesD0<decltype(groupedMCParticles), true>(groupedMCParticles, recoMcIds); // fills mc particles
fillMCTruthCollisions(groupedCollisions, groupedMCParticles); // fills the reco collisions for mc collision
fillMCTruthParticlesD0<decltype(groupedMCParticles), true>(groupedMCParticles); // fills mc particles
}
}
PROCESS_SWITCH(FemtoUniverseProducerTask, processTrackD0MC, "Provide MC data for track D0 analysis", false);
Expand Down
8 changes: 4 additions & 4 deletions PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackD0.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1252,9 +1252,9 @@ struct FemtoUniversePairTaskTrackD0 {
if (pdgCode == o2::constants::physics::Pdg::kD0) {
if (std::abs(hfFlagMcGen) == o2::hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiK) {
mcTruthRegistry.fill(HIST("hMcGenD0"), part.pt(), part.eta());
if (part.mAntiLambda() > 0) {
if (part.mAntiLambda() == 1) {
mcTruthRegistry.fill(HIST("hMcGenD0Prompt"), part.pt(), part.eta());
} else {
} else if (part.mAntiLambda() == 2) {
mcTruthRegistry.fill(HIST("hMcGenD0NonPrompt"), part.pt(), part.eta());
}
}
Expand All @@ -1278,9 +1278,9 @@ struct FemtoUniversePairTaskTrackD0 {
if (pdgCode == o2::constants::physics::Pdg::kD0Bar) {
if (std::abs(hfFlagMcGen) == o2::hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiK) {
mcTruthRegistry.fill(HIST("hMcGenD0bar"), part.pt(), part.eta());
if (part.mAntiLambda() > 0) {
if (part.mAntiLambda() == 1) {
mcTruthRegistry.fill(HIST("hMcGenD0barPrompt"), part.pt(), part.eta());
} else {
} else if (part.mAntiLambda() == 2) {
mcTruthRegistry.fill(HIST("hMcGenD0barNonPrompt"), part.pt(), part.eta());
}
}
Expand Down
Loading
Loading