Skip to content

Commit 9c3f8c0

Browse files
fgrosafgrosa
committed
Add generated ct and resonant channel info for c-deuteron analysis
1 parent bf425d5 commit 9c3f8c0

1 file changed

Lines changed: 72 additions & 46 deletions

File tree

‎PWGHF/D2H/Tasks/taskCd.cxx‎

Lines changed: 72 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -115,6 +115,8 @@ DECLARE_SOA_COLUMN(FlagMc, flagMc, int8_t); //! MC match
115115
DECLARE_SOA_COLUMN(OriginMcRec, originMcRec, int8_t); //! MC origin for reconstructed candidates
116116
DECLARE_SOA_COLUMN(FlagMcDecayChanRec, flagMcDecayChanRec, int8_t); //! Resonant MC decay channel for reconstructed candidates
117117
DECLARE_SOA_COLUMN(OriginMcGen, originMcGen, int8_t); //! MC origin for generated particles
118+
DECLARE_SOA_COLUMN(FlagMcDecayChanGen, flagMcDecayChanGen, int8_t); //! Resonant MC decay channel for generated candidates
119+
DECLARE_SOA_COLUMN(CtGen, ctGen, float); //! Generated ct computed wrt to c-deuteron production vertex, which can be either PV (prompt) or B-hadron decay vertex (non-prompt)
118120
DECLARE_SOA_COLUMN(Cent, cent, float); //! Centrality
119121
DECLARE_SOA_COLUMN(VtxZ, vtxZ, float); //! Vertex Z
120122
DECLARE_SOA_COLUMN(GIndexCol, gIndexCol, int); //! Global index for the collision
@@ -147,6 +149,7 @@ DECLARE_SOA_TABLE(HfCandCdLite, "AOD", "HFCANDCDLITE",
147149
full::FlagMc,
148150
full::OriginMcRec,
149151
full::FlagMcDecayChanRec,
152+
full::CtGen,
150153
full::Cent);
151154

152155
// full table for local Rotation & Event Mixing
@@ -179,6 +182,7 @@ DECLARE_SOA_TABLE(HfCandCdFull, "AOD", "HFCANDCDFULL",
179182
full::FlagMc,
180183
full::OriginMcRec,
181184
full::FlagMcDecayChanRec,
185+
full::CtGen,
182186
full::Cent,
183187
full::VtxZ,
184188
full::GIndexCol,
@@ -191,6 +195,8 @@ DECLARE_SOA_TABLE(HfCandCdGen, "AOD", "HFCANDCDGEN",
191195
full::Y,
192196
full::FlagMc,
193197
full::OriginMcGen,
198+
full::FlagMcDecayChanGen,
199+
full::CtGen,
194200
full::Cent,
195201
full::VtxZ,
196202
full::McCollisionId);
@@ -243,10 +249,11 @@ struct HfTaskCd {
243249
ConfigurableAxis thnAxisRapidity{"thnAxisRapidity", {20, -1, 1}, "Cand. rapidity bins"};
244250
ConfigurableAxis thnConfigAxisGenPtB{"thnConfigAxisGenPtB", {1000, 0, 100}, "Gen Pt B"};
245251
ConfigurableAxis thnConfigAxisNumPvContr{"thnConfigAxisNumPvContr", {200, -0.5, 199.5}, "Number of PV contributors"};
252+
ConfigurableAxis thnConfigAxisCt{"thnConfigAxisCt", {500, 0., 5000.}, ""};
246253

247-
constexpr static double MassCDeuteron = 3.226;
248254
constexpr static std::string_view SignalFolders[] = {"signal", "prompt", "nonprompt"};
249255
constexpr static std::string_view SignalSuffixes[] = {"", "Prompt", "NonPrompt"};
256+
const float cmToMum = 1.e4;
250257

251258
enum SignalClasses : int {
252259
Signal = 0,
@@ -257,30 +264,30 @@ struct HfTaskCd {
257264
HistogramRegistry registry{
258265
"registry",
259266
{/// mass candidate
260-
{"Data/hMass", "3-prong candidates;inv. mass (de K #pi) (GeV/#it{c}^{2})", {HistType::kTH1F, {{400, 2.4, 4.4}}}},
267+
{"Data/hMass", "3-prong candidates;inv. mass (de K #pi) (GeV/#it{c}^{2})", {HistType::kTH1D, {{400, 2.4, 4.4}}}},
261268
/// pT
262-
{"Data/hPt", "3-prong candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{360, 0., 36.}}}},
263-
{"Data/hPtProng0", "3-prong candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{360, 0., 36.}}}},
264-
{"Data/hPtProng1", "3-prong candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{360, 0., 36.}}}},
265-
{"Data/hPtProng2", "3-prong candidates;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{360, 0., 36.}}}},
269+
{"Data/hPt", "3-prong candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1D, {{360, 0., 36.}}}},
270+
{"Data/hPtProng0", "3-prong candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1D, {{360, 0., 36.}}}},
271+
{"Data/hPtProng1", "3-prong candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1D, {{360, 0., 36.}}}},
272+
{"Data/hPtProng2", "3-prong candidates;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1D, {{360, 0., 36.}}}},
266273
/// DCAxy to prim. vertex prongs
267-
{"Data/hd0Prong0", "3-prong candidates;prong 0 DCAxy to prim. vertex (cm);entries", {HistType::kTH1F, {{600, -0.4, 0.4}}}},
268-
{"Data/hd0Prong1", "3-prong candidates;prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH1F, {{600, -0.4, 0.4}}}},
269-
{"Data/hd0Prong2", "3-prong candidates;prong 2 DCAxy to prim. vertex (cm);entries", {HistType::kTH1F, {{600, -0.4, 0.4}}}},
274+
{"Data/hd0Prong0", "3-prong candidates;prong 0 DCAxy to prim. vertex (cm);entries", {HistType::kTH1D, {{600, -0.4, 0.4}}}},
275+
{"Data/hd0Prong1", "3-prong candidates;prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH1D, {{600, -0.4, 0.4}}}},
276+
{"Data/hd0Prong2", "3-prong candidates;prong 2 DCAxy to prim. vertex (cm);entries", {HistType::kTH1D, {{600, -0.4, 0.4}}}},
270277
/// decay length candidate
271-
{"Data/hDecLength", "3-prong candidates;decay length (cm);entries", {HistType::kTH1F, {{400, 0., 1.}}}},
278+
{"Data/hDecLength", "3-prong candidates;decay length (cm);entries", {HistType::kTH1D, {{400, 0., 1.}}}},
272279
/// decay length xy candidate
273-
{"Data/hDecLengthxy", "3-prong candidates;decay length xy (cm);entries", {HistType::kTH1F, {{400, 0., 1.}}}},
280+
{"Data/hDecLengthxy", "3-prong candidates;decay length xy (cm);entries", {HistType::kTH1D, {{400, 0., 1.}}}},
274281
/// cosine of pointing angle
275-
{"Data/hCPA", "3-prong candidates;cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}},
282+
{"Data/hCPA", "3-prong candidates;cosine of pointing angle;entries", {HistType::kTH1D, {{110, -1.1, 1.1}}}},
276283
/// cosine of pointing angle xy
277-
{"Data/hCPAxy", "3-prong candidates;cosine of pointing angle xy;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}},
284+
{"Data/hCPAxy", "3-prong candidates;cosine of pointing angle xy;entries", {HistType::kTH1D, {{110, -1.1, 1.1}}}},
278285
/// Chi 2 PCA to sec. vertex
279-
{"Data/hDca2", "3-prong candidates;prong Chi2PCA to sec. vertex (cm);entries", {HistType::kTH1F, {{400, 0., 20.}}}},
286+
{"Data/hDca2", "3-prong candidates;prong Chi2PCA to sec. vertex (cm);entries", {HistType::kTH1D, {{400, 0., 20.}}}},
280287
/// eta
281-
{"Data/hEta", "3-prong candidates;#it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}},
288+
{"Data/hEta", "3-prong candidates;#it{#eta};entries", {HistType::kTH1D, {{100, -2., 2.}}}},
282289
/// phi
283-
{"Data/hPhi", "3-prong candidates;#it{#Phi};entries", {HistType::kTH1F, {{100, 0., 6.3}}}}}};
290+
{"Data/hPhi", "3-prong candidates;#it{#Phi};entries", {HistType::kTH1D, {{100, 0., 6.3}}}}}};
284291

285292
HistogramRegistry qaRegistry{"QAHistos", {}, OutputObjHandlingPolicy::AnalysisObject};
286293

@@ -308,28 +315,28 @@ struct HfTaskCd {
308315
}
309316
};
310317

311-
addHistogramsRec("hMass", "inv. mass (de K #pi) (GeV/#it{c}^{2})", "", {HistType::kTH1F, {{400, 2.4, 4.4}}});
312-
addHistogramsRec("hPt", "#it{p}_{T}^{rec.} (GeV/#it{c})", "entries", {HistType::kTH1F, {{360, 0., 36.}}});
313-
addHistogramsGen("hPt", "#it{p}_{T}^{gen.} (GeV/#it{c})", "entries", {HistType::kTH1F, {{360, 0., 36.}}});
318+
addHistogramsRec("hMass", "inv. mass (de K #pi) (GeV/#it{c}^{2})", "", {HistType::kTH1D, {{400, 2.4, 4.4}}});
319+
addHistogramsRec("hPt", "#it{p}_{T}^{rec.} (GeV/#it{c})", "entries", {HistType::kTH1D, {{360, 0., 36.}}});
320+
addHistogramsGen("hPt", "#it{p}_{T}^{gen.} (GeV/#it{c})", "entries", {HistType::kTH1D, {{360, 0., 36.}}});
314321
if (!isData) {
315-
registry.add("MC/generated/signal/hPtGenSig", "3-prong candidates (matched);#it{p}_{T}^{gen.} (GeV/#it{c});entries", {HistType::kTH1F, {{360, 0., 36.}}});
322+
registry.add("MC/generated/signal/hPtGenSig", "3-prong candidates (matched);#it{p}_{T}^{gen.} (GeV/#it{c});entries", {HistType::kTH1D, {{360, 0., 36.}}});
316323
}
317-
addHistogramsRec("hPtProng0", "prong 0 #it{p}_{T} (GeV/#it{c})", "entries", {HistType::kTH1F, {{360, 0., 36.}}});
318-
addHistogramsRec("hPtProng1", "prong 1 #it{p}_{T} (GeV/#it{c})", "entries", {HistType::kTH1F, {{360, 0., 36.}}});
319-
addHistogramsRec("hPtProng2", "prong 2 #it{p}_{T} (GeV/#it{c})", "entries", {HistType::kTH1F, {{360, 0., 36.}}});
320-
addHistogramsRec("hd0Prong0", "prong 0 DCAxy to prim. vertex (cm)", "entries", {HistType::kTH1F, {{600, -0.4, 0.4}}});
321-
addHistogramsRec("hd0Prong1", "prong 1 DCAxy to prim. vertex (cm)", "entries", {HistType::kTH1F, {{600, -0.4, 0.4}}});
322-
addHistogramsRec("hd0Prong2", "prong 2 DCAxy to prim. vertex (cm)", "entries", {HistType::kTH1F, {{600, -0.4, 0.4}}});
323-
addHistogramsRec("hDecLength", "decay length (cm)", "entries", {HistType::kTH1F, {{400, 0., 1.}}});
324-
addHistogramsRec("hDecLengthxy", "decay length xy (cm)", "entries", {HistType::kTH1F, {{400, 0., 1.}}});
325-
addHistogramsRec("hCPA", "cosine of pointing angle", "entries", {HistType::kTH1F, {{110, -1.1, 1.1}}});
326-
addHistogramsRec("hCPAxy", "cosine of pointing angle xy", "entries", {HistType::kTH1F, {{110, -1.1, 1.1}}});
327-
addHistogramsRec("hDca2", "prong Chi2PCA to sec. vertex (cm)", "entries", {HistType::kTH1F, {{400, 0., 20.}}});
328-
addHistogramsRec("hEta", "#it{#eta}", "entries", {HistType::kTH1F, {{100, -2., 2.}}});
329-
addHistogramsGen("hEta", "#it{#eta}", "entries", {HistType::kTH1F, {{100, -2., 2.}}});
330-
addHistogramsGen("hY", "#it{y}", "entries", {HistType::kTH1F, {{100, -2., 2.}}});
331-
addHistogramsRec("hPhi", "#it{#Phi}", "entries", {HistType::kTH1F, {{100, 0., 6.3}}});
332-
addHistogramsGen("hPhi", "#it{#Phi}", "entries", {HistType::kTH1F, {{100, 0., 6.3}}});
324+
addHistogramsRec("hPtProng0", "prong 0 #it{p}_{T} (GeV/#it{c})", "entries", {HistType::kTH1D, {{360, 0., 36.}}});
325+
addHistogramsRec("hPtProng1", "prong 1 #it{p}_{T} (GeV/#it{c})", "entries", {HistType::kTH1D, {{360, 0., 36.}}});
326+
addHistogramsRec("hPtProng2", "prong 2 #it{p}_{T} (GeV/#it{c})", "entries", {HistType::kTH1D, {{360, 0., 36.}}});
327+
addHistogramsRec("hd0Prong0", "prong 0 DCAxy to prim. vertex (cm)", "entries", {HistType::kTH1D, {{600, -0.4, 0.4}}});
328+
addHistogramsRec("hd0Prong1", "prong 1 DCAxy to prim. vertex (cm)", "entries", {HistType::kTH1D, {{600, -0.4, 0.4}}});
329+
addHistogramsRec("hd0Prong2", "prong 2 DCAxy to prim. vertex (cm)", "entries", {HistType::kTH1D, {{600, -0.4, 0.4}}});
330+
addHistogramsRec("hDecLength", "decay length (cm)", "entries", {HistType::kTH1D, {{400, 0., 1.}}});
331+
addHistogramsRec("hDecLengthxy", "decay length xy (cm)", "entries", {HistType::kTH1D, {{400, 0., 1.}}});
332+
addHistogramsRec("hCPA", "cosine of pointing angle", "entries", {HistType::kTH1D, {{110, -1.1, 1.1}}});
333+
addHistogramsRec("hCPAxy", "cosine of pointing angle xy", "entries", {HistType::kTH1D, {{110, -1.1, 1.1}}});
334+
addHistogramsRec("hDca2", "prong Chi2PCA to sec. vertex (cm)", "entries", {HistType::kTH1D, {{400, 0., 20.}}});
335+
addHistogramsRec("hEta", "#it{#eta}", "entries", {HistType::kTH1D, {{100, -2., 2.}}});
336+
addHistogramsGen("hEta", "#it{#eta}", "entries", {HistType::kTH1D, {{100, -2., 2.}}});
337+
addHistogramsGen("hY", "#it{y}", "entries", {HistType::kTH1D, {{100, -2., 2.}}});
338+
addHistogramsRec("hPhi", "#it{#Phi}", "entries", {HistType::kTH1D, {{100, 0., 6.3}}});
339+
addHistogramsGen("hPhi", "#it{#Phi}", "entries", {HistType::kTH1D, {{100, 0., 6.3}}});
333340

334341
/// mass candidate
335342
if (isData) {
@@ -398,10 +405,11 @@ struct HfTaskCd {
398405
const AxisSpec thnAxisCentrality{thnConfigAxisCentrality, "centrality (FT0C)"};
399406
const AxisSpec thnAxisY{thnAxisRapidity, "rapidity"};
400407
const AxisSpec thnAxisPtB{thnConfigAxisGenPtB, "#it{p}_{T}^{B} (GeV/#it{c})"};
408+
const AxisSpec thnAxisCt{thnConfigAxisCt, "#it{ct} (#mum)"};
401409
const AxisSpec thnAxisTracklets{thnConfigAxisNumPvContr, "Number of PV contributors"};
402410

403411
std::vector axesStd{thnAxisMass, thnAxisPt, thnAxisPtProng0, thnAxisPtProng1, thnAxisPtProng2, thnAxisChi2PCA, thnAxisDecLength, thnAxisCPA, thnAxisCentrality};
404-
std::vector axesGen{thnAxisPt, thnAxisCentrality, thnAxisY, thnAxisTracklets, thnAxisPtB};
412+
std::vector axesGen{thnAxisPt, thnAxisCentrality, thnAxisY, thnAxisTracklets, thnAxisCt, thnAxisPtB};
405413
registry.add("hnCdVars", isData ? "THn for Reconstructed Cd candidates for data" : "THn for Reconstructed Cd candidates for MC", HistType::kTHnSparseF, axesStd);
406414
if (!isData) {
407415
registry.add("hnCdVarsGen", "THn for Generated Cd", HistType::kTHnSparseF, axesGen);
@@ -476,14 +484,28 @@ struct HfTaskCd {
476484
const int64_t timeStamp = bc.timestamp();
477485

478486
for (const auto& candidate : groupedCdCandidates) {
487+
if (candidate.flagMcMatchRec() == 0) { // we skip combinatorial background
488+
continue;
489+
}
479490
if (!TESTBIT(candidate.hfflag(), aod::hf_cand_3prong::DecayType::CdToDeKPi)) {
480491
continue;
481492
}
482-
const auto yCd = RecoDecay::y(candidate.pVector(), MassCDeuteron);
493+
const auto yCd = RecoDecay::y(candidate.pVector(), o2::constants::physics::MassCDeuteron);
483494
if (yCandRecoMax >= 0. && std::abs(yCd) > yCandRecoMax) {
484495
continue;
485496
}
486497

498+
float ctGen{-1.f}, ptGen{-1.f};
499+
int pdgCodeProng0{0};
500+
if (candidate.flagMcMatchRec() == hf_decay::hf_cand_3prong::DecayChannelMain::CDeuteronToDeKPi) {
501+
const auto& mcParticleProng0 = candidate.template prong0_as<HFTracksMc>().template mcParticle_as<CandCdMcGen>();
502+
pdgCodeProng0 = std::abs(mcParticleProng0.pdgCode());
503+
const auto indexMother = RecoDecay::getMother(mcParticles, mcParticleProng0, o2::constants::physics::Pdg::kCDeuteron, true);
504+
const auto particleMother = mcParticles.rawIteratorAt(indexMother);
505+
ctGen = RecoDecay::ct(std::array{particleMother.px(), particleMother.py(), particleMother.pz()}, RecoDecay::distance(std::array{particleMother.vx(), particleMother.vy(), particleMother.vz()}, std::array{mcParticleProng0.vx(), mcParticleProng0.vy(), mcParticleProng0.vz()}), o2::constants::physics::MassCDeuteron) * cmToMum;
506+
ptGen = particleMother.pt();
507+
}
508+
487509
if (fillCandLiteTree || fillCandFullTree) {
488510
float invMassCd = 0.f;
489511
float invMassLc = 0.f;
@@ -571,6 +593,7 @@ struct HfTaskCd {
571593
candidate.flagMcMatchRec(),
572594
candidate.originMcRec(),
573595
candidate.flagMcDecayChanRec(),
596+
ctGen,
574597
o2::hf_centrality::getCentralityColl(collision));
575598
}
576599

@@ -604,6 +627,7 @@ struct HfTaskCd {
604627
candidate.flagMcMatchRec(),
605628
candidate.originMcRec(),
606629
candidate.flagMcDecayChanRec(),
630+
ctGen,
607631
o2::hf_centrality::getCentralityColl(collision),
608632
collision.posZ(),
609633
collision.globalIndex(),
@@ -615,11 +639,7 @@ struct HfTaskCd {
615639
continue;
616640
}
617641

618-
const auto& mcParticleProng0 = candidate.template prong0_as<HFTracksMc>().template mcParticle_as<McParticles3ProngMatched>();
619-
const auto pdgCodeProng0 = std::abs(mcParticleProng0.pdgCode());
620-
const auto indexMother = RecoDecay::getMother(mcParticles, mcParticleProng0, o2::constants::physics::Pdg::kCDeuteron, true);
621-
const auto particleMother = mcParticles.rawIteratorAt(indexMother);
622-
registry.fill(HIST("MC/generated/signal/hPtGenSig"), particleMother.pt());
642+
registry.fill(HIST("MC/generated/signal/hPtGenSig"), ptGen);
623643

624644
fillHistogramsRecSig<Signal>(candidate);
625645
if (candidate.originMcRec() == RecoDecay::OriginType::Prompt) {
@@ -661,10 +681,10 @@ struct HfTaskCd {
661681
void fillHistosMcGen(CandCdMcGen const& mcParticles, Coll const& recoCollisions)
662682
{
663683
for (const auto& particle : mcParticles) {
664-
if (std::abs(particle.flagMcMatchGen()) != hf_decay::hf_cand_3prong::DecayChannelMain::CDeuteronToDeKPi) {
684+
if (std::abs(particle.flagMcMatchGen()) == hf_decay::hf_cand_3prong::DecayChannelMain::CDeuteronToDeKPi) {
665685
continue;
666686
}
667-
const auto yGen = RecoDecay::y(particle.pVector(), MassCDeuteron);
687+
const auto yGen = RecoDecay::y(particle.pVector(), o2::constants::physics::MassCDeuteron);
668688
if (yCandGenMax >= 0. && std::abs(yGen) > yCandGenMax) {
669689
continue;
670690
}
@@ -677,6 +697,8 @@ struct HfTaskCd {
677697
}
678698
const float cent = o2::hf_centrality::getCentralityGenColl(recoCollsPerMcColl);
679699
const float ptGenB = particle.originMcGen() == RecoDecay::OriginType::Prompt ? -1.f : mcParticles.rawIteratorAt(particle.idxBhadMotherPart()).pt();
700+
const auto firstDau = particle.template daughters_as<CandCdMcGen>().begin();
701+
const float ctGen = RecoDecay::ct(std::array{particle.px(), particle.py(), particle.pz()}, RecoDecay::distance(std::array{particle.vx(), particle.vy(), particle.vz()}, std::array{firstDau.vx(), firstDau.vy(), firstDau.vz()}), o2::constants::physics::MassCDeuteron) * cmToMum;
680702

681703
fillHistogramsGen<Signal>(particle, yGen);
682704
if (particle.originMcGen() == RecoDecay::OriginType::Prompt) {
@@ -686,7 +708,7 @@ struct HfTaskCd {
686708
}
687709

688710
if (fillTHn) {
689-
std::vector<double> valuesToFill{particle.pt(), cent, yGen, static_cast<double>(numPvContributors), ptGenB};
711+
std::vector<double> valuesToFill{particle.pt(), cent, yGen, static_cast<double>(numPvContributors), ctGen, ptGenB};
690712
registry.get<THnSparse>(HIST("hnCdVarsGen"))->Fill(valuesToFill.data());
691713
}
692714

@@ -697,6 +719,8 @@ struct HfTaskCd {
697719
yGen,
698720
particle.flagMcMatchGen(),
699721
particle.originMcGen(),
722+
particle.flagMcDecayChanGen(),
723+
ctGen,
700724
cent,
701725
vtxZ,
702726
particle.mcCollision().globalIndex());
@@ -921,6 +945,7 @@ struct HfTaskCd {
921945
0,
922946
0,
923947
-1,
948+
-1.f,
924949
cent);
925950
}
926951

@@ -955,6 +980,7 @@ struct HfTaskCd {
955980
0,
956981
0,
957982
-1,
983+
-1.f,
958984
cent,
959985
collision.posZ(),
960986
collision.globalIndex(),

0 commit comments

Comments
 (0)