Skip to content

Commit b9847c6

Browse files
mhemmer-cernmhemmer-cern
committed
[PWGEM] emcalPhotonMcTask: Add three new TruthClasses for Bremsstrahlung gammas:
Also added purity for both clusters to the output table for possible later selection. - PhotonBSPhotonPair: photon + photon from Bremsstrahlung - ElectronBSPhotonPair: one cluster from Bremsstrahlung and one electron cluster except case BSPhotonElectron - BSPhotonPair: both photons from Bremsstrahlung
1 parent ede5910 commit b9847c6

2 files changed

Lines changed: 46 additions & 17 deletions

File tree

‎PWGEM/PhotonMeson/DataModel/ConversionMl.h‎

Lines changed: 12 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -39,15 +39,26 @@ DECLARE_SOA_COLUMN(Time1, time1, float); //! time of first cluster in
3939
DECLARE_SOA_COLUMN(Time2, time2, float); //! time of second cluster in ns
4040
DECLARE_SOA_COLUMN(NCells1, nCells1, uint8_t); //! NCells of first cluster
4141
DECLARE_SOA_COLUMN(NCells2, nCells2, uint8_t); //! NCells of second cluster
42-
DECLARE_SOA_COLUMN(TruthLabel, truthLabel, int8_t); //! truth label for ML training (0 == conversion, 1 == from Pi0, 2 == background)
42+
DECLARE_SOA_COLUMN(TruthLabel, truthLabel, int8_t); //! truth label for ML training -- see TruthClass enum for the mapping
4343
DECLARE_SOA_COLUMN(CentOrMult, centOrMult, float); //! centrality or multiplicity value of the collision
44+
DECLARE_SOA_COLUMN(Purity1, purity1, float); // leading-contributor amplitude fraction, cluster 1
45+
DECLARE_SOA_COLUMN(Purity2, purity2, float); // leading-contributor amplitude fraction, cluster 2
46+
4447
} // namespace convtag
4548
DECLARE_SOA_TABLE(ConvTagCandidates, "AOD", "CONVTAGCAND",
4649
convtag::PMEvent,
4750
convtag::Minv, convtag::HarmonicET, convtag::DeltaEta, convtag::DeltaPhi, convtag::Phiv,
4851
convtag::E1, convtag::E2, convtag::M021, convtag::M022,
4952
convtag::Time1, convtag::Time2, convtag::NCells1, convtag::NCells2,
5053
convtag::TruthLabel, convtag::CentOrMult);
54+
55+
DECLARE_SOA_TABLE_VERSIONED(ConvTagCandidates_001, "AOD", "CONVTAGCAND", 1,
56+
convtag::PMEvent,
57+
convtag::Minv, convtag::HarmonicET, convtag::DeltaEta, convtag::DeltaPhi, convtag::Phiv,
58+
convtag::E1, convtag::E2, convtag::M021, convtag::M022,
59+
convtag::Time1, convtag::Time2, convtag::NCells1, convtag::NCells2, convtag::Purity1, convtag::Purity2,
60+
convtag::TruthLabel, convtag::CentOrMult);
61+
5162
} // namespace o2::aod
5263

5364
#endif // PWGEM_PHOTONMESON_DATAMODEL_CONVERSIONML_H_

‎PWGEM/PhotonMeson/Tasks/emcalPhotonMcTask.cxx‎

Lines changed: 34 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -94,16 +94,19 @@ enum class TruthClass {
9494
PhotonElectronSamePi0, // 4: photon + electron cluster, same Pi0
9595
PhotonElectronDiffPi0, // 5: photon + electron cluster, different Pi0s
9696
PhotonElectronOnePi0, // 6: photon + electron cluster, only one from a Pi0
97-
PhotonElectronBS, // 7: photon + electron cluster, from Bremsstrahlung
97+
BSPhotonElectron, // 7: photon + electron cluster, from Bremsstrahlung
9898

9999
ElectronPairSamePi0, // 8: e+e cluster pair, same Pi0 (conversion and/or Dalitz)
100100
ElectronPairDiffPi0, // 9: e+e cluster pair, different Pi0s
101101
ElectronPairOnePi0, // 10: e+e cluster pair, only one from a Pi0
102102

103-
SplitPhotonCluster, // 11: one photon producing two clusters
104-
SplitLeptonCluster, // 12: one lepton producing two clusters
103+
SplitPhotonCluster, // 11: one photon producing two clusters
104+
SplitLeptonCluster, // 12: one lepton producing two clusters
105+
PhotonBSPhotonPair, // 13: photon + photon from Bremsstrahlung
106+
ElectronBSPhotonPair, // 14: one cluster from Bremsstrahlung and one electron cluster except case BSPhotonElectron
107+
BSPhotonPair, // 15: both photons from Bremsstrahlung
105108

106-
Background, // 13: else / uncorrelated
109+
Background, // 16: else / uncorrelated
107110

108111
NClasses
109112
};
@@ -132,15 +135,17 @@ struct ClusterMcInfo {
132135
bool isFromBremsstrahlung = false;
133136
int convMotherId = -1;
134137
int photonId = -1;
138+
float purity = 0;
135139
};
136140

137141
template <o2::soa::is_iterator TGroup, o2::soa::is_iterator TIter, o2::soa::is_table McParticles>
138142
ClusterMcInfo classifyCluster(const TGroup& g, TIter& mcCluster, TIter& mcClusterLooper, TIter& mcClusterLooper2, McParticles const& mcParticles)
139143
{
140144
ClusterMcInfo info;
141145
mcCluster.setCursor(g.emmcparticleIds()[0]);
142-
info.isFromBremsstrahlung = isFromBremsstrahlung(mcCluster, mcClusterLooper);
146+
info.isFromBremsstrahlung = isFromBremsstrahlung(mcCluster, mcClusterLooper); // particle has to be a photon and it has to have a e+ or e- as mother!
143147
float leadingAmplitude = g.amplitude()[0];
148+
info.purity = leadingAmplitude;
144149
if (std::abs(mcCluster.pdgCode()) == PDG_t::kElectron) {
145150
info.isLepton = true;
146151
info.convMotherId = getMotherIndexFromChain(mcCluster, mcClusterLooper, PDG_t::kGamma);
@@ -178,18 +183,18 @@ struct EmcalPhotonMcTask {
178183

179184
static constexpr std::array<const char*, static_cast<size_t>(TruthClass::NClasses)> kTruthClassNames = {
180185
"Conversion", "PhotonPairSamePi0", "PhotonPairDiffPi0", "PhotonPairOnePi0",
181-
"PhotonElectronSamePi0", "PhotonElectronDiffPi0", "PhotonElectronOnePi0", "PhotonElectronBS",
186+
"PhotonElectronSamePi0", "PhotonElectronDiffPi0", "PhotonElectronOnePi0", "BSPhotonElectron",
182187
"ElectronPairSamePi0", "ElectronPairDiffPi0", "ElectronPairOnePi0",
183188
"SplitPhotonCluster", "SplitLeptonCluster", "Background"};
184189

185-
Produces<aod::ConvTagCandidates> convTagCandidates;
190+
Produces<aod::ConvTagCandidates_001> convTagCandidates;
186191

187192
Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
188193
Configurable<std::string> grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"};
189194
Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
190195
Configurable<bool> skipGRPOquery{"skipGRPOquery", true, "skip grpo query"};
191196
Configurable<bool> writeTable{"writeTable", true, "write table for ML."};
192-
Configurable<std::vector<int>> classPrescale{"classPrescale", {1, 1, 700, 25, 1, 350, 15, 1, 1, 35, 2, 1, 1, 1000}, "prescale factor per TruthClass, indexed 0..10 matching the enum order"};
197+
Configurable<std::vector<int>> classPrescale{"classPrescale", {1, 1, 700, 25, 1, 350, 15, 1, 1, 35, 2, 1, 1, 1, 1, 1, 1000}, "prescale factor per TruthClass, indexed 0..10 matching the enum order"};
193198
Configurable<uint32_t> bkgPrescaleSeed{"bkgPrescaleSeed", 42, "seed for the background-prescale RNG"};
194199

195200
// configurable axis
@@ -334,6 +339,7 @@ struct EmcalPhotonMcTask {
334339
o2::aod::pwgem::photonmeson::utils::eventhistogram::addEventHistograms(&registry);
335340

336341
const AxisSpec thnAxisERec{thnConfigAxisPt, "#it{E}_{Rec} (GeV)"};
342+
const AxisSpec thnAxisPtRec{thnConfigAxisPt, "#it{p}_{T} (GeV/#it{c})"};
337343
const AxisSpec thnAxisInvMass{thnConfigAxisInvMass, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"};
338344

339345
const AxisSpec thnAxisrConvRec{100, 0, 500, "#it{R}_{rec}"};
@@ -354,7 +360,7 @@ struct EmcalPhotonMcTask {
354360
thnAxisCentOrMult = {thnConfigAxisMult, "FT0C Multiplicity"};
355361
}
356362

357-
auto hTruthLabel = registry.add<TH1>("hTruthLabel", "Truth label distribution;;Counts", HistType::kTH1D, {{static_cast<int>(TruthClass::NClasses), -0.5, static_cast<double>(TruthClass::NClasses) - 0.5}});
363+
auto hTruthLabel = registry.add<TH2>("hTruthLabel", "Truth label distribution;;Counts", HistType::kTH2D, {{static_cast<int>(TruthClass::NClasses), -0.5, static_cast<double>(TruthClass::NClasses) - 0.5}, thnAxisPtRec});
358364

359365
// set bin labels once at init, so histogram is human-readable without decoding the enum
360366
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::Conversion) + 1, "Conversion");
@@ -364,12 +370,15 @@ struct EmcalPhotonMcTask {
364370
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::PhotonElectronSamePi0) + 1, "PhotonElectronSamePi0");
365371
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::PhotonElectronDiffPi0) + 1, "PhotonElectronDiffPi0");
366372
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::PhotonElectronOnePi0) + 1, "PhotonElectronOnePi0");
367-
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::PhotonElectronBS) + 1, "PhotonElectronBS");
373+
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::BSPhotonElectron) + 1, "BSPhotonElectron");
368374
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::ElectronPairSamePi0) + 1, "ElectronPairSamePi0");
369375
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::ElectronPairDiffPi0) + 1, "ElectronPairDiffPi0");
370376
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::ElectronPairOnePi0) + 1, "ElectronPairOnePi0");
371377
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::SplitPhotonCluster) + 1, "SplitPhotonCluster");
372378
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::SplitLeptonCluster) + 1, "SplitLeptonCluster");
379+
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::PhotonBSPhotonPair) + 1, "PhotonBSPhotonPair");
380+
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::ElectronBSPhotonPair) + 1, "ElectronBSPhotonPair");
381+
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::BSPhotonPair) + 1, "BSPhotonPair");
373382
hTruthLabel->GetXaxis()->SetBinLabel(static_cast<int>(TruthClass::Background) + 1, "Background");
374383

375384
auto hPi0BothResolvedLost = registry.add<TH1>("EMCal/hPi0BothResolvedLost", "Confusion matrix for conversion tagging", HistType::kTH1D, {{2, -0.5, 1.5}});
@@ -582,7 +591,7 @@ struct EmcalPhotonMcTask {
582591
bool areConversionLegs = false;
583592
bool areSplitPhotonCluster = false;
584593
bool areSplitLeptonCluster = false;
585-
bool arePhotonElectronBS = false;
594+
bool areBSPhotonElectron = false;
586595

587596
auto c1 = classifyCluster(g1, mcCluster1, mcClusterLooper, mcClusterLooper2, mcParticles);
588597
auto c2 = classifyCluster(g2, mcCluster2, mcClusterLooper, mcClusterLooper2, mcParticles);
@@ -599,6 +608,9 @@ struct EmcalPhotonMcTask {
599608
}
600609
}
601610

611+
const bool isAnyBSPhoton = c1.isFromBremsstrahlung || c2.isFromBremsstrahlung;
612+
const bool areBSPhotons = c1.isFromBremsstrahlung && c2.isFromBremsstrahlung;
613+
602614
// if they are not a split cluster check for proper conversion pair
603615
if (!isSameDominantParticle && c1.isFromConv && c2.isFromConv && c1.convMotherId == c2.convMotherId) {
604616
emcFlagsFromTrueConversion.set(g1.globalIndex());
@@ -616,7 +628,7 @@ struct EmcalPhotonMcTask {
616628
// bremsstrahlung: one side is a photon born from the other side's lepton lineage
617629
const bool photonIsBS = (c1.isPhoton && c1.isFromBremsstrahlung) || (c2.isPhoton && c2.isFromBremsstrahlung);
618630
if (photonIsBS && ((c1.isLepton && c2.isPhoton) || (c2.isLepton && c1.isPhoton))) {
619-
arePhotonElectronBS = true;
631+
areBSPhotonElectron = true;
620632
} else {
621633
areFromSamePi0 = true;
622634
emcFlagsFromTrueMesonSameGamma.set(g1.globalIndex());
@@ -637,8 +649,14 @@ struct EmcalPhotonMcTask {
637649
bTruthLabel = static_cast<int8_t>(TruthClass::SplitLeptonCluster);
638650
} else if (areConversionLegs) {
639651
bTruthLabel = static_cast<int8_t>(TruthClass::Conversion);
640-
} else if (arePhotonElectronBS) {
641-
bTruthLabel = static_cast<int8_t>(TruthClass::PhotonElectronBS);
652+
} else if (areBSPhotonElectron) {
653+
bTruthLabel = static_cast<int8_t>(TruthClass::BSPhotonElectron);
654+
} else if (areBSPhotons && (c1.isFromPi0 || c2.isFromPi0)) {
655+
bTruthLabel = static_cast<int8_t>(TruthClass::BSPhotonPair);
656+
} else if (isAnyBSPhoton && (c1.isFromPi0 || c2.isFromPi0) && ((c1.isPhoton && c2.isLepton) || (c2.isPhoton && c1.isLepton))) {
657+
bTruthLabel = static_cast<int8_t>(TruthClass::ElectronBSPhotonPair);
658+
} else if (isAnyBSPhoton && (c1.isFromPi0 || c2.isFromPi0) && (c1.isPhoton && c2.isPhoton)) {
659+
bTruthLabel = static_cast<int8_t>(TruthClass::PhotonBSPhotonPair);
642660
} else if (areFromSamePi0) {
643661
if ((c1.isLepton && c2.isPhoton) || (c2.isLepton && c1.isPhoton)) {
644662
bTruthLabel = static_cast<int8_t>(TruthClass::PhotonElectronSamePi0);
@@ -665,13 +683,13 @@ struct EmcalPhotonMcTask {
665683
}
666684
}
667685

668-
registry.fill(HIST("hTruthLabel"), bTruthLabel);
686+
registry.fill(HIST("hTruthLabel"), bTruthLabel, vMeson.Pt());
669687

670688
// final tree values plus filling
671689
const int prescale = classPrescale.value[static_cast<uint>(bTruthLabel)];
672690
const bool keepThisRow = (prescale <= 1) || (std::uniform_int_distribution<int>(0, prescale - 1)(mRandGen) == 0);
673691
if (writeTable.value && keepThisRow) {
674-
convTagCandidates(collision.globalIndex(), vMeson.M(), harmonicET, deltaEta, deltaPhi, phiV, g1.e(), g2.e(), g1.m02(), g2.m02(), g1.time(), g2.time(), g1.nCells(), g2.nCells(), bTruthLabel, centOrMult);
692+
convTagCandidates(collision.globalIndex(), vMeson.M(), harmonicET, deltaEta, deltaPhi, phiV, g1.e(), g2.e(), g1.m02(), g2.m02(), g1.time(), g2.time(), g1.nCells(), g2.nCells(), c1.purity, c2.purity, bTruthLabel, centOrMult);
675693
}
676694
} // pair loop
677695
} // collision loop

0 commit comments

Comments
 (0)