From f3f2b3d66d3be528b8a2a05644109d1370a2bf08 Mon Sep 17 00:00:00 2001 From: Francesco Mazzaschi Date: Fri, 7 Aug 2026 16:50:20 +0200 Subject: [PATCH] add process function to look for pe+e- findable candidates --- .../Strangeness/sigmaplusbuilder.cxx | 263 ++++++++++++++++++ 1 file changed, 263 insertions(+) diff --git a/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx b/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx index cbd28655b99..739e7859c0f 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx @@ -36,6 +36,7 @@ #include #include +#include #include #include @@ -184,6 +185,7 @@ struct Sigmaplusbuilder { const AxisSpec axisDiscriminant{200, -1., 1., "discriminant (GeV^{4}/#it{c}^{4})"}; const AxisSpec axisMassSigma{200, 1.0, 1.4, "m_{p#gamma#gamma} (GeV/#it{c}^{2})"}; const AxisSpec axisSigmaPt{100, 0., 6., "#it{p}_{T,#Sigma^{+}} (GeV/#it{c})"}; + const AxisSpec axisMomentum{100, 0., 10., "#it{p} (GeV/#it{c})"}; histos.add("Candidate/hSelectionCounter", "Candidate/hSelectionCounter", kTH1F, {axisCandSel}); auto hCandSel = histos.get(HIST("Candidate/hSelectionCounter")); @@ -271,6 +273,44 @@ struct Sigmaplusbuilder { hPhotonTruthQA->GetXaxis()->SetBinLabel(4, "Gamma's mother is pi0"); hPhotonTruthQA->GetXaxis()->SetBinLabel(5, "Pi0's mother is Sigma+"); } + + if (doprocessFindable) { + const AxisSpec axisDetectorPresence{3, -0.5, 2.5, "detector"}; + const AxisSpec axisDuplicateTrack{2, -0.5, 1.5, "track"}; + const AxisSpec axisV0Presence{2, -0.5, 1.5, "V0 match"}; + const AxisSpec axisTPCClusters{160, -0.5, 159.5, "TPC clusters"}; + const AxisSpec axisPhotonMomResolution{200, -1., 1., "(#it{p}_{reco,#gamma} - #it{p}_{MC,#gamma})/#it{p}_{MC,#gamma}"}; + const AxisSpec axisPhotonCosPA{200, -1., 1., "cosPA_{#gamma}"}; + const AxisSpec axisPairMass{200, 0., 0.1, "m_{e^{+}e^{-}} (GeV/#it{c}^{2})"}; + histos.add("Findable/hSigmaPlusPt", "Findable/hSigmaPlusPt", kTH1F, {axisSigmaPt}); + histos.add("Findable/h2ProtonPtVsSigmaPlusPt", "Findable/h2ProtonPtVsSigmaPlusPt", kTH2F, {axisSigmaPt, axisMomentum}); + histos.add("Findable/hElectronPt", "Findable/hElectronPt", kTH1F, {axisMomentum}); + histos.add("Findable/hPositronPt", "Findable/hPositronPt", kTH1F, {axisMomentum}); + histos.add("Findable/hElectronPositronMass", "Findable/hElectronPositronMass", kTH1F, {axisPairMass}); + histos.add("Findable/hPhotonConversionRadius", "Findable/hPhotonConversionRadius", kTH1F, {axisPhotonRadius}); + histos.add("Findable/hPhotonMomentumResolution", "Findable/hPhotonMomentumResolution", kTH1F, {axisPhotonMomResolution}); + histos.add("Findable/hPhotonCosPA", "Findable/hPhotonCosPA", kTH1F, {axisPhotonCosPA}); + histos.add("Findable/hConversionPairV0Presence", "Findable/hConversionPairV0Presence", kTH1F, {axisV0Presence}); + histos.add("Findable/hDuplicateConversionTrackCounter", "Findable/hDuplicateConversionTrackCounter", kTH1F, {axisDuplicateTrack}); + histos.add("Findable/hDuplicateElectronTPCNClsFound", "Findable/hDuplicateElectronTPCNClsFound", kTH1F, {axisTPCClusters}); + histos.add("Findable/hDuplicatePositronTPCNClsFound", "Findable/hDuplicatePositronTPCNClsFound", kTH1F, {axisTPCClusters}); + histos.add("Findable/hElectronDetectorPresence", "Findable/hElectronDetectorPresence", kTH1F, {axisDetectorPresence}); + histos.add("Findable/hPositronDetectorPresence", "Findable/hPositronDetectorPresence", kTH1F, {axisDetectorPresence}); + auto hConversionPairV0Presence = histos.get(HIST("Findable/hConversionPairV0Presence")); + hConversionPairV0Presence->GetXaxis()->SetBinLabel(1, "valid pair"); + hConversionPairV0Presence->GetXaxis()->SetBinLabel(2, "in V0"); + auto hDuplicateConversionTrackCounter = histos.get(HIST("Findable/hDuplicateConversionTrackCounter")); + hDuplicateConversionTrackCounter->GetXaxis()->SetBinLabel(1, "e^{-}"); + hDuplicateConversionTrackCounter->GetXaxis()->SetBinLabel(2, "e^{+}"); + auto hElectronDetectorPresence = histos.get(HIST("Findable/hElectronDetectorPresence")); + auto hPositronDetectorPresence = histos.get(HIST("Findable/hPositronDetectorPresence")); + hElectronDetectorPresence->GetXaxis()->SetBinLabel(1, "ITS"); + hElectronDetectorPresence->GetXaxis()->SetBinLabel(2, "TPC"); + hElectronDetectorPresence->GetXaxis()->SetBinLabel(3, "TOF"); + hPositronDetectorPresence->GetXaxis()->SetBinLabel(1, "ITS"); + hPositronDetectorPresence->GetXaxis()->SetBinLabel(2, "TPC"); + hPositronDetectorPresence->GetXaxis()->SetBinLabel(3, "TOF"); + } } // photon (PCM) candidate selection @@ -466,6 +506,61 @@ struct Sigmaplusbuilder { return hasProton && hasPi0; } + template + void fillFindableTrackDetectors(const TTrack& track) + { + auto fillDetector = [&](int detector) { + if constexpr (IsElectron) { + histos.fill(HIST("Findable/hElectronDetectorPresence"), detector); + } else { + histos.fill(HIST("Findable/hPositronDetectorPresence"), detector); + } + }; + + if (track.hasITS()) { + fillDetector(0); + } + if (track.hasTPC()) { + fillDetector(1); + } + if (track.hasTOF()) { + fillDetector(2); + } + } + + template + int findSigmaPlusMotherOfConversionElectron(const TMcPart& mcElectron, int& gammaIndex, std::array& conversionVertex) + { + if (std::abs(mcElectron.pdgCode()) != PDG_t::kElectron || mcElectron.getProcess() != TMCProcess::kPPair) { + return -1; + } + + auto const& gammaMothers = mcElectron.template mothers_as(); + if (gammaMothers.empty() || gammaMothers.front().pdgCode() != PDG_t::kGamma) { + return -1; + } + + auto mcGamma = gammaMothers.front(); + if (mcGamma.getProcess() != TMCProcess::kPDecay) { + return -1; + } + + auto const& pi0Mothers = mcGamma.template mothers_as(); + if (pi0Mothers.empty() || std::abs(pi0Mothers.front().pdgCode()) != PDG_t::kPi0) { + return -1; + } + + auto mcPi0 = pi0Mothers.front(); + auto const& sigmaMothers = mcPi0.template mothers_as(); + if (sigmaMothers.empty() || std::abs(sigmaMothers.front().pdgCode()) != PDG_t::kSigmaPlus) { + return -1; + } + + gammaIndex = mcGamma.globalIndex(); + conversionVertex = {mcElectron.vx(), mcElectron.vy(), mcElectron.vz()}; + return sigmaMothers.front().globalIndex(); + } + // Build a Sigma+ -> p pi0 candidate from a proton track and a PCM photon template void buildSigmaPlusCandidate(const TTrack& protonTrack, const TV0& photon, const std::array& pv) @@ -903,6 +998,174 @@ struct Sigmaplusbuilder { } } PROCESS_SWITCH(Sigmaplusbuilder, processMc, "Process MC", false); + + void processFindable(CollisionsFullMC const& collisions, aod::V0Datas const& v0s, TracksFullMC const& tracks, aod::McParticles const&) + { + for (const auto& collision : collisions) { + auto tracksThisCollision = tracks.sliceBy(tracksPerCollisionMC, collision.globalIndex()); + auto v0sThisCollision = v0s.sliceBy(v0PerCollision, collision.globalIndex()); + + for (const auto& protonTrack : tracksThisCollision) { + if (!protonTrack.has_mcParticle()) { + continue; + } + + auto mcProton = protonTrack.template mcParticle_as(); + int sigmaIndex = findSigmaPlusMotherOfProton(mcProton); + if (sigmaIndex < 0) { + continue; + } + + auto const& protonMothers = mcProton.template mothers_as(); + if (protonMothers.empty()) { + continue; + } + auto mcSigmaPlus = protonMothers.front(); + if (std::abs(mcSigmaPlus.y()) > 1) { + continue; + } + histos.fill(HIST("Findable/h2ProtonPtVsSigmaPlusPt"), mcSigmaPlus.pt(), protonTrack.pt()); + std::vector seenElectronMcIds; + std::vector seenElectronTrackIds; + std::vector seenPositronMcIds; + std::vector seenPositronTrackIds; + + for (const auto& electronTrack : tracksThisCollision) { + if (!electronTrack.has_mcParticle()) { + continue; + } + + if (electronTrack.tpcNClsFound() < 90 || electronTrack.sign() > 0) { + continue; + } + + auto mcElectron = electronTrack.template mcParticle_as(); + if (mcElectron.pdgCode() != PDG_t::kElectron) { + continue; + } + + int electronGammaIndex = -1; + std::array electronVertex{}; + if (findSigmaPlusMotherOfConversionElectron(mcElectron, electronGammaIndex, electronVertex) != sigmaIndex) { + continue; + } + + for (const auto& positronTrack : tracksThisCollision) { + if (!positronTrack.has_mcParticle()) { + continue; + } + + if (positronTrack.tpcNClsFound() < 90 || positronTrack.sign() < 0) { + continue; + } + + auto mcPositron = positronTrack.template mcParticle_as(); + if (mcPositron.pdgCode() != PDG_t::kPositron) { + continue; + } + + int positronGammaIndex = -1; + std::array positronVertex{}; + if (findSigmaPlusMotherOfConversionElectron(mcPositron, positronGammaIndex, positronVertex) != sigmaIndex || positronGammaIndex != electronGammaIndex) { + continue; + } + + bool sameConversionPoint = electronVertex[0] == positronVertex[0] && + electronVertex[1] == positronVertex[1] && + electronVertex[2] == positronVertex[2]; + if (!sameConversionPoint) { + continue; + } + + int electronMcId = mcElectron.globalIndex(); + int electronTrackId = electronTrack.globalIndex(); + int positronMcId = mcPositron.globalIndex(); + int positronTrackId = positronTrack.globalIndex(); + + bool seenElectronMc = false; + bool duplicateElectronTrack = false; + for (int iSeen = 0; iSeen < static_cast(seenElectronMcIds.size()); ++iSeen) { + if (seenElectronMcIds[iSeen] == electronMcId) { + seenElectronMc = true; + duplicateElectronTrack |= (seenElectronTrackIds[iSeen] != electronTrackId); + } + } + if (duplicateElectronTrack) { + histos.fill(HIST("Findable/hDuplicateConversionTrackCounter"), 0); + histos.fill(HIST("Findable/hDuplicateElectronTPCNClsFound"), electronTrack.tpcNClsFound()); + } + if (!seenElectronMc) { + seenElectronMcIds.push_back(electronMcId); + seenElectronTrackIds.push_back(electronTrackId); + } + + bool seenPositronMc = false; + bool duplicatePositronTrack = false; + for (int iSeen = 0; iSeen < static_cast(seenPositronMcIds.size()); ++iSeen) { + if (seenPositronMcIds[iSeen] == positronMcId) { + seenPositronMc = true; + duplicatePositronTrack |= (seenPositronTrackIds[iSeen] != positronTrackId); + } + } + if (duplicatePositronTrack) { + histos.fill(HIST("Findable/hDuplicateConversionTrackCounter"), 1); + histos.fill(HIST("Findable/hDuplicatePositronTPCNClsFound"), positronTrack.tpcNClsFound()); + } + if (!seenPositronMc) { + seenPositronMcIds.push_back(positronMcId); + seenPositronTrackIds.push_back(positronTrackId); + } + if (duplicateElectronTrack || duplicatePositronTrack) { + continue; + } + + std::array recoGammaMom{electronTrack.px() + positronTrack.px(), electronTrack.py() + positronTrack.py(), electronTrack.pz() + positronTrack.pz()}; + float recoGammaP = std::sqrt(dot3(recoGammaMom, recoGammaMom)); + constexpr float electronMass = o2::constants::physics::MassElectron; + float electronP2 = electronTrack.px() * electronTrack.px() + electronTrack.py() * electronTrack.py() + electronTrack.pz() * electronTrack.pz(); + float positronP2 = positronTrack.px() * positronTrack.px() + positronTrack.py() * positronTrack.py() + positronTrack.pz() * positronTrack.pz(); + float pairEnergy = std::sqrt(electronP2 + electronMass * electronMass) + std::sqrt(positronP2 + electronMass * electronMass); + float pairMass2 = pairEnergy * pairEnergy - recoGammaP * recoGammaP; + histos.fill(HIST("Findable/hElectronPositronMass"), std::sqrt(std::max(pairMass2, 0.f))); + + auto const& gammaMothers = mcElectron.template mothers_as(); + if (!gammaMothers.empty()) { + auto mcGamma = gammaMothers.front(); + std::array trueGammaMom{mcGamma.px(), mcGamma.py(), mcGamma.pz()}; + float trueGammaP = std::sqrt(dot3(trueGammaMom, trueGammaMom)); + if (trueGammaP > 0.f) { + histos.fill(HIST("Findable/hPhotonMomentumResolution"), (recoGammaP - trueGammaP) / trueGammaP); + } + } + + std::array photonFlightVec{electronVertex[0] - collision.posX(), electronVertex[1] - collision.posY(), electronVertex[2] - collision.posZ()}; + float flightNorm = std::sqrt(dot3(photonFlightVec, photonFlightVec)); + if (flightNorm > 0.f && recoGammaP > 0.f) { + histos.fill(HIST("Findable/hPhotonCosPA"), dot3(photonFlightVec, recoGammaMom) / (flightNorm * recoGammaP)); + } + + histos.fill(HIST("Findable/hSigmaPlusPt"), mcSigmaPlus.pt()); + histos.fill(HIST("Findable/hConversionPairV0Presence"), 0); + for (const auto& v0 : v0sThisCollision) { + auto posTrack = v0.template posTrack_as(); + auto negTrack = v0.template negTrack_as(); + bool sameChargeMatched = posTrack.globalIndex() == positronTrack.globalIndex() && negTrack.globalIndex() == electronTrack.globalIndex(); + if (sameChargeMatched) { + histos.fill(HIST("Findable/hConversionPairV0Presence"), 1); + break; + } + } + histos.fill(HIST("Findable/hElectronPt"), electronTrack.pt()); + histos.fill(HIST("Findable/hPositronPt"), positronTrack.pt()); + histos.fill(HIST("Findable/hPhotonConversionRadius"), std::hypot(electronVertex[0], electronVertex[1])); + fillFindableTrackDetectors(electronTrack); + fillFindableTrackDetectors(positronTrack); + } + } + } + } + } + PROCESS_SWITCH(Sigmaplusbuilder, processFindable, "Process findable MC", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)