Skip to content

Commit f3f2b3d

Browse files
author
Francesco Mazzaschi
committed
add process function to look for pe+e- findable candidates
1 parent cb509c0 commit f3f2b3d

1 file changed

Lines changed: 263 additions & 0 deletions

File tree

‎PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx‎

Lines changed: 263 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@
3636
#include <ReconstructionDataFormats/Track.h>
3737

3838
#include <TF1.h>
39+
#include <TMCProcess.h>
3940
#include <TPDGCode.h>
4041
#include <TVector3.h>
4142

@@ -184,6 +185,7 @@ struct Sigmaplusbuilder {
184185
const AxisSpec axisDiscriminant{200, -1., 1., "discriminant (GeV^{4}/#it{c}^{4})"};
185186
const AxisSpec axisMassSigma{200, 1.0, 1.4, "m_{p#gamma#gamma} (GeV/#it{c}^{2})"};
186187
const AxisSpec axisSigmaPt{100, 0., 6., "#it{p}_{T,#Sigma^{+}} (GeV/#it{c})"};
188+
const AxisSpec axisMomentum{100, 0., 10., "#it{p} (GeV/#it{c})"};
187189

188190
histos.add("Candidate/hSelectionCounter", "Candidate/hSelectionCounter", kTH1F, {axisCandSel});
189191
auto hCandSel = histos.get<TH1>(HIST("Candidate/hSelectionCounter"));
@@ -271,6 +273,44 @@ struct Sigmaplusbuilder {
271273
hPhotonTruthQA->GetXaxis()->SetBinLabel(4, "Gamma's mother is pi0");
272274
hPhotonTruthQA->GetXaxis()->SetBinLabel(5, "Pi0's mother is Sigma+");
273275
}
276+
277+
if (doprocessFindable) {
278+
const AxisSpec axisDetectorPresence{3, -0.5, 2.5, "detector"};
279+
const AxisSpec axisDuplicateTrack{2, -0.5, 1.5, "track"};
280+
const AxisSpec axisV0Presence{2, -0.5, 1.5, "V0 match"};
281+
const AxisSpec axisTPCClusters{160, -0.5, 159.5, "TPC clusters"};
282+
const AxisSpec axisPhotonMomResolution{200, -1., 1., "(#it{p}_{reco,#gamma} - #it{p}_{MC,#gamma})/#it{p}_{MC,#gamma}"};
283+
const AxisSpec axisPhotonCosPA{200, -1., 1., "cosPA_{#gamma}"};
284+
const AxisSpec axisPairMass{200, 0., 0.1, "m_{e^{+}e^{-}} (GeV/#it{c}^{2})"};
285+
histos.add("Findable/hSigmaPlusPt", "Findable/hSigmaPlusPt", kTH1F, {axisSigmaPt});
286+
histos.add("Findable/h2ProtonPtVsSigmaPlusPt", "Findable/h2ProtonPtVsSigmaPlusPt", kTH2F, {axisSigmaPt, axisMomentum});
287+
histos.add("Findable/hElectronPt", "Findable/hElectronPt", kTH1F, {axisMomentum});
288+
histos.add("Findable/hPositronPt", "Findable/hPositronPt", kTH1F, {axisMomentum});
289+
histos.add("Findable/hElectronPositronMass", "Findable/hElectronPositronMass", kTH1F, {axisPairMass});
290+
histos.add("Findable/hPhotonConversionRadius", "Findable/hPhotonConversionRadius", kTH1F, {axisPhotonRadius});
291+
histos.add("Findable/hPhotonMomentumResolution", "Findable/hPhotonMomentumResolution", kTH1F, {axisPhotonMomResolution});
292+
histos.add("Findable/hPhotonCosPA", "Findable/hPhotonCosPA", kTH1F, {axisPhotonCosPA});
293+
histos.add("Findable/hConversionPairV0Presence", "Findable/hConversionPairV0Presence", kTH1F, {axisV0Presence});
294+
histos.add("Findable/hDuplicateConversionTrackCounter", "Findable/hDuplicateConversionTrackCounter", kTH1F, {axisDuplicateTrack});
295+
histos.add("Findable/hDuplicateElectronTPCNClsFound", "Findable/hDuplicateElectronTPCNClsFound", kTH1F, {axisTPCClusters});
296+
histos.add("Findable/hDuplicatePositronTPCNClsFound", "Findable/hDuplicatePositronTPCNClsFound", kTH1F, {axisTPCClusters});
297+
histos.add("Findable/hElectronDetectorPresence", "Findable/hElectronDetectorPresence", kTH1F, {axisDetectorPresence});
298+
histos.add("Findable/hPositronDetectorPresence", "Findable/hPositronDetectorPresence", kTH1F, {axisDetectorPresence});
299+
auto hConversionPairV0Presence = histos.get<TH1>(HIST("Findable/hConversionPairV0Presence"));
300+
hConversionPairV0Presence->GetXaxis()->SetBinLabel(1, "valid pair");
301+
hConversionPairV0Presence->GetXaxis()->SetBinLabel(2, "in V0");
302+
auto hDuplicateConversionTrackCounter = histos.get<TH1>(HIST("Findable/hDuplicateConversionTrackCounter"));
303+
hDuplicateConversionTrackCounter->GetXaxis()->SetBinLabel(1, "e^{-}");
304+
hDuplicateConversionTrackCounter->GetXaxis()->SetBinLabel(2, "e^{+}");
305+
auto hElectronDetectorPresence = histos.get<TH1>(HIST("Findable/hElectronDetectorPresence"));
306+
auto hPositronDetectorPresence = histos.get<TH1>(HIST("Findable/hPositronDetectorPresence"));
307+
hElectronDetectorPresence->GetXaxis()->SetBinLabel(1, "ITS");
308+
hElectronDetectorPresence->GetXaxis()->SetBinLabel(2, "TPC");
309+
hElectronDetectorPresence->GetXaxis()->SetBinLabel(3, "TOF");
310+
hPositronDetectorPresence->GetXaxis()->SetBinLabel(1, "ITS");
311+
hPositronDetectorPresence->GetXaxis()->SetBinLabel(2, "TPC");
312+
hPositronDetectorPresence->GetXaxis()->SetBinLabel(3, "TOF");
313+
}
274314
}
275315

276316
// photon (PCM) candidate selection
@@ -466,6 +506,61 @@ struct Sigmaplusbuilder {
466506
return hasProton && hasPi0;
467507
}
468508

509+
template <bool IsElectron, typename TTrack>
510+
void fillFindableTrackDetectors(const TTrack& track)
511+
{
512+
auto fillDetector = [&](int detector) {
513+
if constexpr (IsElectron) {
514+
histos.fill(HIST("Findable/hElectronDetectorPresence"), detector);
515+
} else {
516+
histos.fill(HIST("Findable/hPositronDetectorPresence"), detector);
517+
}
518+
};
519+
520+
if (track.hasITS()) {
521+
fillDetector(0);
522+
}
523+
if (track.hasTPC()) {
524+
fillDetector(1);
525+
}
526+
if (track.hasTOF()) {
527+
fillDetector(2);
528+
}
529+
}
530+
531+
template <typename TMcPart>
532+
int findSigmaPlusMotherOfConversionElectron(const TMcPart& mcElectron, int& gammaIndex, std::array<float, 3>& conversionVertex)
533+
{
534+
if (std::abs(mcElectron.pdgCode()) != PDG_t::kElectron || mcElectron.getProcess() != TMCProcess::kPPair) {
535+
return -1;
536+
}
537+
538+
auto const& gammaMothers = mcElectron.template mothers_as<aod::McParticles>();
539+
if (gammaMothers.empty() || gammaMothers.front().pdgCode() != PDG_t::kGamma) {
540+
return -1;
541+
}
542+
543+
auto mcGamma = gammaMothers.front();
544+
if (mcGamma.getProcess() != TMCProcess::kPDecay) {
545+
return -1;
546+
}
547+
548+
auto const& pi0Mothers = mcGamma.template mothers_as<aod::McParticles>();
549+
if (pi0Mothers.empty() || std::abs(pi0Mothers.front().pdgCode()) != PDG_t::kPi0) {
550+
return -1;
551+
}
552+
553+
auto mcPi0 = pi0Mothers.front();
554+
auto const& sigmaMothers = mcPi0.template mothers_as<aod::McParticles>();
555+
if (sigmaMothers.empty() || std::abs(sigmaMothers.front().pdgCode()) != PDG_t::kSigmaPlus) {
556+
return -1;
557+
}
558+
559+
gammaIndex = mcGamma.globalIndex();
560+
conversionVertex = {mcElectron.vx(), mcElectron.vy(), mcElectron.vz()};
561+
return sigmaMothers.front().globalIndex();
562+
}
563+
469564
// Build a Sigma+ -> p pi0 candidate from a proton track and a PCM photon
470565
template <bool IsMC, typename TTracks, typename TTrack, typename TV0>
471566
void buildSigmaPlusCandidate(const TTrack& protonTrack, const TV0& photon, const std::array<float, 3>& pv)
@@ -903,6 +998,174 @@ struct Sigmaplusbuilder {
903998
}
904999
}
9051000
PROCESS_SWITCH(Sigmaplusbuilder, processMc, "Process MC", false);
1001+
1002+
void processFindable(CollisionsFullMC const& collisions, aod::V0Datas const& v0s, TracksFullMC const& tracks, aod::McParticles const&)
1003+
{
1004+
for (const auto& collision : collisions) {
1005+
auto tracksThisCollision = tracks.sliceBy(tracksPerCollisionMC, collision.globalIndex());
1006+
auto v0sThisCollision = v0s.sliceBy(v0PerCollision, collision.globalIndex());
1007+
1008+
for (const auto& protonTrack : tracksThisCollision) {
1009+
if (!protonTrack.has_mcParticle()) {
1010+
continue;
1011+
}
1012+
1013+
auto mcProton = protonTrack.template mcParticle_as<aod::McParticles>();
1014+
int sigmaIndex = findSigmaPlusMotherOfProton(mcProton);
1015+
if (sigmaIndex < 0) {
1016+
continue;
1017+
}
1018+
1019+
auto const& protonMothers = mcProton.template mothers_as<aod::McParticles>();
1020+
if (protonMothers.empty()) {
1021+
continue;
1022+
}
1023+
auto mcSigmaPlus = protonMothers.front();
1024+
if (std::abs(mcSigmaPlus.y()) > 1) {
1025+
continue;
1026+
}
1027+
histos.fill(HIST("Findable/h2ProtonPtVsSigmaPlusPt"), mcSigmaPlus.pt(), protonTrack.pt());
1028+
std::vector<int> seenElectronMcIds;
1029+
std::vector<int> seenElectronTrackIds;
1030+
std::vector<int> seenPositronMcIds;
1031+
std::vector<int> seenPositronTrackIds;
1032+
1033+
for (const auto& electronTrack : tracksThisCollision) {
1034+
if (!electronTrack.has_mcParticle()) {
1035+
continue;
1036+
}
1037+
1038+
if (electronTrack.tpcNClsFound() < 90 || electronTrack.sign() > 0) {
1039+
continue;
1040+
}
1041+
1042+
auto mcElectron = electronTrack.template mcParticle_as<aod::McParticles>();
1043+
if (mcElectron.pdgCode() != PDG_t::kElectron) {
1044+
continue;
1045+
}
1046+
1047+
int electronGammaIndex = -1;
1048+
std::array<float, 3> electronVertex{};
1049+
if (findSigmaPlusMotherOfConversionElectron(mcElectron, electronGammaIndex, electronVertex) != sigmaIndex) {
1050+
continue;
1051+
}
1052+
1053+
for (const auto& positronTrack : tracksThisCollision) {
1054+
if (!positronTrack.has_mcParticle()) {
1055+
continue;
1056+
}
1057+
1058+
if (positronTrack.tpcNClsFound() < 90 || positronTrack.sign() < 0) {
1059+
continue;
1060+
}
1061+
1062+
auto mcPositron = positronTrack.template mcParticle_as<aod::McParticles>();
1063+
if (mcPositron.pdgCode() != PDG_t::kPositron) {
1064+
continue;
1065+
}
1066+
1067+
int positronGammaIndex = -1;
1068+
std::array<float, 3> positronVertex{};
1069+
if (findSigmaPlusMotherOfConversionElectron(mcPositron, positronGammaIndex, positronVertex) != sigmaIndex || positronGammaIndex != electronGammaIndex) {
1070+
continue;
1071+
}
1072+
1073+
bool sameConversionPoint = electronVertex[0] == positronVertex[0] &&
1074+
electronVertex[1] == positronVertex[1] &&
1075+
electronVertex[2] == positronVertex[2];
1076+
if (!sameConversionPoint) {
1077+
continue;
1078+
}
1079+
1080+
int electronMcId = mcElectron.globalIndex();
1081+
int electronTrackId = electronTrack.globalIndex();
1082+
int positronMcId = mcPositron.globalIndex();
1083+
int positronTrackId = positronTrack.globalIndex();
1084+
1085+
bool seenElectronMc = false;
1086+
bool duplicateElectronTrack = false;
1087+
for (int iSeen = 0; iSeen < static_cast<int>(seenElectronMcIds.size()); ++iSeen) {
1088+
if (seenElectronMcIds[iSeen] == electronMcId) {
1089+
seenElectronMc = true;
1090+
duplicateElectronTrack |= (seenElectronTrackIds[iSeen] != electronTrackId);
1091+
}
1092+
}
1093+
if (duplicateElectronTrack) {
1094+
histos.fill(HIST("Findable/hDuplicateConversionTrackCounter"), 0);
1095+
histos.fill(HIST("Findable/hDuplicateElectronTPCNClsFound"), electronTrack.tpcNClsFound());
1096+
}
1097+
if (!seenElectronMc) {
1098+
seenElectronMcIds.push_back(electronMcId);
1099+
seenElectronTrackIds.push_back(electronTrackId);
1100+
}
1101+
1102+
bool seenPositronMc = false;
1103+
bool duplicatePositronTrack = false;
1104+
for (int iSeen = 0; iSeen < static_cast<int>(seenPositronMcIds.size()); ++iSeen) {
1105+
if (seenPositronMcIds[iSeen] == positronMcId) {
1106+
seenPositronMc = true;
1107+
duplicatePositronTrack |= (seenPositronTrackIds[iSeen] != positronTrackId);
1108+
}
1109+
}
1110+
if (duplicatePositronTrack) {
1111+
histos.fill(HIST("Findable/hDuplicateConversionTrackCounter"), 1);
1112+
histos.fill(HIST("Findable/hDuplicatePositronTPCNClsFound"), positronTrack.tpcNClsFound());
1113+
}
1114+
if (!seenPositronMc) {
1115+
seenPositronMcIds.push_back(positronMcId);
1116+
seenPositronTrackIds.push_back(positronTrackId);
1117+
}
1118+
if (duplicateElectronTrack || duplicatePositronTrack) {
1119+
continue;
1120+
}
1121+
1122+
std::array<float, 3> recoGammaMom{electronTrack.px() + positronTrack.px(), electronTrack.py() + positronTrack.py(), electronTrack.pz() + positronTrack.pz()};
1123+
float recoGammaP = std::sqrt(dot3(recoGammaMom, recoGammaMom));
1124+
constexpr float electronMass = o2::constants::physics::MassElectron;
1125+
float electronP2 = electronTrack.px() * electronTrack.px() + electronTrack.py() * electronTrack.py() + electronTrack.pz() * electronTrack.pz();
1126+
float positronP2 = positronTrack.px() * positronTrack.px() + positronTrack.py() * positronTrack.py() + positronTrack.pz() * positronTrack.pz();
1127+
float pairEnergy = std::sqrt(electronP2 + electronMass * electronMass) + std::sqrt(positronP2 + electronMass * electronMass);
1128+
float pairMass2 = pairEnergy * pairEnergy - recoGammaP * recoGammaP;
1129+
histos.fill(HIST("Findable/hElectronPositronMass"), std::sqrt(std::max(pairMass2, 0.f)));
1130+
1131+
auto const& gammaMothers = mcElectron.template mothers_as<aod::McParticles>();
1132+
if (!gammaMothers.empty()) {
1133+
auto mcGamma = gammaMothers.front();
1134+
std::array<float, 3> trueGammaMom{mcGamma.px(), mcGamma.py(), mcGamma.pz()};
1135+
float trueGammaP = std::sqrt(dot3(trueGammaMom, trueGammaMom));
1136+
if (trueGammaP > 0.f) {
1137+
histos.fill(HIST("Findable/hPhotonMomentumResolution"), (recoGammaP - trueGammaP) / trueGammaP);
1138+
}
1139+
}
1140+
1141+
std::array<float, 3> photonFlightVec{electronVertex[0] - collision.posX(), electronVertex[1] - collision.posY(), electronVertex[2] - collision.posZ()};
1142+
float flightNorm = std::sqrt(dot3(photonFlightVec, photonFlightVec));
1143+
if (flightNorm > 0.f && recoGammaP > 0.f) {
1144+
histos.fill(HIST("Findable/hPhotonCosPA"), dot3(photonFlightVec, recoGammaMom) / (flightNorm * recoGammaP));
1145+
}
1146+
1147+
histos.fill(HIST("Findable/hSigmaPlusPt"), mcSigmaPlus.pt());
1148+
histos.fill(HIST("Findable/hConversionPairV0Presence"), 0);
1149+
for (const auto& v0 : v0sThisCollision) {
1150+
auto posTrack = v0.template posTrack_as<TracksFullMC>();
1151+
auto negTrack = v0.template negTrack_as<TracksFullMC>();
1152+
bool sameChargeMatched = posTrack.globalIndex() == positronTrack.globalIndex() && negTrack.globalIndex() == electronTrack.globalIndex();
1153+
if (sameChargeMatched) {
1154+
histos.fill(HIST("Findable/hConversionPairV0Presence"), 1);
1155+
break;
1156+
}
1157+
}
1158+
histos.fill(HIST("Findable/hElectronPt"), electronTrack.pt());
1159+
histos.fill(HIST("Findable/hPositronPt"), positronTrack.pt());
1160+
histos.fill(HIST("Findable/hPhotonConversionRadius"), std::hypot(electronVertex[0], electronVertex[1]));
1161+
fillFindableTrackDetectors<true>(electronTrack);
1162+
fillFindableTrackDetectors<false>(positronTrack);
1163+
}
1164+
}
1165+
}
1166+
}
1167+
}
1168+
PROCESS_SWITCH(Sigmaplusbuilder, processFindable, "Process findable MC", false);
9061169
};
9071170

9081171
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)