Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
22 changes: 17 additions & 5 deletions ALICE3/TableProducer/alice3-decayfinder.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -151,10 +151,10 @@

HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};

Partition<aod::McParticles> trueD = aod::mcparticle::pdgCode == 421;

Check failure on line 154 in ALICE3/TableProducer/alice3-decayfinder.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code]

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
Partition<aod::McParticles> trueDbar = aod::mcparticle::pdgCode == -421;

Check failure on line 155 in ALICE3/TableProducer/alice3-decayfinder.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code]

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
Partition<aod::McParticles> trueLc = aod::mcparticle::pdgCode == 4122;

Check failure on line 156 in ALICE3/TableProducer/alice3-decayfinder.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code]

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
Partition<aod::McParticles> trueLcbar = aod::mcparticle::pdgCode == -4122;

Check failure on line 157 in ALICE3/TableProducer/alice3-decayfinder.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code]

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.

// filter expressions for D mesons
static constexpr uint32_t trackSelectionPiPlusFromD = 1 << kInnerTOFPion | 1 << kOuterTOFPion | 1 << kRICHPion | 1 << kTruePiPlusFromD;
Expand Down Expand Up @@ -284,12 +284,12 @@
dmeson.Ndaug[2] = negP[2];

// return mass and kinematic variables
dmeson.mass = RecoDecay::m(array{array{posP[0], posP[1], posP[2]}, array{negP[0], negP[1], negP[2]}}, array{posMass, negMass});

Check failure on line 287 in ALICE3/TableProducer/alice3-decayfinder.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix]

Use std:: prefix for names from the std namespace.
dmeson.pt = std::hypot(posP[0] + negP[0], posP[1] + negP[1]);
dmeson.ptdaugPos = std::hypot(posP[0], posP[1]);
dmeson.ptdaugNeg = std::hypot(negP[0], negP[1]);
dmeson.phi = RecoDecay::phi(array{posP[0] + negP[0], posP[1] + negP[1]});

Check failure on line 291 in ALICE3/TableProducer/alice3-decayfinder.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix]

Use std:: prefix for names from the std namespace.
dmeson.eta = RecoDecay::eta(array{posP[0] + negP[0], posP[1] + negP[1], posP[2] + negP[2]});

Check failure on line 292 in ALICE3/TableProducer/alice3-decayfinder.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix]

Use std:: prefix for names from the std namespace.
dmeson.y = RecoDecay::y(std::array{posP[0] + negP[0], posP[1] + negP[1], posP[2] + negP[2]}, dmeson.mass);
const auto posSV = fitter2prongs.getPCACandidate();
dmeson.posSV[0] = posSV[0];
Expand Down Expand Up @@ -406,14 +406,14 @@
mCandidate3Prong.errorImpactParameterZ2 = impactParameter2.getSigmaZ2();

// return mass
mCandidate3Prong.mass = RecoDecay::m(array{array{P0[0], P0[1], P0[2]},

Check failure on line 409 in ALICE3/TableProducer/alice3-decayfinder.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix]

Use std:: prefix for names from the std namespace.
array{P1[0], P1[1], P1[2]},
array{P2[0], P2[1], P2[2]}},
daughtersMasses3Prong);

mCandidate3Prong.pt = std::hypot(P0[0] + P1[0] + P2[0], P0[1] + P1[1] + P2[1]);
mCandidate3Prong.phi = RecoDecay::phi(array{P0[0] + P1[0] + P2[0], P0[1] + P1[1] + P2[1]});

Check failure on line 415 in ALICE3/TableProducer/alice3-decayfinder.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix]

Use std:: prefix for names from the std namespace.
mCandidate3Prong.eta = RecoDecay::eta(array{P0[0] + P1[0] + P2[0], P0[1] + P1[1] + P2[1], P0[2] + P1[2] + P2[2]});

Check failure on line 416 in ALICE3/TableProducer/alice3-decayfinder.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix]

Use std:: prefix for names from the std namespace.
mCandidate3Prong.Pdaug0[0] = P0[0];
mCandidate3Prong.Pdaug0[1] = P0[1];
mCandidate3Prong.Pdaug0[2] = P0[2];
Expand Down Expand Up @@ -661,7 +661,7 @@
if (doprocessFindLc) {
for (auto const& mcParticle : mcParticles) {
if (std::abs(mcParticle.pdgCode()) != motherPdgCode) {
mcGenFlags(-1, -1, -1);
mcGenFlags(-1, 0, 0);
continue;
}
std::vector<int> idxBhadMothers{};
Expand All @@ -671,7 +671,7 @@
auto bHadMother = mcParticles.rawIteratorAt(idxBhadMothers[0]);
ptBMotherGen = bHadMother.pt();
}
mcGenFlags(origin, ptBMotherGen, mcParticle.pdgCode() ? charmHadFlag : -charmHadFlag);
mcGenFlags(origin, ptBMotherGen, mcParticle.pdgCode() > 0 ? charmHadFlag : -charmHadFlag);
if (mcParticle.pdgCode() > 0) {
histos.fill(HIST("h2dGen3Prong"), mcParticle.pt(), mcParticle.eta());
} else {
Expand Down Expand Up @@ -1040,9 +1040,9 @@
mCandidate3Prong.eta,
mCandidate3Prong.phi,
mCandidate3Prong.pt,
mCandidate3Prong.Pdaug2[0], mCandidate3Prong.Pdaug2[1], mCandidate3Prong.Pdaug2[2],
mCandidate3Prong.Pdaug1[0], mCandidate3Prong.Pdaug1[1], mCandidate3Prong.Pdaug1[2],
mCandidate3Prong.Pdaug0[0], mCandidate3Prong.Pdaug0[1], mCandidate3Prong.Pdaug0[2],
mCandidate3Prong.Pdaug1[0], mCandidate3Prong.Pdaug1[1], mCandidate3Prong.Pdaug1[2],
mCandidate3Prong.Pdaug2[0], mCandidate3Prong.Pdaug2[1], mCandidate3Prong.Pdaug2[2],
mCandidate3Prong.impactParameterY0, mCandidate3Prong.impactParameterY1, mCandidate3Prong.impactParameterY2,
std::sqrt(mCandidate3Prong.errorImpactParameterY0),
std::sqrt(mCandidate3Prong.errorImpactParameterY1),
Expand All @@ -1066,15 +1066,27 @@

void processFindLc(aod::Collision const& collision,
aod::McParticles const& mcParticles,
Alice3TracksWPid const&)
Alice3TracksWPid const& tracks)
{
LOG(debug) << "Processing Lc candidates for collision " << collision.globalIndex() << " with " << tracks.size() << " tracks";
for (auto const& track : tracks) {
if (track.has_mcParticle()) {
LOG(debug) << " - track index: " << track.globalIndex() << ", pT: " << track.pt() << " (MC pt " << track.mcParticle().pt() << "), PID: " << track.mcParticle().pdgCode();
}
}

auto tracksPiPlus = tracksPiPlusFromLc->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
LOG(debug) << " - found " << tracksPiPlus.size() << " pi+ from Lc";
auto tracksKaPlus = tracksKaPlusFromLc->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
LOG(debug) << " - found " << tracksKaPlus.size() << " K+ from Lc";
auto tracksPrPlus = tracksPrPlusFromLc->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
LOG(debug) << " - found " << tracksPrPlus.size() << " p+ from Lc";
auto tracksPiMinus = tracksPiMinusFromLc->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
LOG(debug) << " - found " << tracksPiMinus.size() << " pi- from Lc";
auto tracksKaMinus = tracksKaMinusFromLc->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
LOG(debug) << " - found " << tracksKaMinus.size() << " K- from Lc";
auto tracksPrMinus = tracksPrMinusFromLc->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
LOG(debug) << " - found " << tracksPrMinus.size() << " p- from Lc";

if (doDCAplots3Prong) {
for (auto const& track : tracksPiPlus)
Expand Down
20 changes: 9 additions & 11 deletions ALICE3/TableProducer/alice3HfSelector3Prong.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -313,8 +313,6 @@ struct Alice3HfSelector3Prong {
template <CharmHadAlice3 CharmHad, typename CandType>
void runSelect3Prong(CandType const& cands)
{
bool isSelMassHypo0{false};
bool isSelMassHypo1{false};
std::vector<float> outputMl{-1.f, -1.f, -1.f};
uint32_t pidMask = 0;

Expand All @@ -324,21 +322,21 @@ struct Alice3HfSelector3Prong {
outputMl = {-1.f, -1.f, -1.f};
pidMask = 0;

auto ptCand = cand.pt();
int const ptBin = findBin(binsPt, ptCand);
const float ptCand = cand.pt();
const int ptBin = findBin(binsPt, ptCand);
if (ptBin == -1) {
candSelFlags(isSelMassHypo0, isSelMassHypo1, pidMask);
candSelFlags(false, false, pidMask);
if (applyMl) {
candMlScores(outputMl[0], outputMl[1], outputMl[2]);
}
continue;
}

// Here all cands pass the cut on the mass selection
bool const selMassHypo0 = selectionCandidateMass<CharmHad, false>(ptBin, cand);
bool const selMassHypo1 = selectionCandidateMass<CharmHad, true>(ptBin, cand);
const bool selMassHypo0 = selectionCandidateMass<CharmHad, false>(ptBin, cand);
const bool selMassHypo1 = selectionCandidateMass<CharmHad, true>(ptBin, cand);
if (!selMassHypo0 && !selMassHypo1) {
candSelFlags(isSelMassHypo0, isSelMassHypo1, pidMask);
candSelFlags(false, false, pidMask);
if (applyMl) {
candMlScores(outputMl[0], outputMl[1], outputMl[2]);
}
Expand All @@ -348,7 +346,7 @@ struct Alice3HfSelector3Prong {

// Topological selection (TODO: track quality selection)
if (!selectionTopol(cand, ptCand)) {
candSelFlags(isSelMassHypo0, isSelMassHypo1, pidMask);
candSelFlags(false, false, pidMask);
if (applyMl) {
candMlScores(outputMl[0], outputMl[1], outputMl[2]);
}
Expand All @@ -359,7 +357,7 @@ struct Alice3HfSelector3Prong {
// PID selection
configurePidMask<CharmHad>(cand, pidMask);
if (pidMask == 0) {
candSelFlags(isSelMassHypo0, isSelMassHypo1, pidMask);
candSelFlags(false, false, pidMask);
if (applyMl) {
candMlScores(outputMl[0], outputMl[1], outputMl[2]);
}
Expand All @@ -375,7 +373,7 @@ struct Alice3HfSelector3Prong {
isSelectedMl = mlResponse.isSelectedMl(inputFeaturesMassHypo0, ptCand, outputMl);
candMlScores(outputMl[0], outputMl[1], outputMl[2]);
if (!isSelectedMl) {
candSelFlags(isSelMassHypo0, isSelMassHypo1, pidMask);
candSelFlags(false, false, pidMask);
continue;
}

Expand Down
80 changes: 58 additions & 22 deletions ALICE3/Tasks/alice3HfTask3Prong.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -87,14 +87,15 @@ struct Alice3HfTask3Prong {
HistogramRegistry registry{"registry", {}};

// Names of folders and suffixes for MC signal histograms
constexpr static std::string_view SignalFolders[] = {"signal", "prompt", "nonprompt"};
constexpr static std::string_view SignalSuffixes[] = {"", "Prompt", "NonPrompt"};

enum SignalClasses : int {
Signal = 0,
Prompt,
NonPrompt
NonPrompt,
Bkg,
NSignalClasses
};
constexpr static std::string_view SignalFolders[SignalClasses::NSignalClasses] = {"signal", "prompt", "nonprompt", "background"};
constexpr static std::string_view SignalSuffixes[SignalClasses::NSignalClasses] = {"", "Prompt", "NonPrompt", "Bkg"};

void init(InitContext&)
{
Expand All @@ -109,9 +110,11 @@ struct Alice3HfTask3Prong {
}

auto addHistogramsRec = [&](const std::string& histoName, const std::string& xAxisTitle, const std::string& yAxisTitle, const HistogramConfigSpec& configSpec) {
registry.add(("MC/rec/signal/" + histoName + "RecSig").c_str(), ("3-prong cands (matched);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
registry.add(("MC/rec/prompt/" + histoName + "RecSigPrompt").c_str(), ("3-prong cands (matched, prompt);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
registry.add(("MC/rec/nonprompt/" + histoName + "RecSigNonPrompt").c_str(), ("3-prong cands (matched, non-prompt);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
const char* basePath = "MC/rec";
registry.add(Form("%s/signal/%sRecSig%s", basePath, histoName.c_str(), SignalSuffixes[SignalClasses::Signal].data()), ("3-prong cands (matched);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
registry.add(Form("%s/prompt/%sRecSig%s", basePath, histoName.c_str(), SignalSuffixes[SignalClasses::Prompt].data()), ("3-prong cands (matched, prompt);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
registry.add(Form("%s/nonprompt/%sRecSig%s", basePath, histoName.c_str(), SignalSuffixes[SignalClasses::NonPrompt].data()), ("3-prong cands (matched, non-prompt);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
registry.add(Form("%s/background/%sRecSig%s", basePath, histoName.c_str(), SignalSuffixes[SignalClasses::Bkg].data()), ("3-prong cands (unmatched);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
};

auto addHistogramsGen = [&](const std::string& histoName, const std::string& xAxisTitle, const std::string& yAxisTitle, const HistogramConfigSpec& configSpec) {
Expand Down Expand Up @@ -174,9 +177,16 @@ struct Alice3HfTask3Prong {
auto h2 = registry.add<TH2>("hSelectionStatus", "3-prong cands;selection status;entries", {HistType::kTH2F, {{5, -0.5, 4.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}});
h2->GetXaxis()->SetBinLabel(1, "mass hypo 0");
h2->GetXaxis()->SetBinLabel(2, "mass hypo 1");
auto h = registry.add<TH1>("MC/rec/hCandidateCounter", "Candidate counter;entries", {HistType::kTH1D, {{2, -0.5, 1.5}}});
auto h = registry.add<TH1>("MC/rec/hCandidateCounter", "Candidate counter;entries", {HistType::kTH1D, {{4, -0.5, 3.5}}});
h->GetXaxis()->SetBinLabel(1, "Calls");
h->GetXaxis()->SetBinLabel(2, "Candidates");
h->GetXaxis()->SetBinLabel(3, "Passed Y cut");
h->GetXaxis()->SetBinLabel(4, "Has MC match");

registry.add("MC/rec/hPtDeltaProng0", ";prong 0 (#it{p}_{T}-#it{p}_{T, gen}) (GeV/#it{c});entries", {HistType::kTH1F, {{100, -5, 5.}}});
registry.add("MC/rec/hPtDeltaProng1", ";prong 1 (#it{p}_{T}-#it{p}_{T, gen}) (GeV/#it{c});entries", {HistType::kTH1F, {{100, -5, 5.}}});
registry.add("MC/rec/hPtDeltaProng2", ";prong 2 (#it{p}_{T}-#it{p}_{T, gen}) (GeV/#it{c});entries", {HistType::kTH1F, {{100, -5, 5.}}});

// Number of events processed
h = registry.add<TH1>("hNEventsProcessed", "number of events processed;entries;", {HistType::kTH1F, {{2, 0.5, 2.5}}});
h->GetXaxis()->SetBinLabel(1, "Generated");
Expand Down Expand Up @@ -207,7 +217,7 @@ struct Alice3HfTask3Prong {
/// Helper function for filling MC reconstructed histograms for prompt, nonpromt and common (signal)
/// \param candidate is a reconstructed candidate
/// \tparam SignalType is an enum defining which histogram in which folder (signal, prompt or nonpromt) to fill
template <CharmHadAlice3 CharmHad, int SignalType, typename CandidateType>
template <CharmHadAlice3 CharmHad, SignalClasses SignalType, typename CandidateType>
void fillHistogramsRecSig(CandidateType const& candidate, float mass, bool isSwapped = false)
{
static constexpr auto histoPrefix = HIST("MC/rec/") + HIST(SignalFolders[SignalType]) + HIST("/");
Expand Down Expand Up @@ -268,33 +278,50 @@ struct Alice3HfTask3Prong {
if (yCandRecoMax >= 0. && std::abs(hfHelper.getCandY<CharmHad>(candidate)) > yCandRecoMax) {
continue;
}
auto mcParticle = allParticles.iteratorAt(candidate.particleMcRec());
if (candidate.particleMcRec() > 0) {
registry.fill(HIST("MC/rec/hCandidateCounter"), 2.);
if (candidate.particleMcRec() >= 0) {
registry.fill(HIST("MC/rec/hCandidateCounter"), 3.);
auto mcParticle = allParticles.iteratorAt(candidate.particleMcRec());
if (mcParticle.has_daughters()) {
auto daughters = mcParticle.daughtersIds();
LOG(info) << "Reco candidate matched to MC particle with PDG " << mcParticle.pdgCode() << " daughters: " << daughters.size();
for (auto dauId : daughters) {
LOG(debug) << "Reco candidate matched to MC particle with PDG " << mcParticle.pdgCode() << " daughters: " << daughters.size();
int prongIdx = 0;
for (int dauId = daughters[0]; dauId <= daughters[1]; dauId++) {
auto dau = allParticles.iteratorAt(dauId);
LOG(info) << " dauId: " << dauId << " PDG: " << dau.pdgCode();
LOG(debug) << " dauId: " << dauId << " PDG: " << dau.pdgCode() << " with pT: " << dau.pt();
switch (prongIdx) {
case 0:
registry.fill(HIST("MC/rec/hPtDeltaProng0"), candidate.ptProng0() - dau.pt());
break;
case 1:
registry.fill(HIST("MC/rec/hPtDeltaProng1"), candidate.ptProng1() - dau.pt());
break;
case 2:
registry.fill(HIST("MC/rec/hPtDeltaProng2"), candidate.ptProng2() - dau.pt());
break;
default:
break;
}
prongIdx++;
}
}
}

if (candidate.flagMcRec() != 0) {
// Get the corresponding MC particle.
if (candidate.flagMcRec() != 0) { // Particle is matched to MC truth

// Get the corresponding MC particle.
const auto pt = candidate.pt();
const auto originType = candidate.originMcRec();

if (candidate.isSelMassHypo0()) {
registry.fill(HIST("hSelectionStatus"), 0., pt);
double mass = hfHelper.getCandMass<CharmHad, false>(candidate);
/// Fill histograms
fillHistogramsRecSig<CharmHad, Signal>(candidate, mass, false);
fillHistogramsRecSig<CharmHad, SignalClasses::Signal>(candidate, mass, false);
if (originType == RecoDecay::OriginType::Prompt) {
fillHistogramsRecSig<CharmHad, Prompt>(candidate, mass, false);
fillHistogramsRecSig<CharmHad, SignalClasses::Prompt>(candidate, mass, false);
} else if (originType == RecoDecay::OriginType::NonPrompt) {
fillHistogramsRecSig<CharmHad, NonPrompt>(candidate, mass, false);
fillHistogramsRecSig<CharmHad, SignalClasses::NonPrompt>(candidate, mass, false);
}
if (fillThn) {
std::vector<double> valuesToFill{mass, pt};
Expand All @@ -312,11 +339,11 @@ struct Alice3HfTask3Prong {
registry.fill(HIST("hSelectionStatus"), 1., pt);
double mass = hfHelper.getCandMass<CharmHad, true>(candidate);
/// Fill histograms
fillHistogramsRecSig<CharmHad, Signal>(candidate, mass, true);
fillHistogramsRecSig<CharmHad, SignalClasses::Signal>(candidate, mass, true);
if (originType == RecoDecay::OriginType::Prompt) {
fillHistogramsRecSig<CharmHad, Prompt>(candidate, mass, true);
fillHistogramsRecSig<CharmHad, SignalClasses::Prompt>(candidate, mass, true);
} else if (originType == RecoDecay::OriginType::NonPrompt) {
fillHistogramsRecSig<CharmHad, NonPrompt>(candidate, mass, true);
fillHistogramsRecSig<CharmHad, SignalClasses::NonPrompt>(candidate, mass, true);
}
if (fillThn) {
std::vector<double> valuesToFill{mass, pt};
Expand All @@ -330,6 +357,15 @@ struct Alice3HfTask3Prong {
registry.get<THnSparse>(HIST("hSparseRec"))->Fill(valuesToFill.data());
}
}
} else { // Background
if (candidate.isSelMassHypo0()) {
double mass = hfHelper.getCandMass<CharmHad, false>(candidate);
fillHistogramsRecSig<CharmHad, SignalClasses::Bkg>(candidate, mass, false);
}
if (candidate.isSelMassHypo1()) {
double mass = hfHelper.getCandMass<CharmHad, true>(candidate);
fillHistogramsRecSig<CharmHad, SignalClasses::Bkg>(candidate, mass, true);
}
}
}
}
Expand Down
Loading