Skip to content

Commit 4653bdb

Browse files
Mingyu3360715Mingyu3360715
authored
[PWGHF] small bug fixed in sigmaC task, add Thns for gen level (#10842)
1 parent 8dcabfb commit 4653bdb

File tree

1 file changed

+84
-24
lines changed

1 file changed

+84
-24
lines changed

‎PWGHF/D2H/Tasks/taskSigmac.cxx‎

Lines changed: 84 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -37,11 +37,14 @@ struct HfTaskSigmac {
3737
/// Properly normalize your results to provide a cross section
3838
/// OR
3939
/// consider the new parametrization of the fiducial acceptance (to be seen for reco signal in MC)
40-
Configurable<float> yCandMax{"yCandMax", -1, "Sc rapidity"};
40+
Configurable<float> yCandGenMax{"yCandGenMax", -1, "Maximum generated Sc rapidity"};
41+
Configurable<float> yCandRecoMax{"yCandRecoMax", -1, "Maximum Sc candidate rapidity"};
4142

4243
/// THn for candidate Λc+ and Σc0,++ cut variation
4344
Configurable<bool> enableTHn{"enableTHn", false, "enable the usage of THn for Λc+ and Σc0,++"};
4445
ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {16, 0, 16}, ""};
46+
ConfigurableAxis thnConfigAxisGenPt{"thnConfigAxisGenPt", {240, 0, 24}, "Gen pt prompt"};
47+
ConfigurableAxis thnConfigAxisGenPtB{"thnConfigAxisGenPtB", {800, 0, 80}, "Gen pt non-prompt"};
4548
ConfigurableAxis thnConfigAxisDecLength{"thnConfigAxisDecLength", {10, 0, 0.05}, ""};
4649
ConfigurableAxis thnConfigAxisDecLengthXY{"thnConfigAxisDecLengthXY", {10, 0, 0.05}, ""};
4750
ConfigurableAxis thnConfigAxisCPA{"thnConfigAxisCPA", {20, 0.8, 1}, ""};
@@ -255,12 +258,36 @@ struct HfTaskSigmac {
255258
const AxisSpec thnAxisChannel{4, -0.5, 3.5, "0: direct 1,2,3: resonant"};
256259
const AxisSpec thnAxisBdtScoreLcBkg{thnConfigAxisBdtScoreLcBkg, "BDT bkg score (Lc)"};
257260
const AxisSpec thnAxisBdtScoreLcNonPrompt{thnConfigAxisBdtScoreLcNonPrompt, "BDT non-prompt score (Lc)"};
258-
if (doprocessDataWithMl || doprocessMcWithMl) {
259-
registry.add("hnLambdaC", "THn for Lambdac", HistType::kTHnSparseF, {thnAxisPtLambdaC, thnAxisMassLambdaC, thnAxisBdtScoreLcBkg, thnAxisBdtScoreLcNonPrompt, thnAxisOriginMc, thnAxisChannel});
260-
registry.add("hnSigmaC", "THn for Sigmac", HistType::kTHnSparseF, {thnAxisPtLambdaC, axisDeltaMassSigmaC, thnAxisBdtScoreLcBkg, thnAxisBdtScoreLcNonPrompt, thnAxisOriginMc, thnAxisChannel, thnAxisPtSigmaC, thnAxisChargeSigmaC});
261+
const AxisSpec thnAxisGenPtLambdaC{thnConfigAxisGenPt, "#it{p}_{T}^{gen}(#Lambda_{c}^{+}) (GeV/#it{c})"};
262+
const AxisSpec thnAxisGenPtSigmaC{thnConfigAxisGenPt, "#it{p}_{T}^{gen}(#Sigma_{c}^{0,++}) (GeV/#it{c})"};
263+
const AxisSpec thnAxisGenPtLambdaCBMother{thnConfigAxisGenPtB, "#it{p}_{T}^{gen}(#Lambda_{c}^{+} B mother) (GeV/#it{c})"};
264+
const AxisSpec thnAxisGenPtSigmaCBMother{thnConfigAxisGenPtB, "#it{p}_{T}^{gen}(#Sigma_{c}^{0,++} B mother) (GeV/#it{c})"};
265+
std::vector<AxisSpec> axesLambdaCWithMl = {thnAxisPtLambdaC, thnAxisMassLambdaC, thnAxisBdtScoreLcBkg, thnAxisBdtScoreLcNonPrompt, thnAxisOriginMc, thnAxisChannel};
266+
std::vector<AxisSpec> axesSigmaCWithMl = {thnAxisPtLambdaC, axisDeltaMassSigmaC, thnAxisBdtScoreLcBkg, thnAxisBdtScoreLcNonPrompt, thnAxisOriginMc, thnAxisChannel, thnAxisPtSigmaC, thnAxisChargeSigmaC};
267+
std::vector<AxisSpec> axesLambdaCWoMl = {thnAxisPtLambdaC, thnAxisMassLambdaC, thnAxisDecLength, thnAxisDecLengthXY, thnAxisCPA, thnAxisCPAXY, thnAxisOriginMc, thnAxisChannel};
268+
std::vector<AxisSpec> axesSigmaCWoMl = {thnAxisPtLambdaC, axisDeltaMassSigmaC, thnAxisDecLength, thnAxisDecLengthXY, thnAxisCPA, thnAxisCPAXY, thnAxisOriginMc, thnAxisChannel, thnAxisPtSigmaC, thnAxisChargeSigmaC};
269+
if (isMc) {
270+
registry.add("hnLambdaCGen", "THn for Lambdac gen", HistType::kTHnSparseF, {thnAxisGenPtLambdaC, thnAxisGenPtLambdaCBMother, thnAxisOriginMc, thnAxisChannel});
271+
registry.add("hnSigmaCGen", "THn for Sigmac gen", HistType::kTHnSparseF, {thnAxisGenPtSigmaC, thnAxisGenPtSigmaCBMother, thnAxisOriginMc, thnAxisChannel, thnAxisGenPtLambdaC, thnAxisChargeSigmaC});
272+
if (doprocessMcWithMl) {
273+
axesLambdaCWithMl.push_back(thnAxisGenPtLambdaCBMother);
274+
axesSigmaCWithMl.push_back(thnAxisGenPtSigmaCBMother);
275+
registry.add("hnLambdaC", "THn for Lambdac", HistType::kTHnSparseF, axesLambdaCWithMl);
276+
registry.add("hnSigmaC", "THn for Sigmac", HistType::kTHnSparseF, axesSigmaCWithMl);
277+
} else {
278+
axesLambdaCWoMl.push_back(thnAxisGenPtLambdaCBMother);
279+
axesSigmaCWoMl.push_back(thnAxisGenPtSigmaCBMother);
280+
registry.add("hnLambdaC", "THn for Lambdac", HistType::kTHnSparseF, axesLambdaCWoMl);
281+
registry.add("hnSigmaC", "THn for Sigmac", HistType::kTHnSparseF, axesSigmaCWoMl);
282+
}
261283
} else {
262-
registry.add("hnLambdaC", "THn for Lambdac", HistType::kTHnSparseF, {thnAxisPtLambdaC, thnAxisMassLambdaC, thnAxisDecLength, thnAxisDecLengthXY, thnAxisCPA, thnAxisCPAXY, thnAxisOriginMc, thnAxisChannel});
263-
registry.add("hnSigmaC", "THn for Sigmac", HistType::kTHnSparseF, {thnAxisPtLambdaC, axisDeltaMassSigmaC, thnAxisDecLength, thnAxisDecLengthXY, thnAxisCPA, thnAxisCPAXY, thnAxisOriginMc, thnAxisChannel, thnAxisPtSigmaC, thnAxisChargeSigmaC});
284+
if (doprocessDataWithMl) {
285+
registry.add("hnLambdaC", "THn for Lambdac", HistType::kTHnSparseF, axesLambdaCWithMl);
286+
registry.add("hnSigmaC", "THn for Sigmac", HistType::kTHnSparseF, axesSigmaCWithMl);
287+
} else {
288+
registry.add("hnLambdaC", "THn for Lambdac", HistType::kTHnSparseF, axesLambdaCWoMl);
289+
registry.add("hnSigmaC", "THn for Sigmac", HistType::kTHnSparseF, axesSigmaCWoMl);
290+
}
264291
}
265292
}
266293

@@ -522,7 +549,7 @@ struct HfTaskSigmac {
522549
aod::TracksWMc const&)
523550
{
524551

525-
/// MC generated particles
552+
/// loop over Sc generated particles
526553
for (const auto& particle : mcParticlesSc) {
527554

528555
/// reject immediately particles different from Σc0,++
@@ -542,14 +569,15 @@ struct HfTaskSigmac {
542569
OR
543570
consider the new parametrization of the fiducial acceptance (to be seen for reco signal in MC)
544571
*/
545-
if (yCandMax >= 0. && std::abs(RecoDecay::y(particle.pVector(), o2::constants::physics::MassSigmaC0)) > yCandMax) {
572+
if (yCandGenMax >= 0. && std::abs(RecoDecay::y(particle.pVector(), o2::constants::physics::MassSigmaC0)) > yCandGenMax) {
546573
continue;
547574
}
548575

549576
/// Get the kinematic information of Σc0,++ and the daughters
550577
/// Get information about origin (prompt, non-prompt)
551578
/// Get information about decay Λc+ channel (direct, resonant)
552579
double ptGenSc(particle.pt()), etaGenSc(particle.eta()), phiGenSc(particle.phi());
580+
double ptGenScBMother(-1.);
553581
auto arrayDaughtersIds = particle.daughtersIds();
554582
if (arrayDaughtersIds.size() != 2) {
555583
/// This should never happen
@@ -561,7 +589,8 @@ struct HfTaskSigmac {
561589
double phiGenLc(-1.), phiGenSoftPi(-1.);
562590
int origin = -1;
563591
int8_t channel = -1;
564-
if (std::abs(arrayDaughtersIds[0]) == o2::constants::physics::Pdg::kLambdaCPlus) {
592+
auto daughter0 = mcParticles.rawIteratorAt(arrayDaughtersIds[0]);
593+
if (std::abs(daughter0.pdgCode()) == o2::constants::physics::Pdg::kLambdaCPlus) {
565594
/// daughter 0 is the Λc+, daughter 1 the soft π
566595
auto daugLc = mcParticlesLc.rawIteratorAt(arrayDaughtersIds[0]);
567596
auto daugSoftPi = mcParticles.rawIteratorAt(arrayDaughtersIds[1]);
@@ -573,7 +602,7 @@ struct HfTaskSigmac {
573602
ptGenSoftPi = daugSoftPi.pt();
574603
etaGenSoftPi = daugSoftPi.eta();
575604
phiGenSoftPi = daugSoftPi.phi();
576-
} else if (std::abs(arrayDaughtersIds[0]) == kPiPlus) {
605+
} else if (std::abs(daughter0.pdgCode()) == kPiPlus) {
577606
/// daughter 0 is the soft π, daughter 1 the Λc+
578607
auto daugLc = mcParticlesLc.rawIteratorAt(arrayDaughtersIds[1]);
579608
auto daugSoftPi = mcParticles.rawIteratorAt(arrayDaughtersIds[0]);
@@ -609,6 +638,12 @@ struct HfTaskSigmac {
609638
registry.fill(HIST("MC/generated/hPtGenLcFromSc0PlusPlusSig"), ptGenLc, origin, channel);
610639
registry.fill(HIST("MC/generated/hEtaGenLcFromSc0PlusPlusSig"), etaGenLc, origin, channel);
611640
registry.fill(HIST("MC/generated/hPhiGenLcFromSc0PlusPlusSig"), phiGenLc, origin, channel); /// Generated Λc+ ← Σc0,++ signal
641+
if (origin == RecoDecay::OriginType::Prompt) {
642+
registry.fill(HIST("hnSigmaCGen"), ptGenSc, ptGenScBMother, origin, channel, ptGenLc, 0);
643+
} else {
644+
ptGenScBMother = mcParticlesSc.rawIteratorAt(particle.idxBhadMotherPart()).pt();
645+
registry.fill(HIST("hnSigmaCGen"), ptGenSc, ptGenScBMother, origin, channel, ptGenLc, 0);
646+
}
612647
} else if (isScPlusPlusGen) {
613648
/// Generated Σc++ and Λc+ ← Σc++ signals
614649
registry.fill(HIST("MC/generated/hPtGenScPlusPlusSig"), ptGenSc, origin, channel);
@@ -630,9 +665,34 @@ struct HfTaskSigmac {
630665
registry.fill(HIST("MC/generated/hPtGenLcFromSc0PlusPlusSig"), ptGenLc, origin, channel);
631666
registry.fill(HIST("MC/generated/hEtaGenLcFromSc0PlusPlusSig"), etaGenLc, origin, channel);
632667
registry.fill(HIST("MC/generated/hPhiGenLcFromSc0PlusPlusSig"), phiGenLc, origin, channel); /// Generated Λc+ ← Σc0,++ signal
668+
if (origin == RecoDecay::OriginType::Prompt) {
669+
registry.fill(HIST("hnSigmaCGen"), ptGenSc, ptGenScBMother, origin, channel, ptGenLc, 2);
670+
} else {
671+
ptGenScBMother = mcParticlesSc.rawIteratorAt(particle.idxBhadMotherPart()).pt();
672+
registry.fill(HIST("hnSigmaCGen"), ptGenSc, ptGenScBMother, origin, channel, ptGenLc, 2);
673+
}
633674
}
634675

635-
} /// end loop over generated particles
676+
} /// end loop over Sc generated particles
677+
678+
/// loop over Lc generated particles
679+
for (const auto& particle : mcParticlesLc) {
680+
if (std::abs(particle.flagMcMatchGen()) != 1 << aod::hf_cand_3prong::DecayType::LcToPKPi) {
681+
continue;
682+
}
683+
if (yCandGenMax >= 0. && std::abs(RecoDecay::y(particle.pVector(), o2::constants::physics::MassLambdaCPlus)) > yCandGenMax) {
684+
continue;
685+
}
686+
double ptGenLc(-1.), ptGenLcBMother(-1.);
687+
int origin = particle.originMcGen();
688+
int channel = particle.flagMcDecayChanGen();
689+
if (origin == RecoDecay::OriginType::Prompt) {
690+
registry.fill(HIST("hnLambdaCGen"), ptGenLc, ptGenLcBMother, origin, channel);
691+
} else {
692+
ptGenLcBMother = mcParticlesLc.rawIteratorAt(particle.idxBhadMotherPart()).pt();
693+
registry.fill(HIST("hnLambdaCGen"), ptGenLc, ptGenLcBMother, origin, channel);
694+
}
695+
} /// end loop over Lc generated particles
636696

637697
/// reconstructed Σc0,++ matched to MC
638698
for (const auto& candSc : candidatesSc) {
@@ -642,7 +702,7 @@ struct HfTaskSigmac {
642702
continue;
643703
}
644704
/// rapidity selection on Σc0,++
645-
if (yCandMax >= 0. && std::abs(hfHelper.ySc0(candSc)) > yCandMax && std::abs(hfHelper.yScPlusPlus(candSc)) > yCandMax) {
705+
if (yCandRecoMax >= 0. && std::abs(hfHelper.ySc0(candSc)) > yCandRecoMax && std::abs(hfHelper.yScPlusPlus(candSc)) > yCandRecoMax) {
646706
continue;
647707
}
648708

@@ -748,10 +808,10 @@ struct HfTaskSigmac {
748808
outputMl.at(0) = candidateLc.mlProbLcToPKPi()[0]; /// bkg score
749809
outputMl.at(1) = candidateLc.mlProbLcToPKPi()[2]; /// non-prompt score
750810
}
751-
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, outputMl.at(0), outputMl.at(1), origin, channel, ptSc, std::abs(chargeSc));
811+
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, outputMl.at(0), outputMl.at(1), origin, channel, ptSc, std::abs(chargeSc), candSc.ptBhadMotherPart());
752812
} else {
753813
/// fill w/o BDT information
754-
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, decLengthLc, decLengthXYLc, cpaLc, cpaXYLc, origin, channel, ptSc, std::abs(chargeSc));
814+
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, decLengthLc, decLengthXYLc, cpaLc, cpaXYLc, origin, channel, ptSc, std::abs(chargeSc), candSc.ptBhadMotherPart());
755815
}
756816
}
757817

@@ -822,10 +882,10 @@ struct HfTaskSigmac {
822882
outputMl.at(0) = candidateLc.mlProbLcToPiKP()[0]; /// bkg score
823883
outputMl.at(1) = candidateLc.mlProbLcToPiKP()[2]; /// non-prompt score
824884
}
825-
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, outputMl.at(0), outputMl.at(1), origin, channel, ptSc, std::abs(chargeSc));
885+
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, outputMl.at(0), outputMl.at(1), origin, channel, ptSc, std::abs(chargeSc), candSc.ptBhadMotherPart());
826886
} else {
827887
/// fill w/o BDT information
828-
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, decLengthLc, decLengthXYLc, cpaLc, cpaXYLc, origin, channel, ptSc, std::abs(chargeSc));
888+
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, decLengthLc, decLengthXYLc, cpaLc, cpaXYLc, origin, channel, ptSc, std::abs(chargeSc), candSc.ptBhadMotherPart());
829889
}
830890
}
831891

@@ -923,10 +983,10 @@ struct HfTaskSigmac {
923983
outputMl.at(0) = candidateLc.mlProbLcToPKPi()[0]; /// bkg score
924984
outputMl.at(1) = candidateLc.mlProbLcToPKPi()[2]; /// non-prompt score
925985
}
926-
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, outputMl.at(0), outputMl.at(1), origin, channel, ptSc, std::abs(chargeSc));
986+
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, outputMl.at(0), outputMl.at(1), origin, channel, ptSc, std::abs(chargeSc), candSc.ptBhadMotherPart());
927987
} else {
928988
/// fill w/o BDT information
929-
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, decLengthLc, decLengthXYLc, cpaLc, cpaXYLc, origin, channel, ptSc, std::abs(chargeSc));
989+
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, decLengthLc, decLengthXYLc, cpaLc, cpaXYLc, origin, channel, ptSc, std::abs(chargeSc), candSc.ptBhadMotherPart());
930990
}
931991
}
932992

@@ -995,10 +1055,10 @@ struct HfTaskSigmac {
9951055
outputMl.at(0) = candidateLc.mlProbLcToPiKP()[0]; /// bkg score
9961056
outputMl.at(1) = candidateLc.mlProbLcToPiKP()[2]; /// non-prompt score
9971057
}
998-
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, outputMl.at(0), outputMl.at(1), origin, channel, ptSc, std::abs(chargeSc));
1058+
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, outputMl.at(0), outputMl.at(1), origin, channel, ptSc, std::abs(chargeSc), candSc.ptBhadMotherPart());
9991059
} else {
10001060
/// fill w/o BDT information
1001-
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, decLengthLc, decLengthXYLc, cpaLc, cpaXYLc, origin, channel, ptSc, std::abs(chargeSc));
1061+
registry.get<THnSparse>(HIST("hnSigmaC"))->Fill(ptLc, deltaMass, decLengthLc, decLengthXYLc, cpaLc, cpaXYLc, origin, channel, ptSc, std::abs(chargeSc), candSc.ptBhadMotherPart());
10021062
}
10031063
}
10041064

@@ -1034,10 +1094,10 @@ struct HfTaskSigmac {
10341094
outputMl.at(0) = candidateLc.mlProbLcToPKPi()[0]; /// bkg score
10351095
outputMl.at(1) = candidateLc.mlProbLcToPKPi()[2]; /// non-prompt score
10361096
}
1037-
registry.get<THnSparse>(HIST("hnLambdaC"))->Fill(ptLc, massLc, outputMl.at(0), outputMl.at(1), origin, channel);
1097+
registry.get<THnSparse>(HIST("hnLambdaC"))->Fill(ptLc, massLc, outputMl.at(0), outputMl.at(1), origin, channel, candidateLc.ptBhadMotherPart());
10381098
} else {
10391099
/// fill w/o BDT information
1040-
registry.get<THnSparse>(HIST("hnLambdaC"))->Fill(ptLc, massLc, decLengthLc, decLengthXYLc, cpaLc, cpaXYLc, origin, channel);
1100+
registry.get<THnSparse>(HIST("hnLambdaC"))->Fill(ptLc, massLc, decLengthLc, decLengthXYLc, cpaLc, cpaXYLc, origin, channel, candidateLc.ptBhadMotherPart());
10411101
}
10421102
}
10431103
if (candidateLc.isSelLcToPiKP() >= 1 && pdgAbs == kPiPlus) {
@@ -1050,10 +1110,10 @@ struct HfTaskSigmac {
10501110
outputMl.at(0) = candidateLc.mlProbLcToPiKP()[0]; /// bkg score
10511111
outputMl.at(1) = candidateLc.mlProbLcToPiKP()[2]; /// non-prompt score
10521112
}
1053-
registry.get<THnSparse>(HIST("hnLambdaC"))->Fill(ptLc, massLc, outputMl.at(0), outputMl.at(1), origin, channel);
1113+
registry.get<THnSparse>(HIST("hnLambdaC"))->Fill(ptLc, massLc, outputMl.at(0), outputMl.at(1), origin, channel, candidateLc.ptBhadMotherPart());
10541114
} else {
10551115
/// fill w/o BDT information
1056-
registry.get<THnSparse>(HIST("hnLambdaC"))->Fill(ptLc, massLc, decLengthLc, decLengthXYLc, cpaLc, cpaXYLc, origin, channel);
1116+
registry.get<THnSparse>(HIST("hnLambdaC"))->Fill(ptLc, massLc, decLengthLc, decLengthXYLc, cpaLc, cpaXYLc, origin, channel, candidateLc.ptBhadMotherPart());
10571117
}
10581118
}
10591119
}

0 commit comments

Comments
 (0)