Skip to content

Commit 30b9e1c

Browse files
EloviyoShirajum Monira
Eloviyo
and
Shirajum Monira
authored
process functions for MC efficiency on derived data (#7048)
Co-authored-by: Shirajum Monira <shirajum.monira@cernch>
1 parent e1d9a7f commit 30b9e1c

1 file changed

Lines changed: 174 additions & 0 deletions

File tree

‎PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackV0Extended.cxx‎

Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,16 +28,20 @@
2828
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseContainer.h"
2929
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseDetaDphiStar.h"
3030
#include "PWGCF/FemtoUniverse/Core/FemtoUtils.h"
31+
#include "Framework/O2DatabasePDGPlugin.h"
3132

3233
using namespace o2;
3334
using namespace o2::soa;
3435
using namespace o2::framework;
3536
using namespace o2::framework::expressions;
3637
using namespace o2::analysis::femtoUniverse;
3738
using namespace o2::aod::pidutils;
39+
using namespace o2::track;
3840

3941
struct femtoUniversePairTaskTrackV0Extended {
4042

43+
Service<o2::framework::O2DatabasePDG> pdgMC;
44+
4145
SliceCache cache;
4246
using FemtoFullParticles = soa::Join<aod::FDParticles, aod::FDExtParticles>;
4347
Preslice<FemtoFullParticles> perCol = aod::femtouniverseparticle::fdCollisionId;
@@ -122,6 +126,11 @@ struct femtoUniversePairTaskTrackV0Extended {
122126
ConfigurableAxis ConfmTBins3D{"ConfmTBins3D", {VARIABLE_WIDTH, 1.02f, 1.14f, 1.20f, 1.26f, 1.38f, 1.56f, 1.86f, 4.50f}, "mT Binning for the 3Dimensional plot: k* vs multiplicity vs mT (set <<ConfUse3D>> to true in order to use)"};
123127
ConfigurableAxis ConfmultBins3D{"ConfMultBins3D", {VARIABLE_WIDTH, 0.0f, 20.0f, 30.0f, 40.0f, 99999.0f}, "multiplicity Binning for the 3Dimensional plot: k* vs multiplicity vs mT (set <<ConfUse3D>> to true in order to use)"};
124128

129+
/// MC
130+
131+
Configurable<float> ConfHPtMC{"ConfHPtMC", 4.0f, "higher limit for pt"};
132+
Configurable<float> ConfLPtMC{"ConfLPtMC", 0.3f, "lower limit for pt"};
133+
125134
static constexpr UInt_t V0ChildTable[][2] = {{0, 1}, {1, 0}, {1, 1}}; // Table to select the V0 children
126135

127136
FemtoUniverseContainer<femtoUniverseContainer::EventType::same, femtoUniverseContainer::Observable::kstar> sameEventCont;
@@ -133,6 +142,8 @@ struct femtoUniversePairTaskTrackV0Extended {
133142
/// Histogram output
134143
HistogramRegistry qaRegistry{"TrackQA", {}, OutputObjHandlingPolicy::AnalysisObject};
135144
HistogramRegistry resultRegistry{"Correlations", {}, OutputObjHandlingPolicy::AnalysisObject};
145+
HistogramRegistry registryMCtruth{"MCtruthHistos", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
146+
HistogramRegistry registryMCreco{"MCrecoHistos", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
136147

137148
bool IsNSigmaCombined(float mom, float nsigmaTPCParticle, float nsigmaTOFParticle)
138149
{
@@ -206,6 +217,46 @@ struct femtoUniversePairTaskTrackV0Extended {
206217
posChildV0Type2.init(&qaRegistry, ConfChildTempFitVarpTBins, ConfChildTempFitVarBins, false, 0, true, "posChildV0Type2");
207218
negChildV0Type2.init(&qaRegistry, ConfChildTempFitVarpTBins, ConfChildTempFitVarBins, false, 0, true, "negChildV0Type2");
208219

220+
// MC truth
221+
registryMCtruth.add("plus/MCtruthLambda", "MC truth Lambdas;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
222+
registryMCtruth.add("minus/MCtruthLambda", "MC truth Lambdas;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
223+
224+
registryMCtruth.add("plus/MCtruthAllPt", "MC truth all;#it{p}_{T} (GeV/c); #eta", {HistType::kTH1F, {{500, 0, 5}}});
225+
registryMCtruth.add("minus/MCtruthAllPt", "MC truth all;#it{p}_{T} (GeV/c); #eta", {HistType::kTH1F, {{500, 0, 5}}});
226+
227+
registryMCtruth.add("plus/MCtruthPi", "MC truth pions;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
228+
registryMCtruth.add("plus/MCtruthPr", "MC truth protons;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
229+
230+
registryMCtruth.add("minus/MCtruthPi", "MC truth pions;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
231+
registryMCtruth.add("minus/MCtruthPr", "MC truth protons;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
232+
233+
registryMCtruth.add("plus/MCtruthPiPt", "MC truth pions;#it{p}_{T} (GeV/c)", {HistType::kTH1F, {{500, 0, 5}}});
234+
registryMCtruth.add("plus/MCtruthPrPt", "MC truth protons;#it{p}_{T} (GeV/c)", {HistType::kTH1F, {{500, 0, 5}}});
235+
registryMCtruth.add("minus/MCtruthPiPt", "MC truth pions;#it{p}_{T} (GeV/c)", {HistType::kTH1F, {{500, 0, 5}}});
236+
registryMCtruth.add("minus/MCtruthPrPt", "MC truth protons;#it{p}_{T} (GeV/c)", {HistType::kTH1F, {{500, 0, 5}}});
237+
238+
// MC reco
239+
registryMCreco.add("plus/MCrecoLambda", "MC reco Lambdas;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
240+
registryMCreco.add("plus/MCrecoLambdaChildPr", "MC reco Lambdas;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
241+
registryMCreco.add("plus/MCrecoLambdaChildPi", "MC reco Lambdas;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
242+
registryMCreco.add("minus/MCrecoLambda", "MC reco Lambdas;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
243+
registryMCreco.add("minus/MCrecoLambdaChildPr", "MC reco Lambdas;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
244+
registryMCreco.add("minus/MCrecoLambdaChildPi", "MC reco Lambdas;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
245+
246+
registryMCreco.add("plus/MCrecoAllPt", "MC reco all;#it{p}_{T} (GeV/c); #eta", {HistType::kTH1F, {{500, 0, 5}}});
247+
registryMCreco.add("minus/MCrecoAllPt", "MC reco all;#it{p}_{T} (GeV/c); #eta", {HistType::kTH1F, {{500, 0, 5}}});
248+
249+
registryMCreco.add("plus/MCrecoPi", "MC reco pions;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
250+
registryMCreco.add("plus/MCrecoPr", "MC reco protons;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
251+
252+
registryMCreco.add("minus/MCrecoPi", "MC reco pions;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
253+
registryMCreco.add("minus/MCrecoPr", "MC reco protons;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
254+
255+
registryMCreco.add("plus/MCrecoPiPt", "MC reco pions;#it{p}_{T} (GeV/c)", {HistType::kTH1F, {{500, 0, 5}}});
256+
registryMCreco.add("plus/MCrecoPrPt", "MC reco protons;#it{p}_{T} (GeV/c)", {HistType::kTH1F, {{500, 0, 5}}});
257+
registryMCreco.add("minus/MCrecoPiPt", "MC reco pions;#it{p}_{T} (GeV/c)", {HistType::kTH1F, {{500, 0, 5}}});
258+
registryMCreco.add("minus/MCrecoPrPt", "MC reco protons;#it{p}_{T} (GeV/c)", {HistType::kTH1F, {{500, 0, 5}}});
259+
209260
sameEventCont.init(&resultRegistry, ConfkstarBins, ConfMultBins, ConfkTBins, ConfmTBins, ConfmultBins3D, ConfmTBins3D, ConfEtaBins, ConfPhiBins, ConfIsMC, ConfUse3D);
210261
sameEventCont.setPDGCodes(ConfTrkPDGCodePartOne, ConfV0PDGCodePartTwo);
211262
mixedEventCont.init(&resultRegistry, ConfkstarBins, ConfMultBins, ConfkTBins, ConfmTBins, ConfmultBins3D, ConfmTBins3D, ConfEtaBins, ConfPhiBins, ConfIsMC, ConfUse3D);
@@ -459,6 +510,129 @@ struct femtoUniversePairTaskTrackV0Extended {
459510
}
460511

461512
PROCESS_SWITCH(femtoUniversePairTaskTrackV0Extended, processMixedEventV0, "Enable processing mixed events for V0 - V0", false);
513+
514+
///--------------------------------------------MC-------------------------------------------------///
515+
516+
using FemtoRecoParticles = soa::Join<aod::FDParticles, aod::FDExtParticles, aod::FDMCLabels>;
517+
518+
/// This function fills MC truth particles from derived MC table
519+
void processMCTruth(aod::FDParticles const& parts)
520+
{
521+
for (auto& part : parts) {
522+
if (!(part.partType() == uint8_t(aod::femtouniverseparticle::ParticleType::kMCTruthTrack)))
523+
continue;
524+
if (TMath::Abs(part.eta()) > ConfEta || part.pt() < ConfLPtMC || part.pt() > ConfHPtMC)
525+
continue;
526+
527+
int pdgCode = static_cast<int>(part.pidcut());
528+
const auto& pdgParticle = pdgMC->GetParticle(pdgCode);
529+
if (!pdgParticle) {
530+
continue;
531+
}
532+
533+
if (pdgCode == 3122) {
534+
registryMCtruth.fill(HIST("plus/MCtruthLambda"), part.pt(), part.eta());
535+
continue;
536+
} else if (pdgCode == -3122) {
537+
registryMCtruth.fill(HIST("minus/MCtruthLambda"), part.pt(), part.eta());
538+
continue;
539+
}
540+
541+
if (pdgParticle->Charge() > 0) {
542+
registryMCtruth.fill(HIST("plus/MCtruthAllPt"), part.pt());
543+
}
544+
if (pdgCode == 211) {
545+
registryMCtruth.fill(HIST("plus/MCtruthPi"), part.pt(), part.eta());
546+
registryMCtruth.fill(HIST("plus/MCtruthPiPt"), part.pt());
547+
}
548+
if (pdgCode == 2212) {
549+
registryMCtruth.fill(HIST("plus/MCtruthPr"), part.pt(), part.eta());
550+
registryMCtruth.fill(HIST("plus/MCtruthPrPt"), part.pt());
551+
}
552+
553+
if (pdgParticle->Charge() < 0) {
554+
registryMCtruth.fill(HIST("minus/MCtruthAllPt"), part.pt());
555+
}
556+
if (pdgCode == -211) {
557+
registryMCtruth.fill(HIST("minus/MCtruthPi"), part.pt(), part.eta());
558+
registryMCtruth.fill(HIST("minus/MCtruthPiPt"), part.pt());
559+
}
560+
if (pdgCode == -2212) {
561+
registryMCtruth.fill(HIST("minus/MCtruthPr"), part.pt(), part.eta());
562+
registryMCtruth.fill(HIST("minus/MCtruthPrPt"), part.pt());
563+
}
564+
}
565+
}
566+
567+
PROCESS_SWITCH(femtoUniversePairTaskTrackV0Extended, processMCTruth, "Process MC truth data", false);
568+
569+
void processMCReco(FemtoRecoParticles const& parts, aod::FDMCParticles const& mcparts)
570+
{
571+
for (auto& part : parts) {
572+
auto mcPartId = part.fdMCParticleId();
573+
if (mcPartId == -1)
574+
continue; // no MC particle
575+
const auto& mcpart = mcparts.iteratorAt(mcPartId);
576+
//
577+
if (part.partType() == aod::femtouniverseparticle::ParticleType::kV0) {
578+
if (mcpart.pdgMCTruth() == 3122) {
579+
const auto& posChild = parts.iteratorAt(part.globalIndex() - 2);
580+
const auto& negChild = parts.iteratorAt(part.globalIndex() - 1);
581+
/// Daughters that do not pass this condition are not selected
582+
if (IsParticleTPC(posChild, 0) && IsParticleTPC(negChild, 1)) {
583+
registryMCreco.fill(HIST("plus/MCrecoLambda"), mcpart.pt(), mcpart.eta()); // lambda
584+
if (auto mcpartIdChild = posChild.fdMCParticleId(); mcpartIdChild != -1) {
585+
const auto& mcpartChild = mcparts.iteratorAt(mcpartIdChild);
586+
registryMCreco.fill(HIST("plus/MCrecoLambdaChildPr"), mcpartChild.pt(), mcpartChild.eta()); // lambda proton child
587+
}
588+
if (auto mcpartIdChild = negChild.fdMCParticleId(); mcpartIdChild != -1) {
589+
const auto& mcpartChild = mcparts.iteratorAt(mcpartIdChild);
590+
registryMCreco.fill(HIST("plus/MCrecoLambdaChildPi"), mcpartChild.pt(), mcpartChild.eta()); // lambda pion child
591+
}
592+
}
593+
} else if (mcpart.pdgMCTruth() == -3122) {
594+
const auto& posChild = parts.iteratorAt(part.globalIndex() - 2);
595+
const auto& negChild = parts.iteratorAt(part.globalIndex() - 1);
596+
/// Daughters that do not pass this condition are not selected
597+
if (IsParticleTPC(posChild, 1) && IsParticleTPC(negChild, 0)) {
598+
registryMCreco.fill(HIST("minus/MCrecoLambda"), mcpart.pt(), mcpart.eta()); // anti-lambda
599+
if (auto mcpartIdChild = posChild.fdMCParticleId(); mcpartIdChild != -1) {
600+
const auto& mcpartChild = mcparts.iteratorAt(mcpartIdChild);
601+
registryMCreco.fill(HIST("minus/MCrecoLambdaChildPi"), mcpartChild.pt(), mcpartChild.eta()); // anti-lambda pion child
602+
}
603+
if (auto mcpartIdChild = negChild.fdMCParticleId(); mcpartIdChild != -1) {
604+
const auto& mcpartChild = mcparts.iteratorAt(mcpartIdChild);
605+
registryMCreco.fill(HIST("minus/MCrecoLambdaChildPr"), mcpartChild.pt(), mcpartChild.eta()); // anti-lambda proton child
606+
}
607+
}
608+
}
609+
} else if (part.partType() == aod::femtouniverseparticle::ParticleType::kTrack) {
610+
if (part.sign() > 0) {
611+
registryMCreco.fill(HIST("plus/MCrecoAllPt"), mcpart.pt());
612+
if (mcpart.pdgMCTruth() == 211 && IsNSigmaCombined(part.p(), unPackInTable<aod::pidtpc_tiny::binning>(part.tpcNSigmaStorePi()), unPackInTable<aod::pidtof_tiny::binning>(part.tofNSigmaStorePi()))) {
613+
registryMCreco.fill(HIST("plus/MCrecoPi"), mcpart.pt(), mcpart.eta());
614+
registryMCreco.fill(HIST("plus/MCrecoPiPt"), mcpart.pt());
615+
} else if (mcpart.pdgMCTruth() == 2212 && IsNSigmaCombined(part.p(), unPackInTable<aod::pidtpc_tiny::binning>(part.tpcNSigmaStorePr()), unPackInTable<aod::pidtof_tiny::binning>(part.tofNSigmaStorePr()))) {
616+
registryMCreco.fill(HIST("plus/MCrecoPr"), mcpart.pt(), mcpart.eta());
617+
registryMCreco.fill(HIST("plus/MCrecoPrPt"), mcpart.pt());
618+
}
619+
}
620+
621+
if (part.sign() < 0) {
622+
registryMCreco.fill(HIST("minus/MCrecoAllPt"), part.pt());
623+
if (mcpart.pdgMCTruth() == -211 && IsNSigmaCombined(part.p(), unPackInTable<aod::pidtpc_tiny::binning>(part.tpcNSigmaStorePi()), unPackInTable<aod::pidtof_tiny::binning>(part.tofNSigmaStorePi()))) {
624+
registryMCreco.fill(HIST("plus/MCrecoPi"), mcpart.pt(), mcpart.eta());
625+
registryMCreco.fill(HIST("plus/MCrecoPiPt"), mcpart.pt());
626+
} else if (mcpart.pdgMCTruth() == -2212 && IsNSigmaCombined(part.p(), unPackInTable<aod::pidtpc_tiny::binning>(part.tpcNSigmaStorePr()), unPackInTable<aod::pidtof_tiny::binning>(part.tofNSigmaStorePr()))) {
627+
registryMCreco.fill(HIST("plus/MCrecoPr"), mcpart.pt(), mcpart.eta());
628+
registryMCreco.fill(HIST("plus/MCrecoPrPt"), mcpart.pt());
629+
}
630+
}
631+
} // partType
632+
}
633+
}
634+
635+
PROCESS_SWITCH(femtoUniversePairTaskTrackV0Extended, processMCReco, "Process MC reco data", false);
462636
};
463637

464638
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)