From 2b5f4faa90d76ad2e3646880584589527087f64b Mon Sep 17 00:00:00 2001 From: Stefanie Mrozinski Date: Thu, 12 Nov 2026 11:23:32 +0100 Subject: [PATCH 1/2] Add new bremsstrahlungstemplate histograms + centrality pool in mixer --- .../TableProducer/associateMCinfoPhoton.cxx | 152 ++++++++++++++- .../photonMixedEventProducer.cxx | 23 ++- .../TableProducer/photonconversionbuilder.h | 103 ++++++++++ PWGEM/PhotonMeson/Tasks/photonhbt.cxx | 182 +++++++++++++++++- 4 files changed, 448 insertions(+), 12 deletions(-) diff --git a/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx b/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx index 7248b009e04..2d4051d2a47 100644 --- a/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx +++ b/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx @@ -25,11 +25,11 @@ #include #include #include -#include #include #include #include #include +#include #include #include @@ -48,6 +48,7 @@ #include #include #include +#include #include using namespace o2; @@ -97,6 +98,10 @@ struct AssociateMCInfoPhoton { Configurable useIsSelected{"useIsSelected", true, "flag to enable the usage of isSelected from our EM event selection. Disabling this will increase the derived data size!"}; Configurable useIsEOI{"useIsEOI", true, "flag to enable the usage of the EOI table"}; + Configurable fillBremsChain{"fillBremsChain", false, "MC truth histograms of the bremsstrahlung chain: photon conversion, radiation by a leg, conversion of the radiated photon (BremsChain/)"}; + Configurable bremsChainMaxR{"bremsChainMaxR", 90.f, "max conversion radius (cm) for the bremsstrahlung chain histograms"}; + Configurable bremsChainMinX{"bremsChainMinX", 0.05f, "min energy fraction E(radiated photon) / E(parent photon) for the bremsstrahlung chain histograms"}; + HistogramRegistry registry{"EMMCEvent"}; void init(o2::framework::InitContext&) @@ -133,6 +138,19 @@ struct AssociateMCInfoPhoton { registry.add(Form("Generated/h2PtY_Accepted_%s", ParticleNames[i].data()), Form("Generated %s", ParticleNames[i].data()), kTH2F, {axisPt, axisRapidity}, true); } + if (fillBremsChain) { + const AxisSpec axR{45, 0.f, 90.f, "#it{R}_{conv}(A) (cm)"}; + const AxisSpec axPhi{360, 0.f, 360.f, "#varphi_{conv}(A) (deg)"}; + const AxisSpec axPtA{{0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 10.0}, "#it{p}_{T}(A) (GeV/#it{c})"}; + const AxisSpec axX{{0.05, 0.1, 0.3, 0.5, 1.0}, "#it{E}(B) / #it{E}(A)"}; + registry.add("BremsChain/hGenPhiPt", "generated primary photons, |#eta| < cfg_max_eta_gamma;#varphi of the momentum (deg);#it{p}_{T} (GeV/#it{c})", kTH2F, {{360, 0.f, 360.f}, axPtA}, false); + registry.add("BremsChain/hConvA", "true conversions of primary photons", kTH3F, {axR, axPhi, axPtA}, false); + registry.add("BremsChain/hConvARec", "true conversions of primary photons, photon reconstructed", kTH3F, {axR, axPhi, axPtA}, false); + registry.add("BremsChain/hChain", "photons radiated by a leg of A, cumulative stages", kTHnSparseF, {axR, axPhi, axPtA, axX, {4, -0.5f, 3.5f, "stage: 0 radiated (in the MC stack), 1 B converts, 2 B reconstructed, 3 B and A reconstructed"}}, false); + registry.add("BremsChain/hGeometry", "radiated photons that convert: where, relative to A", kTHnSparseF, {{18, 0.f, 90.f, "#it{R}_{conv}(A) (cm)"}, {40, -5.f, 95.f, "#it{R}_{conv}(B) #minus #it{R}_{conv}(A) (cm)"}, {40, -20.f, 20.f, "#varphi_{conv}(B) #minus #varphi_{conv}(A) (deg)"}, axX, {3, 0.5f, 3.5f, "stage: 1 B converts, 2 B reconstructed, 3 B and A reconstructed"}}, false); + registry.add("BremsChain/hEmission", "radiated photons with #it{E}(B) / #it{E}(A) > 0.1: emission point;#it{R}_{conv}(A) (cm);#it{R}_{emission} #minus #it{R}_{conv}(A) (cm)", kTH2F, {axR, {50, 0.f, 100.f}}, false); + } + // reserve space for generated vectors if that process enabled auto hBinFinder = registry.get(HIST("Generated/h2PtY_Pi0")); LOGF(info, "Binned generated processing enabled. Initialising with %i elements...", hBinFinder->GetNcells()); @@ -230,6 +248,108 @@ struct AssociateMCInfoPhoton { } // end of mother chain loop } + // conversion point of a photon: production vertex of its e+ e- daughters; false if the photon did not convert into e+ e- (in the MC stack) + template + bool conversionPoint(TMCParticle const& photon, TMCParticles const& mcParticles, float& vx, float& vy, int64_t& idEle, int64_t& idPos) + { + idEle = -1; + idPos = -1; + if (!photon.has_daughters()) { + return false; + } + for (int id = photon.daughtersIds()[0]; id <= photon.daughtersIds()[1]; ++id) { + if (id < 0 || id >= mcParticles.size()) { + continue; + } + const auto d = mcParticles.rawIteratorAt(id); + if (d.pdgCode() == PDG_t::kElectron && idEle < 0) { + idEle = id; + vx = d.vx(); + vy = d.vy(); + } else if (d.pdgCode() == PDG_t::kPositron && idPos < 0) { + idPos = id; + } + } + return idEle >= 0 && idPos >= 0; + } + + // MC truth of the bremsstrahlung chain for one MC collision (see the histogram definitions in init) + template + void fillBremsChainTruth(TGroupedMCParticles const& groupedMcParticles, TMCParticles const& mcParticles, std::unordered_set const& recoPhotons) + { + static constexpr float RadToDeg = 57.29578f; + static constexpr float FullTurn = 360.f, HalfTurn = 180.f; + auto phiDeg = [](float y, float x) { + const float phi = std::atan2(y, x) * RadToDeg; + return phi < 0.f ? phi + FullTurn : phi; + }; + for (const auto& photon : groupedMcParticles) { + if (photon.pdgCode() != PDG_t::kGamma || !(photon.isPhysicalPrimary() || photon.producedByGenerator()) || std::fabs(photon.eta()) > cfg_max_eta_photon) { + continue; + } + registry.fill(HIST("BremsChain/hGenPhiPt"), phiDeg(photon.py(), photon.px()), photon.pt()); + float vxA = 0.f, vyA = 0.f; + int64_t idEle = -1, idPos = -1; + if (!conversionPoint(photon, mcParticles, vxA, vyA, idEle, idPos)) { + continue; + } + const float rA = std::hypot(vxA, vyA); + if (rA > bremsChainMaxR) { + continue; + } + const float phiA = phiDeg(vyA, vxA); + const bool recA = recoPhotons.contains(photon.globalIndex()); + registry.fill(HIST("BremsChain/hConvA"), rA, phiA, photon.pt()); + if (recA) { + registry.fill(HIST("BremsChain/hConvARec"), rA, phiA, photon.pt()); + } + for (const auto& idLeg : {idEle, idPos}) { + const auto leg = mcParticles.rawIteratorAt(idLeg); + if (!leg.has_daughters()) { + continue; + } + for (int id = leg.daughtersIds()[0]; id <= leg.daughtersIds()[1]; ++id) { + if (id < 0 || id >= mcParticles.size()) { + continue; + } + const auto radiated = mcParticles.rawIteratorAt(id); + if (radiated.pdgCode() != PDG_t::kGamma) { + continue; + } + const float x = std::min(radiated.e() / std::max(photon.e(), 1.e-6f), 0.9999f); + if (x < bremsChainMinX) { + continue; + } + registry.fill(HIST("BremsChain/hChain"), rA, phiA, photon.pt(), x, 0.f); + if (x > 0.1f) { // o2-linter: disable=magic-number (energy fraction of the calculation) + registry.fill(HIST("BremsChain/hEmission"), rA, std::hypot(radiated.vx(), radiated.vy()) - rA); + } + float vxB = 0.f, vyB = 0.f; + int64_t idEleB = -1, idPosB = -1; + if (!conversionPoint(radiated, mcParticles, vxB, vyB, idEleB, idPosB)) { + continue; + } + const float rB = std::hypot(vxB, vyB); + if (rB > bremsChainMaxR) { + continue; + } + float dPhi = phiDeg(vyB, vxB) - phiA; + if (dPhi > HalfTurn) { + dPhi -= FullTurn; + } else if (dPhi < -HalfTurn) { + dPhi += FullTurn; + } + const bool recB = recoPhotons.contains(radiated.globalIndex()); + const int lastStage = recB ? (recA ? 3 : 2) : 1; + for (int stage = 1; stage <= lastStage; ++stage) { + registry.fill(HIST("BremsChain/hChain"), rA, phiA, photon.pt(), x, static_cast(stage)); + registry.fill(HIST("BremsChain/hGeometry"), rA, std::clamp(rB - rA, -4.99f, 94.99f), std::clamp(dPhi, -19.99f, 19.99f), x, static_cast(stage)); + } + } + } + } + } + template void skimmingMC(MyCollisionsMC const& collisions, aod::BCs const&, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, TTracks const& o2tracks, TFwdTracks const&, TPCMs const& v0photons, TPCMLegs const& legs, TPHOSs const&, TEMCs const& emcclusters, TEMPrimaryElectrons const& emprimaryelectrons) { @@ -258,6 +378,33 @@ struct AssociateMCInfoPhoton { // mc collision iter auto mcCollisionIter = mcCollisions.begin(); + // MC photons reconstructed as a stored V0 photon (both legs point to the same photon), for the bremsstrahlung chain histograms + std::unordered_set recoPhotons; + if constexpr (static_cast(system & kPCM) && !std::is_same_v) { + if (fillBremsChain) { + auto legIter = legs.begin(); + auto trackIter = o2tracks.begin(); + for (const auto& v0 : v0photons) { + std::array mother{-1, -2}; + int i = 0; + for (const auto& legId : {v0.posTrackId(), v0.negTrackId()}) { + legIter.setCursor(legId); + trackIter.setCursor(legIter.trackId()); + if (trackIter.has_mcParticle()) { + mcParticleIter.setCursor(trackIter.mcParticleId()); + if (mcParticleIter.has_mothers()) { + mother[i] = mcParticleIter.mothersIds()[0]; + } + } + ++i; + } + if (mother[0] == mother[1] && mother[0] >= 0) { + recoPhotons.insert(mother[0]); + } + } + } + } + for (const auto& collision : collisions) { registry.fill(HIST("hEventCounter"), 1); @@ -320,6 +467,9 @@ struct AssociateMCInfoPhoton { // make an entry for this MC event only if it was not already added to the table if (!fEventLabels.contains(mcCollisionIter.globalIndex())) { + if (fillBremsChain) { + fillBremsChainTruth(groupedMcParticles, mcParticles, recoPhotons); + } mcevents(mcCollisionIter.globalIndex(), mcCollisionIter.generatorsID(), mcCollisionIter.posX(), mcCollisionIter.posY(), mcCollisionIter.posZ(), mcCollisionIter.impactParameter(), mcCollisionIter.eventPlaneAngle()); fEventLabels[mcCollisionIter.globalIndex()] = fCounter.events; fCounter.events++; diff --git a/PWGEM/PhotonMeson/TableProducer/photonMixedEventProducer.cxx b/PWGEM/PhotonMeson/TableProducer/photonMixedEventProducer.cxx index 5e873e46102..07b14054672 100644 --- a/PWGEM/PhotonMeson/TableProducer/photonMixedEventProducer.cxx +++ b/PWGEM/PhotonMeson/TableProducer/photonMixedEventProducer.cxx @@ -52,6 +52,7 @@ #include #include #include +#include #include #include #include @@ -122,6 +123,7 @@ struct PhotonMixedEventProducer { Configurable cfgCentEstimator{"cfgCentEstimator", 2, "centrality estimator: 0 FT0M, 1 FT0A, 2 FT0C"}; Configurable cfgCentMin{"cfgCentMin", 10.f, "min centrality (one class)"}; Configurable cfgCentMax{"cfgCentMax", 20.f, "max centrality"}; + Configurable cfgCentBinWidth{"cfgCentBinWidth", 0.f, "pool centrality bin width (%) inside [cfgCentMin, cfgCentMax); 0 = one bin. Partner events are then matched in centrality, which removes the Delta E tilt of the mixed event (same-event pairs share the spectrum of their event)"}; Configurable cfgZMax{"cfgZMax", 10.f, "max |z_vtx| (cm)"}; Configurable cfgZBinWidth{"cfgZBinWidth", 2.f, "pool z bin width (cm)"}; Configurable cfgOccupancyEstimator{"cfgOccupancyEstimator", 0, "occupancy: 0 trackOccupancyInTimeRange, 1 ft0cOccupancyInTimeRange"}; @@ -139,7 +141,7 @@ struct PhotonMixedEventProducer { HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; - std::map, std::deque> mPools; + std::map, std::deque> mPools; // (z bin, occupancy bin, centrality bin) std::mt19937 mRandom; int mRunNumber = -1; int64_t mEventUid = 0; @@ -173,6 +175,7 @@ struct PhotonMixedEventProducer { registry.add("Event/hOccupancy", "accepted events;occupancy", kTH1D, {{200, 0.f, 20000.f}}); registry.add("Pool/hDepthVsZ", "pool size when used;z_{vtx} (cm);events in the pool", kTH2D, {{nZ, -cfgZMax, cfgZMax}, {cfgPoolDepth + 1, -0.5f, cfgPoolDepth + 0.5f}}); registry.add("Pool/hDepthVsOccBin", "pool size when used;occupancy bin;events in the pool", kTH2D, {{static_cast(cfgOccupancyBins->size()) - 1, -0.5f, cfgOccupancyBins->size() - 1.5f}, {cfgPoolDepth + 1, -0.5f, cfgPoolDepth + 0.5f}}); + registry.add("Pool/hDepthVsCentBin", "pool size when used;centrality bin;events in the pool", kTH2D, {{std::max(1, nCentBins()), -0.5f, std::max(1, nCentBins()) - 0.5f}, {cfgPoolDepth + 1, -0.5f, cfgPoolDepth + 0.5f}}); registry.add("Pool/hPartnersVsN", "partner events available vs needed;N of the current event;events in the pool", kTH2D, {axN, {cfgPoolDepth + 1, -0.5f, cfgPoolDepth + 0.5f}}); registry.add("Pool/hNPooledPhotons", "photons in all pools (memory), per event;photons", kTH1D, {{500, 0.f, 5.e5f}}); registry.add("Pool/hNPools", "filled pools, per event;pools", kTH1D, {{1000, 0.f, 10000.f}}); @@ -224,6 +227,20 @@ struct PhotonMixedEventProducer { return cfgOccupancyEstimator.value == 1 ? collision.ft0cOccupancyInTimeRange() : static_cast(collision.trackOccupancyInTimeRange()); } + int nCentBins() const + { + if (cfgCentBinWidth.value <= 0.f) { + return 1; + } + return std::max(1, static_cast(std::ceil((cfgCentMax.value - cfgCentMin.value) / cfgCentBinWidth.value - 1e-4f))); // o2-linter: disable=magic-number (rounding guard) + } + int centralityBin(const float cent) const + { + if (cfgCentBinWidth.value <= 0.f) { + return 0; + } + return std::clamp(static_cast(std::floor((cent - cfgCentMin.value) / cfgCentBinWidth.value)), 0, nCentBins() - 1); + } int occupancyBin(const float occ) const { const auto& edges = cfgOccupancyBins.value; @@ -590,8 +607,10 @@ struct PhotonMixedEventProducer { const int nPhotons = static_cast(current.photons.size()); registry.fill(HIST("Event/hNPhotons"), nPhotons); - auto& pool = mPools[{zBin, occBin}]; + const int centBin = centralityBin(cent); + auto& pool = mPools[{zBin, occBin, centBin}]; const int depth = static_cast(pool.size()); + registry.fill(HIST("Pool/hDepthVsCentBin"), centBin, depth); registry.fill(HIST("Pool/hPartnersVsN"), nPhotons, depth); registry.fill(HIST("Pool/hDepthVsZ"), collision.posZ(), depth); registry.fill(HIST("Pool/hDepthVsOccBin"), occBin, depth); diff --git a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h index 3c77efa0efd..931ef42e9d4 100644 --- a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h +++ b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h @@ -373,6 +373,7 @@ struct PhotonConversionBuilder { Configurable freeLegRuleMask{"freeLegRuleMask", 5, "with freeLegRun1Rules: which run-1 rules a crossed pairing with a pooled free leg gets: bit 0 AP cut, bit 1 purity gate (pairings with a leg without ITS), bit 2 finder fit options (side-based drift correction, finder PCA; pairings with a leg without ITS); 7 = all run-1 rules, 0 = the lenient crossed fit"}; Configurable freeLegMinPurity{"freeLegMinPurity", -999.f, "with freeLegRun1Rules and bit 1 of freeLegRuleMask: purity threshold for crossed pairings with a pooled free leg; below -998 recovery.finderMinScore (as run 1)"}; Configurable mixedDiagnostics{"mixedDiagnostics", false, "MC diagnostics of run 2 (hCrossedFit, four-leg decisions, matched configurations, mixed pair classes, fake classes, hQSym, swap/brems pointing, alternative pairings): CPU on data, no effect on the stored photons"}; + Configurable fillBremsJoint{"fillBremsJoint", false, "pairs of stored photons after run 2, same and mixed artificial events, data and MC: Mixed/Pairs/hBremsJoint = q_inv x |Delta E| x k_T x centrality x bremsstrahlung pointing (distance category incl. no allowed point, angle category, vertex-anchored score) x MC class (15 without labels); Mixed/Pairs/hBremsJointNorm (every pair below 500 MeV/c without the pointing, for the mixed-event normalisation); Mixed/Pairs/hBremsAngle (measured angle of the pointing against the momenta and extrapolation lengths that set its resolution); Mixed/hNStoredVsCent (stored photons per event vs centrality). For the data-driven bremsstrahlung template: every pair is filled once, no selection"}; Configurable diagFreeLegFits{"diagFreeLegFits", false, "diagnostic: every crossed pairing with a pooled free leg is fitted twice as a trial, with the lenient crossed options and with the run-1 rules, and Mixed/Decision/hFreeLegFit records both results (accepted or the rejection line), split by MC truth; the trial candidates are removed again, the real fit follows the configuration (CPU and QA side effect: V0 QA histograms count the trial fits)"}; Configurable freeLegRun1Rules{"freeLegRun1Rules", true, "crossed pairings of a pooled free leg (sourceV0Id -2) with a stored leg are fitted as run 1 would have fitted them: with a leg without ITS as a finder candidate (side-based drift correction, AP cut, purity gate), with both legs with ITS as an SVertexer V0 (AP cut, no purity gate); stored x stored pairings keep crossedRelaxAP (overlay mixed event)"}; Configurable freeFreeFinderOnly{"freeFreeFinderOnly", false, "with skipFreeFreeCrossed off: a crossed pairing of two pooled free legs is formed only if at least one leg has no ITS (the run-1 finder domain); free x free with both legs with ITS could only have been an SVertexer V0 in run 1 and is skipped (overlay mixed event)"}; @@ -799,6 +800,17 @@ struct PhotonConversionBuilder { registry.add("Mixed/Pairs/hQClean", "photon pairs, all legs TPC shared fraction < 0.4;artificial event;0 before / 1 after run 2;q_{inv} (GeV/c)", kTH3F, {axType, {2, -0.5f, 1.5f}, {500, 0.f, 0.5f}}); registry.add("Mixed/Pairs/hBrems", "photon pairs, bremsstrahlung pointing;artificial event;0 before / 1 after run 2;class (see mixFakeClass);MC relation (0 none, 1 radiated by a leg, 2 siblings, 3 leg daughter of a leg);q_{inv} (GeV/c);log_{10} d_{min} (cm);log_{10} angle to the leg (rad);radius of the closest point (cm);B converts further out;p_{T} of B (GeV/c);all legs TPC shared fraction < 0.4", kTHnSparseF, {axType, {2, -0.5f, 1.5f}, {18, -0.5f, 17.5f}, {4, -0.5f, 3.5f}, {50, 0.f, 0.1f}, {40, -3.f, 2.f}, {40, -4.f, 0.5f}, {45, 0.f, 90.f}, {2, -0.5f, 1.5f}, {20, 0.f, 1.f}, {2, -0.5f, 1.5f}}); registry.add("Mixed/Pairs/hSwapDecision", "full swaps after run 2 (classes 14 and 17);artificial event;0 created in run 2, 1 carried from run 1, 2 one of each;0 true pairings scored better, 1 swap scored better, 2 a true pairing failed the fit;q_{inv} (GeV/c);all legs TPC shared fraction < 0.4", kTHnSparseF, {axType, {3, -0.5f, 2.5f}, {3, -0.5f, 2.5f}, {50, 0.f, 0.1f}, {2, -0.5f, 1.5f}}); + // prototype 2026-11-11: vertex-anchored pointing score next to the d_min tag, on the same pairs (after run 2, q < 50 MeV/c), with the MC check + // of the matched leg + registry.add("Mixed/Pairs/hBremsVtx", "photon pairs after run 2, pointing: d_{min} tag vs vertex-anchored score;artificial event;class (see mixFakeClass);q_{inv} (GeV/c);d_{min}: 0 < 0.1, 1 0.1-0.3, 2 0.3-1 cm, 3 > 1 cm or no point;angle < 0.1 rad;log_{10} tangent miss at B's conversion point (cm);log_{10} angle of B to the line emission point - conversion point (rad);p_{T} of B (GeV/c);d_{min} match: 0 pair not radiated, 1 leg is the true emitter, 2 wrong leg or role;score match: same;score (#sigma)", kTHnSparseF, {axType, {18, -0.5f, 17.5f}, {5, 0.f, 0.05f}, {4, -0.5f, 3.5f}, {2, -0.5f, 1.5f}, {16, -3.f, 1.f}, {16, -4.f, 0.f}, {std::vector{0., 0.3, 0.6, 1., 10.}, "pT B"}, {3, -0.5f, 2.5f}, {3, -0.5f, 2.5f}, {20, 0.f, 10.f}, {std::vector{0., 2., 5., 10., 20., 50., 200.}, "lever arm emission point - conversion point of B (cm)"}}); + // data-driven bremsstrahlung template (mixed_group.fillBremsJoint): every pair of stored photons after run 2 below 100 MeV/c, once + const AxisSpec axCentJoint{std::vector{0., 10., 20., 30., 40., 50., 60., 70., 80., 100.}, "centrality FT0C (%)"}; + const AxisSpec axDEJoint{std::vector{0., 0.01, 0.02, 0.03, 0.05, 0.075, 0.1, 0.15, 0.2, 0.3, 0.5, 10.}, "|#Delta E| (GeV)"}; + const AxisSpec axKtJoint{std::vector{0., 0.3, 0.6, 1., 10.}, "k_{T} (GeV/c)"}; + registry.add("Mixed/Pairs/hBremsJointNorm", "photon pairs after run 2 below 500 MeV/c, no pointing (mixed-event normalisation of hBremsJoint);artificial event;class (see mixFakeClass, 15 without labels);q_{inv} (GeV/c);|#Delta E| (GeV);k_{T} (GeV/c);centrality FT0C (%)", kTHnSparseF, {axType, {18, -0.5f, 17.5f}, {100, 0.f, 0.5f}, axDEJoint, axKtJoint, axCentJoint}); + registry.add("Mixed/Pairs/hBremsAngle", "photon pairs after run 2 below 50 MeV/c with an allowed emission point at d_{min} < 1 cm: measured angle of the pointing and what sets its resolution;artificial event;class (see mixFakeClass, 15 without labels);q_{inv} (GeV/c);log_{10} angle of B to the leg tangent at the closest point (rad);p_{T} of B (GeV/c);p_{T} of the matched leg (GeV/c);closest point -> conversion point of B (cm);arc length from the leg's reference point to the closest point (cm)", kTHnSparseF, {axType, {18, -0.5f, 17.5f}, {10, 0.f, 0.05f}, {30, -3.f, 0.f}, {std::vector{0., 0.2, 0.4, 0.6, 1., 10.}, "pT B"}, {std::vector{0., 0.1, 0.2, 0.4, 1., 10.}, "pT leg"}, {std::vector{0., 2., 5., 10., 20., 50., 200.}, "lever"}, {std::vector{0., 5., 10., 20., 40., 100.}, "arc"}}); + registry.add("Mixed/Pairs/hBremsJoint", "photon pairs after run 2 (data and MC);artificial event;class (see mixFakeClass, 15 without labels);q_{inv} (GeV/c);|#Delta E| (GeV);k_{T} (GeV/c);centrality FT0C (%);d_{min}: 0 < 0.1, 1 0.1-0.3, 2 0.3-1, 3 > 1 cm, 4 no allowed point;angle: 0 < 0.03, 1 0.03-0.1, 2 > 0.1 rad;vertex-anchored score (#sigma), last bin: none", kTHnSparseF, {axType, {18, -0.5f, 17.5f}, {20, 0.f, 0.1f}, {std::vector{0., 0.01, 0.02, 0.03, 0.05, 0.075, 0.1, 0.15, 0.2, 0.3, 0.5, 10.}, "dE"}, {std::vector{0., 0.3, 0.6, 1., 10.}, "kT"}, axCentJoint, {5, -0.5f, 4.5f}, {3, -0.5f, 2.5f}, {std::vector{0., 1., 2., 3., 5., 10., 12.}, "score"}}); + registry.add("Mixed/hNStoredVsCent", "stored photons after run 2;artificial event;centrality FT0C (%);photons per event", kTH3F, {axType, {20, 0.f, 100.f}, {101, -0.5f, 100.5f}}); } void initCCDB(MyBCs::iterator const& bc) @@ -2875,11 +2887,16 @@ struct PhotonConversionBuilder { } struct BremsPointing { float dMin = 999.f, angle = 3.15f, rEmit = -1.f, outward = 0.f, ptB = 0.f; + int order = -1, leg = -1; // 0 = a radiates b, 1 = b radiates a; leg 0 positive, 1 negative + float ptLeg = 0.f, arc = 0.f, leverD = 0.f; + float score = 1.e9f, missT = 999.f, thetaB = 3.15f, ptBV = 0.f, lever = 0.f; + int orderV = -1, legV = -1; }; template BremsPointing mixBremsPointing(MixPairPhoton const& a, MixPairPhoton const& b, TLegs const& legs) const { constexpr float StepCm = 0.25f, MaxArcCm = 100.f, MarginCm = 2.f; + constexpr float MinLeverCm = 0.5f, SigB0 = 0.0066f, SigBPow = -0.3f, MinPtB = 0.1f, SigT0 = 0.004f, SigVtxCm = 0.4f; BremsPointing best; for (int order = 0; order < 2; ++order) { const auto& emitter = order == 0 ? a : b; @@ -2890,7 +2907,9 @@ struct PhotonConversionBuilder { } const std::array u{brems.px / pB, brems.py / pB, brems.pz / pB}; const float rA = std::hypot(emitter.vx, emitter.vy), rB = std::hypot(brems.vx, brems.vy); + int iLeg = -1; for (const int64_t row : {emitter.posRow, emitter.eleRow}) { + ++iLeg; const auto leg = legs.rawIteratorAt(row); const o2::track::TrackPar par(leg.x(), leg.alpha(), std::array{leg.y(), leg.z(), leg.snp(), leg.tgl(), leg.signed1Pt()}); o2::math_utils::CircleXYf_t circle; @@ -2927,6 +2946,31 @@ struct PhotonConversionBuilder { best.rEmit = rH; best.outward = rB > rA ? 1.f : 0.f; best.ptB = std::hypot(brems.px, brems.py); + best.order = order; + best.leg = iLeg; + best.ptLeg = std::fabs(leg.signed1Pt()) > 0.f ? 1.f / std::fabs(leg.signed1Pt()) : 0.f; + best.arc = std::fabs(s); + best.leverD = -t; + } + const float lever = std::hypot(d[0], d[1], d[2]); // emission point -> B's conversion point (d = point - conversion point) + if (lever > MinLeverCm) { + const std::array tanV{-sense * std::sin(psi) / norm, sense * std::cos(psi) / norm, tgl / norm}; + const std::array chord{-d[0] / lever, -d[1] / lever, -d[2] / lever}; + const float thT = std::acos(std::clamp(tanV[0] * chord[0] + tanV[1] * chord[1] + tanV[2] * chord[2], -1.f, 1.f)); + const float thB = std::acos(std::clamp(u[0] * chord[0] + u[1] * chord[1] + u[2] * chord[2], -1.f, 1.f)); + const float ptBrems = std::hypot(brems.px, brems.py); + const float sigB = std::hypot(SigB0 * std::pow(std::max(ptBrems, MinPtB), SigBPow), SigVtxCm / lever); + const float sigT = std::hypot(SigT0, SigVtxCm / lever); + const float sc = std::hypot(thT / sigT, thB / sigB); + if (sc < best.score) { + best.score = sc; + best.missT = lever * std::sin(thT); + best.thetaB = thB; + best.ptBV = ptBrems; + best.lever = lever; + best.orderV = order; + best.legV = iLeg; + } } } } @@ -3018,11 +3062,45 @@ struct PhotonConversionBuilder { return 15.f; } + // data-driven bremsstrahlung template: all pairs of stored photons of one artificial event after run 2 below 500 MeV/c (normalisation), below 100 MeV/c with the + template + void fillBremsJoint(const float type, const float cent, std::vector const& photons, TLegs const& legs) + { + constexpr float MaxQNorm = 0.5f, MaxQ = 0.1f, MaxQAngle = 0.05f, DVTight = 0.1f, DTight = 0.3f, DLoose = 1.f, DNone = 100.f, AngleTight = 0.03f, AngleLoose = 0.1f, ScoreNone = 11.f, Half = 0.5f, LogFloorA = 1e-3f; + for (size_t i = 0; i < photons.size(); ++i) { + for (size_t j = i + 1; j < photons.size(); ++j) { + const auto& a = photons[i]; + const auto& b = photons[j]; + if (skipSameSourcePair(type, a, b)) { + continue; + } + const float q = qInv(a.px, a.py, a.pz, b.px, b.py, b.pz); + if (q >= MaxQNorm) { + continue; + } + const float cls = mixPairClass(a, b); + const float dE = std::fabs(std::hypot(a.px, a.py, a.pz) - std::hypot(b.px, b.py, b.pz)); + const float kT = Half * std::hypot(a.px + b.px, a.py + b.py); + registry.fill(HIST("Mixed/Pairs/hBremsJointNorm"), type, cls, q, dE, kT, cent); + if (q >= MaxQ) { + continue; + } + const auto bp = mixBremsPointing(a, b, legs); + const float dCat = bp.dMin < DVTight ? 0.f : (bp.dMin < DTight ? 1.f : (bp.dMin < DLoose ? 2.f : (bp.dMin < DNone ? 3.f : 4.f))); + const float aCat = bp.angle < AngleTight ? 0.f : (bp.angle < AngleLoose ? 1.f : 2.f); + registry.fill(HIST("Mixed/Pairs/hBremsJoint"), type, cls, q, dE, kT, cent, dCat, aCat, std::min(bp.score, ScoreNone)); + if (q < MaxQAngle && bp.dMin < DLoose) { + registry.fill(HIST("Mixed/Pairs/hBremsAngle"), type, cls, q, std::log10(std::max(bp.angle, LogFloorA)), bp.ptB, bp.ptLeg, bp.leverD, bp.arc); + } + } + } + } template void fillMixSwapBrems(const float type, const float stage, const int64_t colId, std::vector const& photons, std::map, float> const& candScore, TLegs const& legs) { constexpr float MaxQ = 0.1f, MaxShared = 0.4f; + constexpr float MaxQVtx = 0.05f, DTagVTight = 0.1f, DTagTight = 0.3f, DTagLoose = 1.f, AngleTag = 0.1f, ScoreClip = 9.99f, StageAfter = 0.5f, LogFloorD = 1e-3f, LogFloorA = 1e-4f; auto score = [&](int64_t pos, int64_t ele) { const auto it = candScore.find({colId, pos, ele}); return it == candScore.end() ? -1.f : it->second; @@ -3041,6 +3119,25 @@ struct PhotonConversionBuilder { const float cls = mixPairClass(a, b); const float clean = mixPairGeometry(a, b)[2] < MaxShared ? 1.f : 0.f; const auto bp = mixBremsPointing(a, b, legs); + if (stage > StageAfter && q < MaxQVtx) { + auto match = [&](const int order, const int iLeg) { + if (order < 0) { + return 0.f; + } + const auto& emitter = order == 0 ? a : b; + const auto& brems = order == 0 ? b : a; + const auto& y = iLeg == 0 ? emitter.posMC : emitter.eleMC; + for (const auto& x : {&brems.posMC, &brems.eleMC}) { + if (x->evt == y.evt && x->gmid >= 0 && x->gmid == y.pid) { + return 1.f; + } + } + return mixRelKind(a, b) == 1 ? 2.f : 0.f; + }; + const float dCat = bp.dMin < DTagVTight ? 0.f : (bp.dMin < DTagTight ? 1.f : (bp.dMin < DTagLoose ? 2.f : 3.f)); + registry.fill(HIST("Mixed/Pairs/hBremsVtx"), type, cls, q, dCat, bp.angle < AngleTag ? 1.f : 0.f, std::log10(std::max(bp.missT, LogFloorD)), std::log10(std::max(bp.thetaB, LogFloorA)), + bp.ptBV, match(bp.order, bp.leg), match(bp.orderV, bp.legV), std::min(bp.score, ScoreClip), bp.lever); + } registry.fill(HIST("Mixed/Pairs/hBrems"), type, stage, cls, static_cast(mixRelKind(a, b)), q, std::log10(std::max(bp.dMin, 1e-3f)), std::log10(std::max(bp.angle, 1e-4f)), bp.rEmit, bp.outward, bp.ptB, clean); // o2-linter: disable=magic-number (log floors) if (stage < 0.5f || (cls != 14.f && cls != 17.f)) { // o2-linter: disable=magic-number (swap classes) @@ -3806,6 +3903,12 @@ struct PhotonConversionBuilder { registry.fill(HIST("Mixed/hNStoredVsIn"), type, nIn, n[0] + n[1]); registry.fill(HIST("Mixed/hNOriginalLost"), type, nIn - n[0]); registry.fill(HIST("Mixed/hNCrossedStored"), type, n[1]); + if (mixed.fillBremsJoint) { + registry.fill(HIST("Mixed/hNStoredVsCent"), type, collision.centFT0C(), n[0] + n[1]); + if (const auto itJ = after.find(collision.globalIndex()); itJ != after.end()) { + fillBremsJoint(type, collision.centFT0C(), itJ->second, legs); + } + } if (const auto itA = after.find(collision.globalIndex()); mixed.mixedDiagnostics && itA != after.end()) { fillMixPairs(type, itA->second); fillSameSourcePairs(type, itA->second); diff --git a/PWGEM/PhotonMeson/Tasks/photonhbt.cxx b/PWGEM/PhotonMeson/Tasks/photonhbt.cxx index 9c89b99bd94..5903fc5904b 100644 --- a/PWGEM/PhotonMeson/Tasks/photonhbt.cxx +++ b/PWGEM/PhotonMeson/Tasks/photonhbt.cxx @@ -38,7 +38,6 @@ #include #include #include -#include #include #include #include @@ -1328,6 +1327,8 @@ struct Photonhbt { fRegistryPairMC.add((spBase + "hEtaPt").c_str(), "#eta vs p_{T};#eta;p_{T} (GeV/c)", kTH2F, {axisEta, axisPt}, true); fRegistryPairMC.add((spBase + "hLegOriginMatrix").c_str(), "pos-leg origin vs neg-leg origin;pos leg;neg leg", kTH2F, {axisLegOrigin, axisLegOrigin}, true); } + fRegistryPairMC.add("SinglePhotonMC/Response/hResponse", "true reconstructed photons: response (MC)", kTHnSparseF,{{40, 0.f, 4.f, "#it{p}_{T}^{true} (GeV/#it{c})"},{60, 0.3f, 1.5f, "#it{E}^{rec} / #it{E}^{true}"}, {40, 0.f, 0.04f, "angle(rec, true) (rad)"},{18, 0.f, 90.f, "#it{R}_{conv} (cm)"},{21, -0.05f, 1.f, "hardest photon radiated by a leg: #it{E} / #it{E}^{true} (first bin: none)"},{2, -0.5f, 1.5f, "that radiated photon reconstructed in this event"},{3, -0.5f, 2.5f, "this photon radiated by a conversion e^{#pm}: 0 no, 1 yes and parent #gamma lost, 2 yes and parent reconstructed"},{21, -0.05f, 1.f, "#it{E}^{true} / #it{E}(parent #gamma) (first bin: n/a)"},{51, -0.005f, 0.25f, "angle(this, parent #gamma) MC (rad) (first bin: n/a)"}},true); + fRegistryPairMC.add("SinglePhotonMC/Response/hGenPt", "generated primary photons of the selected collisions, |#eta| < cfgMaxEtaV0;#it{p}_{T} (GeV/#it{c})", kTH1F, {{80, 0.f, 4.f}}, true); const AxisSpec axisFakeSubtype{3, -0.5f, 2.5f, "0=ordinary, 1=swap doublet, 2=partial swap"}; fRegistryPairMC.add("Pair/same/MC/hFakeSubtypeVsQinv", "fake pair subtype vs q_{inv};q_{inv} (GeV/c);subtype", kTH2D, {axisQinv, axisFakeSubtype}, true); @@ -2986,6 +2987,88 @@ struct Photonhbt { } } + + template + void fillSinglePhotonResponse(TCollision const& collision, TPhotons const& photonsColl, TMCParticles const& mcParticles, TCut const& cut) + { + if (collision.has_emmcevent()) { + auto gen = mcParticles.sliceBy(perMCCollisionEMMCParts, collision.emmceventId()); + for (const auto& m : gen) { + if (m.pdgCode() != kGamma || !(m.isPhysicalPrimary() || m.producedByGenerator()) || std::fabs(m.eta()) > pcmcuts.cfgMaxEtaV0.value) { + continue; + } + fRegistryPairMC.fill(HIST("SinglePhotonMC/Response/hGenPt"), m.pt()); + } + } + auto angleOf = [](float ax, float ay, float az, float bx, float by, float bz) { + const float n = std::sqrt((ax * ax + ay * ay + az * az) * (bx * bx + by * by + bz * bz)); + return n > 0.f ? std::acos(std::clamp((ax * bx + ay * by + az * bz) / n, -1.f, 1.f)) : 0.f; + }; + // MC photons reconstructed in this collision (true photons after the photon cuts) + std::unordered_set recoMothers; + for (const auto& g : photonsColl) { + if (!cut.template IsSelected(g)) { + continue; + } + const auto mcg = buildPhotonMCInfo(g, mcParticles); + if (mcg.hasMC && mcg.sameMother && mcg.isTruePhoton) { + recoMothers.insert(mcg.motherId); + } + } + for (const auto& g : photonsColl) { + if (!cut.template IsSelected(g)) { + continue; + } + const auto mcg = buildPhotonMCInfo(g, mcParticles); + if (!(mcg.hasMC && mcg.sameMother && mcg.isTruePhoton)) { + continue; + } + const auto mother = mcParticles.iteratorAt(mcg.motherId); + const float eTrue = std::hypot(mother.px(), mother.py(), mother.pz()); + if (eTrue <= 0.f) { + continue; + } + const float eRec = std::hypot(g.px(), g.py(), g.pz()); + const float angRec = std::min(angleOf(g.px(), g.py(), g.pz(), mother.px(), mother.py(), mother.pz()), 0.0399f); + // the hardest photon radiated by one of the two legs, and whether it was reconstructed here + float eRadFrac = -0.025f; + float radReco = 0.f; + for (const auto& legId : {mcg.mcPosId, mcg.mcNegId}) { + const auto leg = mcParticles.iteratorAt(legId); + if (!leg.has_daughters()) { + continue; + } + for (const auto& dId : leg.daughtersIds()) { + const auto dtr = mcParticles.iteratorAt(dId); + if (dtr.pdgCode() != kGamma) { + continue; + } + const float frac = std::hypot(dtr.px(), dtr.py(), dtr.pz()) / eTrue; + if (frac > eRadFrac) { + eRadFrac = std::min(frac, 0.999f); + radReco = recoMothers.contains(dId) ? 1.f : 0.f; + } + } + } + // this photon itself radiated by a conversion electron: mother photon <- electron <- parent photon + float isRad = 0.f, fracParent = -0.025f, angParent = -0.0025f; + if (mother.has_mothers()) { + const auto ele = mcParticles.iteratorAt(mother.mothersIds()[0]); + if (std::abs(ele.pdgCode()) == PDG_t::kElectron && ele.has_mothers()) { + const auto parent = mcParticles.iteratorAt(ele.mothersIds()[0]); + if (parent.pdgCode() == kGamma) { + isRad = recoMothers.contains(parent.globalIndex()) ? 2.f : 1.f; + const float eParent = std::hypot(parent.px(), parent.py(), parent.pz()); + fracParent = eParent > 0.f ? std::min(eTrue / eParent, 0.999f) : 0.999f; + angParent = std::min(angleOf(mother.px(), mother.py(), mother.pz(), parent.px(), parent.py(), parent.pz()), 0.249f); + } + } + } + fRegistryPairMC.fill(HIST("SinglePhotonMC/Response/hResponse"), mother.pt(), std::clamp(eRec / eTrue, 0.301f, 1.499f), angRec, std::min(std::hypot(g.vx(), g.vy()), 89.9f), + eRadFrac, radReco, isRad, fracParent, angParent); + } + } + /*************************************************/ // MAIN LOOP FOR PAIRING MC INFORMATION /*************************************************/ @@ -3047,6 +3130,7 @@ struct Photonhbt { } } } + fillSinglePhotonResponse(collision, photonsColl, mcParticles, cut); } std::unordered_set idsAfterPairCuts; for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photonsColl, photonsColl))) { @@ -4383,7 +4467,7 @@ struct Photonhbt { const AxisSpec axOrigin{4, -0.5f, 3.5f, "stored photons: 0 both SVertexer originals, 1 one finder original, 2 both finder, 3 a crossed pairing"}; const AxisSpec axClean{2, -0.5f, 1.5f, "all four legs below the TPC shared-cluster cut"}; fRegistryPairMC.add("PairClass/Pairs/hPairClass", "photon pairs by MC class", kTHnSparseF, {axSample, axClass, axQ, axOrigin, axClean}); - fRegistryPairMC.add("PairClass/Pairs/hDEQinv", "photon pairs: |#Delta E| vs q_{inv} by MC class", kTHnSparseF, {axSample, axClass, {100, 0.f, 0.5f, "#it{q}_{inv} (GeV/#it{c})"}, {100, 0.f, 0.5f, "|#it{E}_{1} #minus #it{E}_{2}| (GeV)"}, {std::vector{0., 0.2, 0.4, 0.6, 0.9, 1.5, 3.}, "#it{k}_{T} (GeV/#it{c})"}}); + fRegistryPairMC.add("PairClass/Pairs/hDEQinv", "photon pairs: |#Delta E| vs q_{inv} by MC class", kTHnSparseF, {axSample, axClass, {100, 0.f, 0.5f, "#it{q}_{inv} (GeV/#it{c})"}, {100, 0.f, 0.5f, "|#it{E}_{1} #minus #it{E}_{2}| (GeV)"}, {std::vector{0., 0.2, 0.4, 0.6, 0.9, 1.5, 3.}, "#it{k}_{T} (GeV/#it{c})"}, {std::vector{0., 0.02, 0.04, 0.3}, "|#it{q}_{out}| (LCMS) (GeV/#it{c})"}}); fRegistryPairMC.add("PairClass/Pairs/hFreeLegPair", "photon pairs by MC class: photons built with a pooled free leg", kTHnSparseF, {axSample, axClass, axQ, {3, -0.5f, 2.5f, "photons of the pair with a free leg: 0, 1, 2"}, {3, -0.5f, 2.5f, "crossed photons (run-2 pairings) in the pair: 0, 1, 2"}}); fRegistryPairMC.add("PairClass/Pairs/hHalfSwapLegs", "half swaps: run-1 origin of the two legs of the shared true photon", kTHnSparseF, {axSample, axQ, {4, -0.5f, 3.5f, "legs of the shared photon: 0 both free legs, 1 one stored photon (taken apart), 2 two different stored photons, 3 one free one stored"}, {3, -0.5f, 2.5f, "photons of the pair with a free leg: 0, 1, 2"}}); fRegistryPairMC.add("PairClass/Pairs/hKtQ", "photon pairs: k_{T} vs q_{inv} by MC class", kTHnSparseF, {axSample, axClass, axQ, {100, 0.f, 2.f, "#it{k}_{T} (GeV/#it{c})"}}); @@ -4395,16 +4479,18 @@ struct Photonhbt { fRegistryPairMC.add("PairClass/Pairs/hPairQ", "photon pairs: fitted q, four-leg invariant mass, true q by MC class", kTHnSparseF, {axSample, axClass, {100, 0.f, 0.2f, "#it{q}_{inv} of the fitted photons (GeV/#it{c})"}, {100, 0.f, 0.2f, "#it{M} of the four legs (GeV/#it{c}^{2})"}, {100, 0.f, 0.2f, "#it{q}_{inv} of the two true photons (GeV/#it{c})"}, {6, 0.f, 1.5f, "#it{k}_{T} (GeV/#it{c})"}}); fRegistryPairMC.add("PairClass/Pairs/hPairAlpha", "photon pairs: signed energy sharing of the two photons by MC class", kTHnSparseF, {axSample, axClass, {50, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {20, -1.f, 1.f, "#alpha_{1}"}, {20, -1.f, 1.f, "#alpha_{2}"}, {10, 0.f, 1.f, "|#it{E}_{1} #minus #it{E}_{2}| / (#it{E}_{1} + #it{E}_{2}) of the mother photons"}}); fRegistryPairMC.add("PairClass/Pairs/hPairLegs", "photon pairs: leg diagnostics by MC class", kTHnSparseF, {axSample, axClass, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {15, 0.f, 0.3f, "max over the two photons of |tan#lambda_{+} #minus tan#lambda_{#minus}|"}, {10, 0.f, 1.f, "softest leg #it{p} (GeV/#it{c})"}, {5, -0.5f, 4.5f, "legs with ITS (of 4)"}, axOrigin}); + fRegistryPairMC.add("PairClass/Pairs/hDupPieces", "duplicates (classes 1, 2, 3, 12): the two pieces of the shared MC particle", kTHnSparseF, {axSample, {4, -0.5f, 3.5f, "class: 0 = 1 (same photon twice), 1 = 2, 2 = 3, 3 = 12 (same track)"}, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {40, 0.f, 0.2f, "|tan#lambda_{1} #minus tan#lambda_{2}| of the two pieces"}, {20, 0.f, 1.f, "|#it{p}_{T,1} #minus #it{p}_{T,2}| / max(#it{p}_{T}) of the two pieces"}, {3, -0.5f, 2.5f, "charge slots: 0 same (both +/both #minus), 1 opposite, 2 no shared particle found"}, axOrigin, {3, -0.5f, 2.5f, "pieces with ITS (of 2)"}}); fRegistryPairMC.add("PairClass/Pairs/hPairCE", "photon pairs: central-electrode variables by MC class", kTHnSparseF, {axSample, axClass, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {std::vector{0., 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.6, 0.9, 2.}, "min |#eta| of the two photons"}, {std::vector{0., 1., 2., 5., 10., 20., 50., 300.}, "min |z_{conv}^{true}| (cm)"}, {5, -0.5f, 4.5f, "legs crossing the central electrode in the TPC (of 4)"}, {5, -0.5f, 4.5f, "legs with ITS (of 4)"}}); fRegistryPairMC.add("PairClass/Pairs/hPairR", "photon pairs: conversion radii by MC class", kTHnSparseF, {axSample, axClass, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {90, 0.f, 90.f, "#it{R}_{xy} of the inner photon (cm)"}, {90, 0.f, 90.f, "#it{R}_{xy} of the outer photon (cm)"}, {91, -1.f, 90.f, "true #it{R}_{xy} of the inner photon (cm), -1 without MC"}}); fRegistryPairMC.add("PairClass/Pairs/hBrems", "bremsstrahlung pairs (classes 16 and 5, MC)", kTHnSparseF, {axSample, {2, -0.5f, 1.5f, "0 class 16 (both true), 1 class 5 (with a fake)"}, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {5, -0.5f, 4.5f, "radiating electron from: 0 photon conversion, 1 #pi^{0} Dalitz, 2 #eta Dalitz, 3 other, 4 not found"}, {6, -0.5f, 5.5f, "mother of the parent photon: 0 #pi^{0}, 1 #eta, 2 #omega, 3 e^{#pm} (bremsstrahlung), 4 other, 5 none / not a photon"}, {50, 0.f, 0.25f, "#theta(#gamma_{brems}, radiating e^{#pm} rec.) (rad)"}, {50, 0.f, 0.25f, "#theta(#gamma_{brems}, parent #gamma) (rad), MC"}, {60, -40.f, 80.f, "true #it{R}_{xy}(#gamma_{brems}) #minus #it{R}_{xy}(parent) (cm)"}, {20, 0.f, 1.f, "#it{E}(#gamma_{brems}) / #it{E}(parent), MC"}}); fRegistryPairMC.add("PairClass/Pairs/hPhotonQuality", "photons of the pairs by MC class (each photon once per pair)", kTHnSparseF, {axSample, axClass, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {2, -0.5f, 1.5f, "0 photon radiated by an e^{#pm} (MC mother e^{#pm}), 1 other"}, {50, 0.f, 50.f, "KF #chi^{2}/ndf"}, {60, 0.f, 0.3f, "pointing angle to the PV (rad)"}, {50, 0.f, 5.f, "|DCA_{xy}| to the PV (cm)"}, {50, 0.f, 5.f, "|DCA_{z}| to the PV (cm)"}}); { - const char* anc = "0 #pi^{0}, 1 #eta, 2 #omega, 3 #eta', 4 #rho/K*, 5 K^{0}, 6 baryon, 7 none (primary), 8 other"; + const char* anc = "0 #pi^{0}, 1 #eta, 2 #omega, 3 #eta', 4 #rho/K*, 5 K^{0}, 6 baryon, 7 primary (from the collision), 8 other"; fRegistryPairMC.add("PairClass/Pairs/hBremsAncestry", "bremsstrahlung pairs (classes 16 and 5, MC): ancestry of the parent photon (the photon whose electron radiated)", kTHnSparseF, {axSample, {2, -0.5f, 1.5f, "0 class 16 (both true), 1 class 5 (with a fake)"}, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {9, -0.5f, 8.5f, Form("mother of the parent photon: %s", anc)}, {9, -0.5f, 8.5f, Form("mother of that particle: %s", anc)}, {9, -0.5f, 8.5f, Form("root of the decay chain: %s", anc)}}); } fRegistryPairMC.add("PairClass/Pairs/hPhotonKF", "photons of the pairs by MC class (each photon once per pair)", kTHnSparseF, {axSample, axClass, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {2, -0.5f, 1.5f, "0 photon radiated by an e^{#pm} (MC mother e^{#pm}), 1 other"}, {50, 0.f, 25.f, "KF #chi^{2}/ndf of the conversion-vertex fit"}, {60, 0.f, 300.f, "KF #chi^{2} to the primary vertex"}, {40, -0.2f, 0.2f, "#psi_{pair} (rad)"}}); fRegistryPairMC.add("PairClass/Pairs/hHiddenBrems", "true distinct pairs (class 13) with a photon radiated by an e^{#pm} (MC)", kTHnSparseF, {axSample, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {3, -0.5f, 2.5f, "photons of the pair radiated by an e^{#pm}"}, {4, -0.5f, 3.5f, "partner vs the radiating e^{#pm}: 0 unrelated, 1 same mother (e.g. Dalitz #gamma), 2 same root ancestor, 3 partner is itself radiated"}, {6, -0.5f, 5.5f, "radiating e^{#pm} from: 0 photon conversion, 1 #pi^{0} Dalitz, 2 #eta Dalitz, 3 heavy flavour (D, B), 4 other hadron, 5 none / unknown"}, {3, -0.5f, 2.5f, "partner vs the radiating e^{#pm}: 0 different first physical primary, 1 same first physical primary, 2 n/a"}}); + fRegistryPairMC.add("PairClass/Pairs/hPairSep", "photon pairs with q_{inv} < 100 MeV: minimum same-sign leg distance in the TPC by MC class", kTHnSparseF,{axSample, axClass, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {50, 0.f, 10.f, "min #sqrt{#Delta(r#varphi)^{2} + #Delta z^{2}} of the same-sign legs (cm)"}, {25, 0.f, 5.f, "min #Delta(r#varphi) of the same-sign legs (cm)"}, {10, 0.f, 0.5f, "|#it{E}_{1} #minus #it{E}_{2}| (GeV)"}}); fRegistryPairMC.add("PairClass/Photons/hLegTrueDist", "stored photons: distance of the true production vertices of the two legs;0 same event, 1 mixed event;0 SVertexer original, 1 finder original, 2 crossed pairing;distance (cm)", kTH3F, {{2, -0.5f, 1.5f}, {3, -0.5f, 2.5f}, {104, -2.f, 50.f}}); auto hEv = fRegistryPairMC.add("PairClass/Events/hArtificialEvents", "artificial events;;events", kTH1F, {{4, -0.5f, 3.5f}}); hEv->GetXaxis()->SetBinLabel(1, "same event"); @@ -4551,6 +4637,46 @@ struct Photonhbt { } return (a.svertexer ? 0 : 1) + (b.svertexer ? 0 : 1); } + struct PcPairSep { + float dMin = 999.f, dRPhiMin = 999.f; + }; + [[nodiscard]] PcPairSep pcPairSep(PcPhoton const& a, PcPhoton const& b) const + { + PcPairSep s; + if (std::fabs(mBzT) < kMinSigma) { + return s; + } + const std::array la{&a.pos, &a.ele}; + const std::array lb{&b.pos, &b.ele}; + for (int i = 0; i < 2; ++i) { // o2-linter: disable=magic-number (two same-sign leg pairs) + const int charge = i == 0 ? +1 : -1; + const float pt1 = std::hypot(la[i]->px, la[i]->py), pt2 = std::hypot(lb[i]->px, lb[i]->py); + if (pt1 < kMinSigma || pt2 < kMinSigma) { + continue; + } + const float phi1 = std::atan2(la[i]->py, la[i]->px), phi2 = std::atan2(lb[i]->py, lb[i]->px); + const float tgl1 = la[i]->pz / pt1, tgl2 = lb[i]->pz / pt2; + for (size_t ir = 0; ir < kPhiStarRadiiM.size(); ++ir) { + const auto h1 = legHelixAt(a.vx / kCmPerM, a.vy / kCmPerM, phi1, pt1, charge, mBzT, kPhiStarRadiiM[ir]); + const auto h2 = legHelixAt(b.vx / kCmPerM, b.vy / kCmPerM, phi2, pt2, charge, mBzT, kPhiStarRadiiM[ir]); + if (h1.status != 0 || h2.status != 0) { + continue; + } + const float rCm = kPhiStarRadiiM[ir] * kCmPerM; + const float dphi = RecoDecay::constrainAngle(h1.phiStar - h2.phiStar, -o2::constants::math::PI); + const float dRPhi = 2.f * rCm * std::fabs(std::sin(0.5f * dphi)); + const float z1 = a.vz + h1.sxy * kCmPerM * tgl1; + const float z2 = b.vz + h2.sxy * kCmPerM * tgl2; + if (z1 * z2 < 0.f) { + continue; + } + s.dRPhiMin = std::min(s.dRPhiMin, dRPhi); + s.dMin = std::min(s.dMin, std::hypot(dRPhi, std::fabs(z1 - z2))); + } + } + return s; + } + void fillPcPair(const float sample, PcPhoton const& a, PcPhoton const& b) { const float q = o2::pwgem::photonmeson::qInv(a.px, a.py, a.pz, b.px, b.py, b.pz); @@ -4561,6 +4687,11 @@ struct Photonhbt { const auto cls = static_cast(pcPairClass(a, b)); fRegistryPairMC.fill(HIST("PairClass/Pairs/hPairClass"), sample, cls, q, static_cast(pcPairOrigin(a, b)), clean); fRegistryPairMC.fill(HIST("PairClass/Pairs/hKtQ"), sample, cls, q, 0.5f * std::hypot(a.px + b.px, a.py + b.py)); + if (q < 0.1f) { // o2-linter: disable=magic-number (diagnostic range) + const auto sep = pcPairSep(a, b); + const float eAs = std::hypot(a.px, a.py, a.pz), eBs = std::hypot(b.px, b.py, b.pz); + fRegistryPairMC.fill(HIST("PairClass/Pairs/hPairSep"), sample, cls, q, std::min(sep.dMin, 9.99f), std::min(sep.dRPhiMin, 4.99f), std::min(std::fabs(eAs - eBs), 0.499f)); + } if (cls == 4.f) { // o2-linter: disable=magic-number (half swap) for (const auto* x : {&a.pos, &a.ele}) { for (const auto* y : {&b.pos, &b.ele}) { @@ -4572,10 +4703,6 @@ struct Photonhbt { } } } - { - const float eA = std::sqrt(a.px * a.px + a.py * a.py + a.pz * a.pz), eB = std::sqrt(b.px * b.px + b.py * b.py + b.pz * b.pz); - fRegistryPairMC.fill(HIST("PairClass/Pairs/hDEQinv"), sample, cls, q, std::fabs(eA - eB), 0.5f * std::hypot(a.px + b.px, a.py + b.py)); - } fRegistryPairMC.fill(HIST("PairClass/Pairs/hFreeLegPair"), sample, cls, q, static_cast((a.nFree > 0 ? 1 : 0) + (b.nFree > 0 ? 1 : 0)), static_cast((a.original ? 0 : 1) + (b.original ? 0 : 1))); { const ROOT::Math::PxPyPzMVector v1(a.px, a.py, a.pz, 0.f), v2(b.px, b.py, b.pz, 0.f); @@ -4585,6 +4712,8 @@ struct Photonhbt { const float qout = std::fabs(q3.Dot(ROOT::Math::XYZVector(sum.Px(), sum.Py(), 0.).Unit())); const float kt = 0.5f * std::hypot(a.px + b.px, a.py + b.py); fRegistryPairMC.fill(HIST("PairClass/Pairs/hQinvQout"), sample, cls, q, qout, std::min(kt, 2.99f)); + const float eA = std::sqrt(a.px * a.px + a.py * a.py + a.pz * a.pz), eB = std::sqrt(b.px * b.px + b.py * b.py + b.pz * b.pz); + fRegistryPairMC.fill(HIST("PairClass/Pairs/hDEQinv"), sample, cls, q, std::fabs(eA - eB), kt, std::min(qout, 0.299f)); const auto ktUnit = ROOT::Math::XYZVector(sum.Px(), sum.Py(), 0.).Unit(); const float qside = std::fabs(q3.Dot(ROOT::Math::XYZVector(-ktUnit.Y(), ktUnit.X(), 0.))); const float qlong = std::fabs(q3.Z()); @@ -4596,7 +4725,6 @@ struct Photonhbt { fRegistryPairMC.fill(HIST("PairClass/Pairs/hPairOpening"), sample, cls, q, std::min(qside, 0.149f), std::min(qlong, 0.149f), std::min(psi, 1.57f)); } if (q < 0.1f) { // o2-linter: disable=magic-number (diagnostic range) - const float eA = std::hypot(a.px, a.py, a.pz), eB = std::hypot(b.px, b.py, b.pz); PcPhoton const& soft = eA <= eB ? a : b; PcPhoton const& hard = eA <= eB ? b : a; const float rSoft = std::hypot(soft.vx, soft.vy), rHard = std::hypot(hard.vx, hard.vy); @@ -4641,6 +4769,37 @@ struct Photonhbt { const float pmin = std::min({std::hypot(a.pos.px, a.pos.py, a.pos.pz), std::hypot(a.ele.px, a.ele.py, a.ele.pz), std::hypot(b.pos.px, b.pos.py, b.pos.pz), std::hypot(b.ele.px, b.ele.py, b.ele.pz)}); const auto nIts = static_cast((a.pos.its ? 1 : 0) + (a.ele.its ? 1 : 0) + (b.pos.its ? 1 : 0) + (b.ele.its ? 1 : 0)); fRegistryPairMC.fill(HIST("PairClass/Pairs/hPairLegs"), sample, cls, q, std::min(std::max(dtgl(a), dtgl(b)), 0.299f), std::min(pmin, 0.999f), nIts, static_cast(pcPairOrigin(a, b))); + const int icls = static_cast(cls); + if (icls == 1 || icls == 2 || icls == 3 || icls == 12) { // o2-linter: disable=magic-number (duplicate classes) + const std::array la{&a.pos, &a.ele}, lb{&b.pos, &b.ele}; + const PcLeg* x = nullptr; + const PcLeg* y = nullptr; + int slots = 2; + auto shares = [icls](PcLeg const& u, PcLeg const& v) { + return icls == 12 ? (u.evt == v.evt && u.trk >= 0 && u.trk == v.trk) : pcSamePart(u, v); // o2-linter: disable=magic-number (class 12) + }; + for (int i = 0; i < 2 && x == nullptr; ++i) { // o2-linter: disable=magic-number (two legs) + for (int j = 0; j < 2; ++j) { // o2-linter: disable=magic-number (two legs) + if (shares(*la[i], *lb[j])) { + x = la[i]; + y = lb[j]; + slots = i == j ? 0 : 1; + break; + } + } + } + float dTglPieces = 0.199f, dPtPieces = 0.999f, nItsPieces = 0.f; + if (x != nullptr) { + const float pt1 = std::hypot(x->px, x->py), pt2 = std::hypot(y->px, y->py); + if (pt1 > 0.f && pt2 > 0.f) { + dTglPieces = std::min(std::fabs(x->pz / pt1 - y->pz / pt2), 0.199f); + dPtPieces = std::min(std::fabs(pt1 - pt2) / std::max(pt1, pt2), 0.999f); + } + nItsPieces = static_cast((x->its ? 1 : 0) + (y->its ? 1 : 0)); + } + const float clsBin = icls == 12 ? 3.f : static_cast(icls - 1); // o2-linter: disable=magic-number (class 12) + fRegistryPairMC.fill(HIST("PairClass/Pairs/hDupPieces"), sample, clsBin, q, dTglPieces, dPtPieces, static_cast(slots), static_cast(pcPairOrigin(a, b)), nItsPieces); + } constexpr float kTpcRin = 85.f, kTpcRout = 247.f; // o2-linter: disable=magic-number (TPC drift volume, diagnostic) auto crossesCE = [&](PcLeg const& l) { const float pt = std::hypot(l.px, l.py), rConv = std::hypot(l.vx, l.vy); @@ -4764,7 +4923,11 @@ struct Photonhbt { } return 8.f; }; - fRegistryPairMC.fill(HIST("PairClass/Pairs/hBremsAncestry"), sample, cls == 16.f ? 0.f : 1.f, q, ancClass(y->gpdg), ancClass(y->ggpdg), ancClass(y->rootPdg)); // o2-linter: disable=magic-number (class) + // a particle whose mother is the root of the decay chain (the beam record) comes directly from the collision: "primary" (7) + const float ancMother = (y->gmid >= 0 && y->gmid == y->root) ? 7.f : ancClass(y->gpdg); // o2-linter: disable=magic-number (class) + const float ancGrand = (y->ggid >= 0 && y->ggid == y->root) ? 7.f : ancClass(y->ggpdg); // o2-linter: disable=magic-number (class) + + fRegistryPairMC.fill(HIST("PairClass/Pairs/hBremsAncestry"), sample, cls == 16.f ? 0.f : 1.f, q, ancMother, ancGrand, ancClass(y->rootPdg)); // o2-linter: disable=magic-number (class) } } auto alphaOf = [](PcPhoton const& p) { @@ -4857,6 +5020,7 @@ struct Photonhbt { perCollision[photon.collisionId()].push_back(p); } for (const auto& collision : collisions) { + initCCDB(collision); // magnetic field for the leg helices of hPairSep (cached per run) const auto it = perCollision.find(collision.globalIndex()); const std::vector empty; const auto& list = it == perCollision.end() ? empty : it->second; From e199cb25bbc118b5fbc166d5c5b424a4cea97b6f Mon Sep 17 00:00:00 2001 From: Stefanie Mrozinski Date: Thu, 12 Nov 2026 11:31:13 +0100 Subject: [PATCH 2/2] fix formatting --- PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h | 4 ++-- PWGEM/PhotonMeson/Tasks/photonhbt.cxx | 5 ++--- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h index 931ef42e9d4..2b7fea779f2 100644 --- a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h +++ b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h @@ -2887,8 +2887,8 @@ struct PhotonConversionBuilder { } struct BremsPointing { float dMin = 999.f, angle = 3.15f, rEmit = -1.f, outward = 0.f, ptB = 0.f; - int order = -1, leg = -1; // 0 = a radiates b, 1 = b radiates a; leg 0 positive, 1 negative - float ptLeg = 0.f, arc = 0.f, leverD = 0.f; + int order = -1, leg = -1; // 0 = a radiates b, 1 = b radiates a; leg 0 positive, 1 negative + float ptLeg = 0.f, arc = 0.f, leverD = 0.f; float score = 1.e9f, missT = 999.f, thetaB = 3.15f, ptBV = 0.f, lever = 0.f; int orderV = -1, legV = -1; }; diff --git a/PWGEM/PhotonMeson/Tasks/photonhbt.cxx b/PWGEM/PhotonMeson/Tasks/photonhbt.cxx index 5903fc5904b..effde0b149e 100644 --- a/PWGEM/PhotonMeson/Tasks/photonhbt.cxx +++ b/PWGEM/PhotonMeson/Tasks/photonhbt.cxx @@ -1327,7 +1327,7 @@ struct Photonhbt { fRegistryPairMC.add((spBase + "hEtaPt").c_str(), "#eta vs p_{T};#eta;p_{T} (GeV/c)", kTH2F, {axisEta, axisPt}, true); fRegistryPairMC.add((spBase + "hLegOriginMatrix").c_str(), "pos-leg origin vs neg-leg origin;pos leg;neg leg", kTH2F, {axisLegOrigin, axisLegOrigin}, true); } - fRegistryPairMC.add("SinglePhotonMC/Response/hResponse", "true reconstructed photons: response (MC)", kTHnSparseF,{{40, 0.f, 4.f, "#it{p}_{T}^{true} (GeV/#it{c})"},{60, 0.3f, 1.5f, "#it{E}^{rec} / #it{E}^{true}"}, {40, 0.f, 0.04f, "angle(rec, true) (rad)"},{18, 0.f, 90.f, "#it{R}_{conv} (cm)"},{21, -0.05f, 1.f, "hardest photon radiated by a leg: #it{E} / #it{E}^{true} (first bin: none)"},{2, -0.5f, 1.5f, "that radiated photon reconstructed in this event"},{3, -0.5f, 2.5f, "this photon radiated by a conversion e^{#pm}: 0 no, 1 yes and parent #gamma lost, 2 yes and parent reconstructed"},{21, -0.05f, 1.f, "#it{E}^{true} / #it{E}(parent #gamma) (first bin: n/a)"},{51, -0.005f, 0.25f, "angle(this, parent #gamma) MC (rad) (first bin: n/a)"}},true); + fRegistryPairMC.add("SinglePhotonMC/Response/hResponse", "true reconstructed photons: response (MC)", kTHnSparseF, {{40, 0.f, 4.f, "#it{p}_{T}^{true} (GeV/#it{c})"}, {60, 0.3f, 1.5f, "#it{E}^{rec} / #it{E}^{true}"}, {40, 0.f, 0.04f, "angle(rec, true) (rad)"}, {18, 0.f, 90.f, "#it{R}_{conv} (cm)"}, {21, -0.05f, 1.f, "hardest photon radiated by a leg: #it{E} / #it{E}^{true} (first bin: none)"}, {2, -0.5f, 1.5f, "that radiated photon reconstructed in this event"}, {3, -0.5f, 2.5f, "this photon radiated by a conversion e^{#pm}: 0 no, 1 yes and parent #gamma lost, 2 yes and parent reconstructed"}, {21, -0.05f, 1.f, "#it{E}^{true} / #it{E}(parent #gamma) (first bin: n/a)"}, {51, -0.005f, 0.25f, "angle(this, parent #gamma) MC (rad) (first bin: n/a)"}}, true); fRegistryPairMC.add("SinglePhotonMC/Response/hGenPt", "generated primary photons of the selected collisions, |#eta| < cfgMaxEtaV0;#it{p}_{T} (GeV/#it{c})", kTH1F, {{80, 0.f, 4.f}}, true); const AxisSpec axisFakeSubtype{3, -0.5f, 2.5f, "0=ordinary, 1=swap doublet, 2=partial swap"}; @@ -2987,7 +2987,6 @@ struct Photonhbt { } } - template void fillSinglePhotonResponse(TCollision const& collision, TPhotons const& photonsColl, TMCParticles const& mcParticles, TCut const& cut) { @@ -4490,7 +4489,7 @@ struct Photonhbt { } fRegistryPairMC.add("PairClass/Pairs/hPhotonKF", "photons of the pairs by MC class (each photon once per pair)", kTHnSparseF, {axSample, axClass, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {2, -0.5f, 1.5f, "0 photon radiated by an e^{#pm} (MC mother e^{#pm}), 1 other"}, {50, 0.f, 25.f, "KF #chi^{2}/ndf of the conversion-vertex fit"}, {60, 0.f, 300.f, "KF #chi^{2} to the primary vertex"}, {40, -0.2f, 0.2f, "#psi_{pair} (rad)"}}); fRegistryPairMC.add("PairClass/Pairs/hHiddenBrems", "true distinct pairs (class 13) with a photon radiated by an e^{#pm} (MC)", kTHnSparseF, {axSample, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {3, -0.5f, 2.5f, "photons of the pair radiated by an e^{#pm}"}, {4, -0.5f, 3.5f, "partner vs the radiating e^{#pm}: 0 unrelated, 1 same mother (e.g. Dalitz #gamma), 2 same root ancestor, 3 partner is itself radiated"}, {6, -0.5f, 5.5f, "radiating e^{#pm} from: 0 photon conversion, 1 #pi^{0} Dalitz, 2 #eta Dalitz, 3 heavy flavour (D, B), 4 other hadron, 5 none / unknown"}, {3, -0.5f, 2.5f, "partner vs the radiating e^{#pm}: 0 different first physical primary, 1 same first physical primary, 2 n/a"}}); - fRegistryPairMC.add("PairClass/Pairs/hPairSep", "photon pairs with q_{inv} < 100 MeV: minimum same-sign leg distance in the TPC by MC class", kTHnSparseF,{axSample, axClass, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {50, 0.f, 10.f, "min #sqrt{#Delta(r#varphi)^{2} + #Delta z^{2}} of the same-sign legs (cm)"}, {25, 0.f, 5.f, "min #Delta(r#varphi) of the same-sign legs (cm)"}, {10, 0.f, 0.5f, "|#it{E}_{1} #minus #it{E}_{2}| (GeV)"}}); + fRegistryPairMC.add("PairClass/Pairs/hPairSep", "photon pairs with q_{inv} < 100 MeV: minimum same-sign leg distance in the TPC by MC class", kTHnSparseF, {axSample, axClass, {20, 0.f, 0.1f, "#it{q}_{inv} (GeV/#it{c})"}, {50, 0.f, 10.f, "min #sqrt{#Delta(r#varphi)^{2} + #Delta z^{2}} of the same-sign legs (cm)"}, {25, 0.f, 5.f, "min #Delta(r#varphi) of the same-sign legs (cm)"}, {10, 0.f, 0.5f, "|#it{E}_{1} #minus #it{E}_{2}| (GeV)"}}); fRegistryPairMC.add("PairClass/Photons/hLegTrueDist", "stored photons: distance of the true production vertices of the two legs;0 same event, 1 mixed event;0 SVertexer original, 1 finder original, 2 crossed pairing;distance (cm)", kTH3F, {{2, -0.5f, 1.5f}, {3, -0.5f, 2.5f}, {104, -2.f, 50.f}}); auto hEv = fRegistryPairMC.add("PairClass/Events/hArtificialEvents", "artificial events;;events", kTH1F, {{4, -0.5f, 3.5f}}); hEv->GetXaxis()->SetBinLabel(1, "same event");