Skip to content
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.

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

[PWGHF] Add MC Generated ThnSparse in Lc task #8887

Merged
merged 2 commits into from
Dec 9, 2024
Merged
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
61 changes: 46 additions & 15 deletions PWGHF/D2H/Tasks/taskLc.cxx
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.

Check warning on line 1 in PWGHF/D2H/Tasks/taskLc.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pwghf/struct-member-order]

Declare struct members in the conventional order. See the PWGHF coding guidelines.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
Expand Down Expand Up @@ -43,7 +43,7 @@
Configurable<std::vector<double>> binsPt{"binsPt", std::vector<double>{hf_cuts_lc_to_p_k_pi::vecBinsPt}, "pT bin limits"};
// ThnSparse for ML outputScores and Vars
Configurable<bool> fillTHn{"fillTHn", false, "fill THn"};
ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {72, 0, 36}, ""};

Check warning on line 46 in PWGHF/D2H/Tasks/taskLc.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pwghf/struct-member-order]

HfTaskLc: ConfigurableAxis appears too early (before end of PresliceUnsorted<).
ConfigurableAxis thnConfigAxisMass{"thnConfigAxisMass", {300, 1.98, 2.58}, ""};
ConfigurableAxis thnConfigAxisPtProng{"thnConfigAxisPtProng", {100, 0, 20}, ""};
ConfigurableAxis thnConfigAxisCentrality{"thnConfigAxisCentrality", {100, 0, 100}, ""};
Expand All @@ -53,10 +53,13 @@
ConfigurableAxis thnConfigAxisBdtScoreBkg{"thnConfigAxisBdtScoreBkg", {1000, 0., 1.}, ""};
ConfigurableAxis thnConfigAxisBdtScoreSignal{"thnConfigAxisBdtScoreSignal", {100, 0., 1.}, ""};
ConfigurableAxis thnConfigAxisCanType{"thnConfigAxisCanType", {5, 0., 5.}, ""};
ConfigurableAxis thnAxisRapidity{"thnAxisRapidity", {20, -1, 1}, "Cand. rapidity bins"};
ConfigurableAxis thnConfigAxisGenPtB{"thnConfigAxisGenPtB", {1000, 0, 100}, "Gen Pt B"};
ConfigurableAxis thnConfigAxisNumPvContr{"thnConfigAxisNumPvContr", {200, -0.5, 199.5}, "Number of PV contributors"};

HfHelper hfHelper;

using Collisions = soa::Join<aod::Collisions, aod::EvSels>;

Check warning on line 62 in PWGHF/D2H/Tasks/taskLc.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pwghf/struct-member-order]

HfTaskLc: using appears too early (before end of SliceCache).
using CollisionsMc = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels>;
using CollisionsWithFT0C = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs>;
using CollisionsMcWithFT0C = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::CentFT0Cs>;
Expand All @@ -70,7 +73,7 @@
using LcCandidatesMlMc = soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelLc, aod::HfMlLcToPKPi, aod::HfCand3ProngMcRec>>;
using McParticles3ProngMatched = soa::Join<aod::McParticles, aod::HfCand3ProngMcGen>;
Filter filterSelectCandidates = aod::hf_sel_candidate_lc::isSelLcToPKPi >= selectionFlagLc || aod::hf_sel_candidate_lc::isSelLcToPiKP >= selectionFlagLc;
Preslice<aod::McParticles> perMcCollision = aod::mcparticle::mcCollisionId;
PresliceUnsorted<aod::McCollisionLabels> colPerMcCollision = aod::mcparticle::mcCollisionId;

Check warning on line 76 in PWGHF/D2H/Tasks/taskLc.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pwghf/struct-member-order]

HfTaskLc: PresliceUnsorted< appears too early (before end of Preslice<).
Preslice<aod::HfCand3Prong> candLcPerCollision = aod::hf_cand::collisionId;
SliceCache cache;

Expand Down Expand Up @@ -286,11 +289,25 @@
const AxisSpec thnAxisBdtScoreLcPrompt{thnConfigAxisBdtScoreSignal, "BDT prompt score (Lc)"};
const AxisSpec thnAxisBdtScoreLcNonPrompt{thnConfigAxisBdtScoreSignal, "BDT non-prompt score (Lc)"};
const AxisSpec thnAxisCanType{thnConfigAxisCanType, "candidates type"};
const AxisSpec thnAxisY{thnAxisRapidity, "rapidity"};
const AxisSpec thnAxisPtB{thnConfigAxisGenPtB, "#it{p}_{T}^{B} (GeV/#it{c})"};
const AxisSpec thnAxisTracklets{thnConfigAxisNumPvContr, "Number of PV contributors"};

if (doprocessDataWithMl || doprocessDataWithMlWithFT0C || doprocessDataWithMlWithFT0M) {

registry.add("hnLcVarsWithBdt", "THn for Lambdac candidates with BDT scores for data with ML", HistType::kTHnSparseF, {thnAxisMass, thnAxisPt, thnAxisCentrality, thnAxisBdtScoreLcBkg, thnAxisBdtScoreLcPrompt, thnAxisBdtScoreLcNonPrompt, thnAxisTracklets});

} else if (doprocessMcWithMl || doprocessMcWithMlWithFT0C || doprocessMcWithMlWithFT0M) {

registry.add("hnLcVarsWithBdt", "THn for Lambdac candidates with BDT scores for mc with ML", HistType::kTHnSparseF, {thnAxisMass, thnAxisPt, thnAxisCentrality, thnAxisBdtScoreLcBkg, thnAxisBdtScoreLcPrompt, thnAxisBdtScoreLcNonPrompt, thnAxisTracklets, thnAxisPtB, thnAxisCanType});
registry.add("hnLcVarsGen", "THn for Generated Lambdac", HistType::kTHnSparseF, {thnAxisPt, thnAxisY, thnAxisTracklets, thnAxisPtB, thnAxisCanType});

} else if (doprocessDataStd || doprocessDataStdWithFT0C || doprocessDataStdWithFT0M) {
registry.add("hnLcVars", "THn for Reconstructed Lambdac candidates for data without ML", HistType::kTHnSparseF, {thnAxisMass, thnAxisPt, thnAxisCentrality, thnAxisPtProng0, thnAxisPtProng1, thnAxisPtProng2, thnAxisChi2PCA, thnAxisDecLength, thnAxisCPA, thnAxisTracklets});

if (doprocessDataWithMl || doprocessDataWithMlWithFT0C || doprocessDataWithMlWithFT0M || doprocessMcWithMl || doprocessMcWithMlWithFT0C || doprocessMcWithMlWithFT0M) {
registry.add("hnLcVarsWithBdt", "THn for Lambdac candidates with BDT scores", HistType::kTHnSparseF, {thnAxisMass, thnAxisPt, thnAxisCentrality, thnAxisBdtScoreLcBkg, thnAxisBdtScoreLcPrompt, thnAxisBdtScoreLcNonPrompt, thnAxisCanType});
} else {
registry.add("hnLcVars", "THn for Lambdac candidates", HistType::kTHnSparseF, {thnAxisMass, thnAxisPt, thnAxisCentrality, thnAxisPtProng0, thnAxisPtProng1, thnAxisPtProng2, thnAxisChi2PCA, thnAxisDecLength, thnAxisCPA, thnAxisCanType});
registry.add("hnLcVars", "THn for Reconstructed Lambdac candidates for mc without ML", HistType::kTHnSparseF, {thnAxisMass, thnAxisPt, thnAxisCentrality, thnAxisPtProng0, thnAxisPtProng1, thnAxisPtProng2, thnAxisChi2PCA, thnAxisDecLength, thnAxisCPA, thnAxisTracklets, thnAxisPtB, thnAxisCanType});
registry.add("hnLcVarsGen", "THn for Generated Lambdac", HistType::kTHnSparseF, {thnAxisPt, thnAxisY, thnAxisTracklets, thnAxisPtB, thnAxisCanType});
}
}
}
Expand Down Expand Up @@ -341,6 +358,9 @@
auto cpa = candidate.cpa();
auto cpaXY = candidate.cpaXY();
auto originType = candidate.originMcRec();
auto numPvContributors = collision.numContrib();
auto ptRecB = candidate.ptBhadMotherPart();

/// MC reconstructed signal
if ((candidate.isSelLcToPKPi() >= selectionFlagLc) && pdgCodeProng0 == kProton) {
registry.fill(HIST("MC/reconstructed/signal/hMassRecSig"), hfHelper.invMassLcToPKPi(candidate));
Expand Down Expand Up @@ -476,9 +496,9 @@
outputFD = candidate.mlProbLcToPKPi()[2]; /// non-prompt score
}
/// Fill the ML outputScores and variables of candidate
registry.get<THnSparse>(HIST("hnLcVarsWithBdt"))->Fill(massLc, pt, cent, outputBkg, outputPrompt, outputFD, originType);
registry.get<THnSparse>(HIST("hnLcVarsWithBdt"))->Fill(massLc, pt, cent, outputBkg, outputPrompt, outputFD, numPvContributors, ptRecB, originType);
} else {
registry.get<THnSparse>(HIST("hnLcVars"))->Fill(massLc, pt, cent, ptProng0, ptProng1, ptProng2, chi2PCA, decayLength, cpa, originType);
registry.get<THnSparse>(HIST("hnLcVars"))->Fill(massLc, pt, cent, ptProng0, ptProng1, ptProng2, chi2PCA, decayLength, cpa, numPvContributors, ptRecB, originType);
}
}
if ((candidate.isSelLcToPiKP() >= selectionFlagLc) && pdgCodeProng0 == kPiPlus) {
Expand All @@ -491,9 +511,9 @@
outputFD = candidate.mlProbLcToPiKP()[2]; /// non-prompt score
}
/// Fill the ML outputScores and variables of candidate (todo: add multiplicity)
registry.get<THnSparse>(HIST("hnLcVarsWithBdt"))->Fill(massLc, pt, cent, outputBkg, outputPrompt, outputFD, originType);
registry.get<THnSparse>(HIST("hnLcVarsWithBdt"))->Fill(massLc, pt, cent, outputBkg, outputPrompt, outputFD, numPvContributors, ptRecB, originType);
} else {
registry.get<THnSparse>(HIST("hnLcVars"))->Fill(massLc, pt, cent, ptProng0, ptProng1, ptProng2, chi2PCA, decayLength, cpa, originType);
registry.get<THnSparse>(HIST("hnLcVars"))->Fill(massLc, pt, cent, ptProng0, ptProng1, ptProng2, chi2PCA, decayLength, cpa, numPvContributors, ptRecB, originType);
}
}
}
Expand All @@ -503,8 +523,8 @@

/// Fill MC histograms at generated level
/// \tparam fillMl switch to fill ML histograms
template <typename CandLcMcGen>
void fillHistosMcGen(CandLcMcGen const& mcParticles)
template <typename CandLcMcGen, typename Coll>
void fillHistosMcGen(CandLcMcGen const& mcParticles, Coll const& recoCollisions)
{
// MC gen.
for (const auto& particle : mcParticles) {
Expand All @@ -514,6 +534,14 @@
continue;
}
auto ptGen = particle.pt();
auto originType = particle.originMcGen();
auto ptGenB = -1;
unsigned int numPvContributors = 0;
const auto& recoCollsPerMcColl = recoCollisions.sliceBy(colPerMcCollision, particle.mcCollision().globalIndex());
for (const auto& recCol : recoCollsPerMcColl) {
numPvContributors = recCol.numContrib() > numPvContributors ? recCol.numContrib() : numPvContributors;
}

registry.fill(HIST("MC/generated/signal/hPtGen"), ptGen);
registry.fill(HIST("MC/generated/signal/hEtaGen"), particle.eta());
registry.fill(HIST("MC/generated/signal/hYGen"), yGen);
Expand All @@ -523,6 +551,7 @@
registry.fill(HIST("MC/generated/signal/hPhiVsPtGenSig"), particle.phi(), ptGen);

if (particle.originMcGen() == RecoDecay::OriginType::Prompt) {
registry.get<THnSparse>(HIST("hnLcVarsGen"))->Fill(ptGen, yGen, numPvContributors, ptGenB, originType);
registry.fill(HIST("MC/generated/prompt/hPtGenPrompt"), ptGen);
registry.fill(HIST("MC/generated/prompt/hEtaGenPrompt"), particle.eta());
registry.fill(HIST("MC/generated/prompt/hYGenPrompt"), yGen);
Expand All @@ -532,6 +561,8 @@
registry.fill(HIST("MC/generated/prompt/hPhiVsPtGenSigPrompt"), particle.phi(), ptGen);
}
if (particle.originMcGen() == RecoDecay::OriginType::NonPrompt) {
ptGenB = mcParticles.rawIteratorAt(particle.idxBhadMotherPart()).pt();
registry.get<THnSparse>(HIST("hnLcVarsGen"))->Fill(ptGen, yGen, numPvContributors, ptGenB, originType);
registry.fill(HIST("MC/generated/nonprompt/hPtGenNonPrompt"), ptGen);
registry.fill(HIST("MC/generated/nonprompt/hEtaGenNonPrompt"), particle.eta());
registry.fill(HIST("MC/generated/nonprompt/hYGenNonPrompt"), yGen);
Expand Down Expand Up @@ -627,9 +658,9 @@
outputFD = candidate.mlProbLcToPKPi()[2]; /// non-prompt score
}
/// Fill the ML outputScores and variables of candidate
registry.get<THnSparse>(HIST("hnLcVarsWithBdt"))->Fill(massLc, pt, cent, outputBkg, outputPrompt, outputFD, 0);
registry.get<THnSparse>(HIST("hnLcVarsWithBdt"))->Fill(massLc, pt, cent, outputBkg, outputPrompt, outputFD, numPvContributors);
} else {
registry.get<THnSparse>(HIST("hnLcVars"))->Fill(massLc, pt, cent, ptProng0, ptProng1, ptProng2, chi2PCA, decayLength, cpa, 0);
registry.get<THnSparse>(HIST("hnLcVars"))->Fill(massLc, pt, cent, ptProng0, ptProng1, ptProng2, chi2PCA, decayLength, cpa, numPvContributors);
}
}
if (candidate.isSelLcToPiKP() >= selectionFlagLc) {
Expand All @@ -642,9 +673,9 @@
outputFD = candidate.mlProbLcToPiKP()[2]; /// non-prompt score
}
/// Fill the ML outputScores and variables of candidate
registry.get<THnSparse>(HIST("hnLcVarsWithBdt"))->Fill(massLc, pt, cent, outputBkg, outputPrompt, outputFD, 0);
registry.get<THnSparse>(HIST("hnLcVarsWithBdt"))->Fill(massLc, pt, cent, outputBkg, outputPrompt, outputFD, numPvContributors);
} else {
registry.get<THnSparse>(HIST("hnLcVars"))->Fill(massLc, pt, cent, ptProng0, ptProng1, ptProng2, chi2PCA, decayLength, cpa, 0);
registry.get<THnSparse>(HIST("hnLcVars"))->Fill(massLc, pt, cent, ptProng0, ptProng1, ptProng2, chi2PCA, decayLength, cpa, numPvContributors);
}
}
}
Expand Down Expand Up @@ -675,7 +706,7 @@
fillHistosMcRec<fillMl>(collision, candidates, mcParticles);
}
// MC gen.
fillHistosMcGen(mcParticles);
fillHistosMcGen(mcParticles, collisions);
}

void processDataStd(Collisions const& collisions,
Expand Down
Loading