diff --git a/PWGEM/PhotonMeson/Tasks/photonhbt.cxx b/PWGEM/PhotonMeson/Tasks/photonhbt.cxx index bacb008cd47..40a6599bb71 100644 --- a/PWGEM/PhotonMeson/Tasks/photonhbt.cxx +++ b/PWGEM/PhotonMeson/Tasks/photonhbt.cxx @@ -10,11 +10,13 @@ // or submit itself to any jurisdiction. /// \file photonhbt.cxx -/// \brief This code loops over v0 photons and makes pairs for photon HBT analysis. -/// \author Daiki Sekihata, daiki.sekihata@cern.ch and Stefanie Mrozinski stefanie.mrozinski@cern.ch +/// \brief V0 photon HBT analysis. +/// \author Daiki Sekihata, daiki.sekihata@cern.ch +/// Stefanie Mrozinski, stefanie.mrozinski@cern.ch #include "PWGEM/Dilepton/Utils/EMTrack.h" #include "PWGEM/Dilepton/Utils/EventMixingHandler.h" +#include "PWGEM/Dilepton/Utils/MCUtilities.h" #include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" #include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" #include "PWGEM/PhotonMeson/DataModel/EventTables.h" @@ -53,6 +55,7 @@ #include #include #include +#include #include #include @@ -60,8 +63,8 @@ enum class V0Combo : int { Inclusive = 0, ItstpcItstpc = 1, ///< both legs ITS+TPC - ItstpcTpconly = 2, /// one ITS+TPC leg, one TPC-only - TpconlyTpconly = 3, /// both legs TPC-only + ItstpcTpconly = 2, ///< one ITS+TPC leg, one TPC-only + TpconlyTpconly = 3, ///< both legs TPC-only }; /// Photon-pair track-type combo. @@ -82,12 +85,43 @@ using namespace o2::framework::expressions; using namespace o2::soa; using namespace o2::aod::pwgem::dilepton::utils; -using MyCollisions = soa::Join; +using MyCollisions = soa::Join; using MyCollision = MyCollisions::iterator; using MyV0Photons = soa::Join; using MyV0Photon = MyV0Photons::iterator; +using MyMCV0Legs = soa::Join; +using MyMCV0Leg = MyMCV0Legs::iterator; + +// ─── MC truth classification types ──────────────────────────────────────────── + +/// Per-photon MC truth information built from the two V0 legs. +struct PhotonMCInfo { + bool hasMC = false; // both legs have a valid MC label + bool sameMother = false; // both legs share the same MC mother + bool isTruePhoton = false; // mother PDG == 22 + + int mcPosId = -1; // MC particle index of the positive leg + int mcNegId = -1; // MC particle index of the negative leg + int motherId = -1; // MC particle index of the common mother + int motherPdg = 0; + + bool isPhysicalPrimary = false; +}; + +/// Classification of a photon pair at the MC-truth level. +enum class PairTruthType : uint8_t { + Unknown = 0, + TrueTrueDistinct, // both photons are true, from different MC photons + TrueTrueSamePhoton, // both photons are true, same MC photon (clone/split) + SharedMcLeg, // different reco tracks but same MC-level leg + TrueFake, // one photon is true, one is fake + FakeFake, // both photons are fake + Pi0Daughters, // both photons come from the same MC pi0 +}; + struct photonhbt { template @@ -116,13 +150,12 @@ struct photonhbt { return PairCombo::Inclusive; const int lo = std::min(i1, i2); const int hi = std::max(i1, i2); - static constexpr std::array, 4> kTable = {{{0, 0, 0, 0}, {0, 1, 2, 3}, {0, 2, 4, 5}, {0, 3, 5, 6}}}; + static constexpr std::array, 4> kTable = { + {{0, 0, 0, 0}, {0, 1, 2, 3}, {0, 2, 4, 5}, {0, 3, 5, 6}}}; return static_cast(kTable[lo][hi]); } - // --------------------------------------------------------------------------- - // Configurables: histogram axes - // --------------------------------------------------------------------------- + // ─── Configurables: histogram axis bins ─────────────────────────────────── // HBT physics ConfigurableAxis confQBins{"confQBins", {60, 0, +0.3f}, "q bins for output histograms"}; @@ -131,16 +164,29 @@ struct photonhbt { // Single-photon QA ConfigurableAxis confPtBins{"confPtBins", {100, 0.f, 2.f}, "pT bins (GeV/c)"}; ConfigurableAxis confEtaBins{"confEtaBins", {80, -0.8f, 0.8f}, "eta bins"}; - ConfigurableAxis confPhiBins{"confPhiBins", {90, -o2::constants::math::PI, o2::constants::math::PI}, "phi bins (rad)"}; + ConfigurableAxis confPhiBins{"confPhiBins", {90, 0.f, o2::constants::math::TwoPI}, "phi bins (rad) — O2 track phi is in [0, 2pi]"}; - // Pair QA + // Pair angular ConfigurableAxis confDeltaEtaBins{"confDeltaEtaBins", {100, -0.9f, +0.9f}, "Delta-eta bins"}; ConfigurableAxis confDeltaPhiBins{"confDeltaPhiBins", {100, -o2::constants::math::PI, o2::constants::math::PI}, "Delta-phi bins (rad)"}; ConfigurableAxis confEllipseValBins{"confEllipseValBins", {200, 0.f, 10.f}, "ellipse value bins"}; ConfigurableAxis confCosThetaBins{"confCosThetaBins", {100, 0.f, 1.f}, "cos(theta*) bins"}; ConfigurableAxis confOpeningAngleBins{"confOpeningAngleBins", {100, 0.f, o2::constants::math::PI}, "opening angle bins (rad)"}; - // Axis specs + // Pair geometry + ConfigurableAxis confRBins{"confRBins", {100, 0.f, 100.f}, "conversion radius bins (cm)"}; + ConfigurableAxis confDeltaRBins{"confDeltaRBins", {120, 0.f, 30.f}, "|R1-R2| bins (cm)"}; + ConfigurableAxis confDeltaR3DBins{"confDeltaR3DBins", {100, 0.f, 100.f}, "3D distance between conversion points (cm)"}; + ConfigurableAxis confDeltaRxyBins{"confDeltaRxyBins", {100, 0.f, 100.f}, "xy distance between conversion points (cm)"}; + ConfigurableAxis confZConvBins{"confZConvBins", {200, -100.f, 100.f}, "conversion z (cm)"}; + ConfigurableAxis confDeltaZBins{"confDeltaZBins", {200, -100.f, 100.f}, "#Deltaz bins (cm)"}; ///< FIX: was missing + + // Event QA + ConfigurableAxis confOccupancyQA{"confOccupancyQA", {100, 0.f, 50000.f}, "occupancy"}; + ConfigurableAxis confCentQABins{"confCentQABins", {110, 0.f, 110.f}, "centrality (%)"}; + + // ─── Axis specs ──────────────────────────────────────────────────────────── + const AxisSpec axisKt{confKtBins, "k_{T} (GeV/c)"}; const AxisSpec axisQinv{confQBins, "q_{inv} (GeV/c)"}; const AxisSpec axisQabsLcms{confQBins, "|#bf{q}|^{LCMS} (GeV/c)"}; @@ -155,39 +201,51 @@ struct photonhbt { const AxisSpec axisEllipseVal{confEllipseValBins, "(#Delta#eta/#sigma_{#eta})^{2}+(#Delta#phi/#sigma_{#phi})^{2}"}; const AxisSpec axisCosTheta{confCosThetaBins, "cos(#theta*)"}; const AxisSpec axisOpeningAngle{confOpeningAngleBins, "Opening angle (rad)"}; + const AxisSpec axisR{confRBins, "R_{conv} (cm)"}; + const AxisSpec axisDeltaR{confDeltaRBins, "|R_{1}-R_{2}| (cm)"}; + const AxisSpec axisDeltaR3D{confDeltaR3DBins, "|#vec{r}_{1}-#vec{r}_{2}| (cm)"}; + const AxisSpec axisDeltaRxy{confDeltaRxyBins, "#Delta r_{xy} (cm)"}; + const AxisSpec axisZConv{confZConvBins, "z_{conv} (cm)"}; + const AxisSpec axisDeltaZ{confDeltaZBins, "#Delta z (cm)"}; ///< FIX: was missing + const AxisSpec axisOccupancy{confOccupancyQA, "occupancy"}; + const AxisSpec axisCentQA{confCentQABins, "centrality (%)"}; - // --------------------------------------------------------------------------- - // Configurables: QA flags - // --------------------------------------------------------------------------- + // ─── Configurables: QA flags ─────────────────────────────────────────────── struct : ConfigurableGroup { std::string prefix = "qaflags_group"; - Configurable doPairQa{"doPairQa", true, "fill pair QA histograms (Before/After ellipse cut)"}; + Configurable doPairQa{"doPairQa", true, "fill pair QA histograms at each cut step"}; Configurable doSinglePhotonQa{"doSinglePhotonQa", true, "fill single-photon QA histograms (pT, eta, phi)"}; + Configurable cfgMaxQinvForQA{"cfgMaxQinvForQA", 0.1f, + "fill per-step pair QA histograms (hDeltaEta, hDeltaPhi, THnSparses, ...) " + "only when q_inv < this value (GeV/c). " + "Set to the HBT signal region, typically 0.1. " + "Set <= 0 to disable the gate. The CF is always filled regardless."}; + Configurable cfgMaxQinvForFullRange{"cfgMaxQinvForFullRange", 0.3f, + "fill full-range histograms (hDeltaRVsQinv, hSparseDeltaRDeltaZQinv, ...) " + "only when q_inv < this value (GeV/c). " + "Should match the upper edge of confQBins (default 0.3) — " + "fills beyond the axis range only go into the overflow bin. " + "Set <= 0 to disable the gate."}; } qaflags; - // --------------------------------------------------------------------------- - // Configurables: HBT kind - // --------------------------------------------------------------------------- + // ─── Configurables: HBT kind ─────────────────────────────────────────────── - Configurable cfgDo3D{"cfgDo3D", false, "enable 3D analysis"}; - Configurable cfgUseLCMS{"cfgUseLCMS", true, "measure relative momentum in LCMS for 1D"}; + Configurable cfgDo3D{"cfgDo3D", false, "enable 3D (qout,qside,qlong) analysis"}; + Configurable cfgDo2D{"cfgDo2D", false, "enable 2D (qout,qinv) projection (requires cfgDo3D)"}; + Configurable cfgUseLCMS{"cfgUseLCMS", true, "measure 1D relative momentum in LCMS"}; - // --------------------------------------------------------------------------- - // Configurables: events - // --------------------------------------------------------------------------- + // ─── Configurables: events ───────────────────────────────────────────────── Configurable cfgCentMin{"cfgCentMin", -1, "min. centrality"}; Configurable cfgCentMax{"cfgCentMax", 999, "max. centrality"}; Configurable maxY{"maxY", 0.9, "maximum rapidity"}; - // --------------------------------------------------------------------------- - // Configurables: mixed event - // --------------------------------------------------------------------------- + // ─── Configurables: mixed event ──────────────────────────────────────────── Configurable cfgDoMix{"cfgDoMix", true, "flag for event mixing"}; Configurable ndepth{"ndepth", 100, "depth for event mixing"}; - Configurable ndiffBCMix{"ndiffBCMix", 594, "difference in global BC required in mixed events"}; + Configurable ndiffBCMix{"ndiffBCMix", 594, "difference in global BC required for mixed events"}; Configurable cfgEP2EstimatorForMix{"cfgEP2EstimatorForMix", 3, "FT0M:0, FT0A:1, FT0C:2, FV0A:3, BTot:4, BPos:5, BNeg:6"}; Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; Configurable cfgOccupancyEstimator{"cfgOccupancyEstimator", 0, "FT0C:0, Track:1"}; @@ -197,22 +255,25 @@ struct photonhbt { ConfigurableAxis ConfEPBins{"ConfEPBins", {16, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}, "Mixing bins - EP angle"}; ConfigurableAxis ConfOccupancyBins{"ConfOccupancyBins", {VARIABLE_WIDTH, -1, 1e+10}, "Mixing bins - occupancy"}; - // --------------------------------------------------------------------------- - // Configurables: pair cuts - // --------------------------------------------------------------------------- + // ─── Configurables: pair cuts ────────────────────────────────────────────── struct : ConfigurableGroup { std::string prefix = "ggpaircut_group"; - Configurable cfgMinDR_CosOA{"cfgMinDR_CosOA", -1, "min. dr/cosOA for kPCMPCM"}; - Configurable cfgApplyEllipseCut{"cfgApplyEllipseCut", false, "reject pairs inside ellipse in DeltaEta-DeltaPhi"}; + // dr/cosOA cut + Configurable cfgMinDRCosOA{"cfgMinDRCosOA", -1.f, "min. dr/cosOA; <0 = disabled"}; + // R/Z geometry cuts + Configurable cfgDoRCut{"cfgDoRCut", false, "apply |R1-R2| > cfgMinDeltaR cut"}; + Configurable cfgMinDeltaR{"cfgMinDeltaR", 0.f, "minimum |R1-R2| (cm)"}; + Configurable cfgDoZCut{"cfgDoZCut", false, "apply |DeltaZ| > cfgMinDeltaZ cut"}; + Configurable cfgMinDeltaZ{"cfgMinDeltaZ", 0.f, "minimum |DeltaZ| (cm)"}; + // Ellipse cut in (DeltaEta, DeltaPhi) + Configurable cfgDoEllipseCut{"cfgDoEllipseCut", false, "reject pairs inside ellipse in DeltaEta-DeltaPhi"}; Configurable cfgEllipseSigEta{"cfgEllipseSigEta", 0.02f, "sigma_eta for ellipse cut"}; Configurable cfgEllipseSigPhi{"cfgEllipseSigPhi", 0.02f, "sigma_phi for ellipse cut"}; - Configurable cfgEllipseR2{"cfgEllipseR2", 1.0f, "R^2 threshold: reject if value < R^2"}; + Configurable cfgEllipseR2{"cfgEllipseR2", 1.0f, "R^2 threshold: reject if ellipse value < R^2"}; } ggpaircuts; - // --------------------------------------------------------------------------- - // Event cut - // --------------------------------------------------------------------------- + // ─── Event cut ───────────────────────────────────────────────────────────── EMPhotonEventCut fEMEventCut; struct : ConfigurableGroup { @@ -240,40 +301,31 @@ struct photonhbt { Configurable cfgRequireGoodITSLayersAll{"cfgRequireGoodITSLayersAll", false, "all ITS layers OK"}; } eventcuts; - // --------------------------------------------------------------------------- - // PCM cut - // --------------------------------------------------------------------------- + // ─── PCM cut ─────────────────────────────────────────────────────────────── V0PhotonCut fV0PhotonCut; struct : ConfigurableGroup { std::string prefix = "pcmcut_group"; - Configurable cfgRequireV0WithITSTPC{"cfgRequireV0WithITSTPC", false, "select V0s with ITS-TPC tracks"}; Configurable cfgRequireV0WithITSOnly{"cfgRequireV0WithITSOnly", false, "select V0s with ITS-only tracks"}; Configurable cfgRequireV0WithTPCOnly{"cfgRequireV0WithTPCOnly", false, "select V0s with TPC-only tracks"}; - Configurable cfgMinPtV0{"cfgMinPtV0", 0.1, "min pT for V0 photons at PV"}; Configurable cfgMaxEtaV0{"cfgMaxEtaV0", 0.8, "max eta for V0 photons at PV"}; Configurable cfgMinV0Radius{"cfgMinV0Radius", 16.0, "min V0 radius"}; Configurable cfgMaxV0Radius{"cfgMaxV0Radius", 90.0, "max V0 radius"}; - Configurable cfgMaxAlphaAP{"cfgMaxAlphaAP", 0.95, "max alpha for AP cut"}; Configurable cfgMaxQtAP{"cfgMaxQtAP", 0.01, "max qT for AP cut"}; - Configurable cfgMinCosPA{"cfgMinCosPA", 0.997, "min V0 CosPA"}; Configurable cfgMaxPCA{"cfgMaxPCA", 3.0, "max distance between 2 legs"}; Configurable cfgMaxChi2KF{"cfgMaxChi2KF", 1e+10, "max chi2/ndf with KF"}; - Configurable cfgRejectV0OnITSIB{"cfgRejectV0OnITSIB", true, "reject V0s on ITSib"}; Configurable cfgDisableITSOnlyTrack{"cfgDisableITSOnlyTrack", false, "disable ITS-only tracks"}; Configurable cfgDisableTPCOnlyTrack{"cfgDisableTPCOnlyTrack", false, "disable TPC-only tracks"}; - Configurable cfgMinNClusterTPC{"cfgMinNClusterTPC", 70, "min ncluster TPC"}; Configurable cfgMinNCrossedRows{"cfgMinNCrossedRows", 70, "min crossed rows"}; Configurable cfgMaxFracSharedClustersTPC{"cfgMaxFracSharedClustersTPC", 999.f, "max fraction of shared TPC clusters"}; Configurable cfgMaxChi2TPC{"cfgMaxChi2TPC", 4.0, "max chi2/NclsTPC"}; Configurable cfgMaxChi2ITS{"cfgMaxChi2ITS", 36.0, "max chi2/NclsITS"}; - Configurable cfgMinTPCNsigmaEl{"cfgMinTPCNsigmaEl", -3.5, "min TPC nsigma electron"}; Configurable cfgMaxTPCNsigmaEl{"cfgMaxTPCNsigmaEl", +3.5, "max TPC nsigma electron"}; } pcmcuts; @@ -286,7 +338,6 @@ struct photonhbt { emh2 = nullptr; mapMixedEventIdToGlobalBC.clear(); usedPhotonIdsPerCol.clear(); - usedPhotonIdsPerCol.shrink_to_fit(); } HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; @@ -294,16 +345,18 @@ struct photonhbt { std::mt19937 engine; std::uniform_int_distribution dist01; - int mRunNumber; + int mRunNumber{0}; std::vector ztxBinEdges; std::vector centBinEdges; std::vector epBinEgdes; std::vector occBinEdges; + // ─── Pair-cut helpers ────────────────────────────────────────────────────── + inline bool isInsideEllipse(float deta, float dphi) const { - if (!ggpaircuts.cfgApplyEllipseCut.value) + if (!ggpaircuts.cfgDoEllipseCut.value) // .value needed: operator T() is non-const in O2 return false; const float sE = ggpaircuts.cfgEllipseSigEta.value; const float sP = ggpaircuts.cfgEllipseSigPhi.value; @@ -312,6 +365,27 @@ struct photonhbt { return (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP) < ggpaircuts.cfgEllipseR2.value; } + inline bool passRZCut(float deltaR, float deltaZ) const + { + if (ggpaircuts.cfgDoRCut.value && deltaR < ggpaircuts.cfgMinDeltaR.value) + return false; + if (ggpaircuts.cfgDoZCut.value && std::fabs(deltaZ) < ggpaircuts.cfgMinDeltaZ.value) + return false; + return true; + } + + inline bool passQinvQAGate(float qinv) const + { + const float limit = qaflags.cfgMaxQinvForQA.value; + return (limit <= 0.f) || (qinv < limit); + } + + inline bool passQinvFullRangeGate(float qinv) const + { + const float limit = qaflags.cfgMaxQinvForFullRange.value; + return (limit <= 0.f) || (qinv < limit); + } + static inline float computeCosTheta(const ROOT::Math::PtEtaPhiMVector& v1, const ROOT::Math::PtEtaPhiMVector& v2) { @@ -341,16 +415,48 @@ struct photonhbt { } } - /// Clamp bin index to valid range [0, nmax]. static int clampBin(int b, int nmax) { return std::clamp(b, 0, nmax); } - /// Find the bin index for val in a sorted edge vector. static int binOf(const std::vector& edges, float val) { - const int b = static_cast(std::lower_bound(edges.begin(), edges.end(), val) - edges.begin()) - 1; + const int b = static_cast( + std::lower_bound(edges.begin(), edges.end(), val) - edges.begin()) - + 1; return clampBin(b, static_cast(edges.size()) - 2); } + /// ev_id : 0 = same-event, 1 = mixed-event + /// step_id: 0 = Before, 1 = AfterDRCosOA, 2 = AfterRZ, 3 = AfterEllipse + template + static constexpr const char* qaPrefix() + { + if constexpr (ev_id == 0) { + if constexpr (step_id == 0) + return "Pair/same/QA/Before/"; + if constexpr (step_id == 1) + return "Pair/same/QA/AfterDRCosOA/"; + if constexpr (step_id == 2) + return "Pair/same/QA/AfterRZ/"; + return "Pair/same/QA/AfterEllipse/"; + } else { + if constexpr (step_id == 0) + return "Pair/mix/QA/Before/"; + if constexpr (step_id == 1) + return "Pair/mix/QA/AfterDRCosOA/"; + if constexpr (step_id == 2) + return "Pair/mix/QA/AfterRZ/"; + return "Pair/mix/QA/AfterEllipse/"; + } + } + + template + static constexpr const char* fullRangePrefix() + { + if constexpr (ev_id == 0) + return "Pair/same/FullRange/"; + return "Pair/mix/FullRange/"; + } + void init(InitContext& /*context*/) { mRunNumber = 0; @@ -373,7 +479,7 @@ struct photonhbt { dist01 = std::uniform_int_distribution(0, 1); fRegistry.add("Pair/mix/hDiffBC", - "diff. global BC in mixed event;|BC_{current} - BC_{mixed}|", + "diff. global BC in mixed event;|BC_{current}-BC_{mixed}|", kTH1D, {{10001, -0.5, 10000.5}}, true); } @@ -384,51 +490,198 @@ struct photonhbt { return; mRunNumber = collision.runNumber(); } + + // ─── PairQAObservables ───────────────────────────────────────────────────── + /// Plain data struct holding all observables computed from a photon pair. + /// Defined early so all histogram-booking and fill functions can use it. + + struct PairQAObservables { + // photon four-vectors and pair kinematics + ROOT::Math::PtEtaPhiMVector v1; + ROOT::Math::PtEtaPhiMVector v2; + ROOT::Math::PtEtaPhiMVector k12; + // conversion-point coordinates + float x1 = 0.f, y1 = 0.f, z1 = 0.f; + float x2 = 0.f, y2 = 0.f, z2 = 0.f; + // conversion-point radii and distances + float r1 = 0.f, r2 = 0.f; + float dx = 0.f, dy = 0.f, dz = 0.f; + float deltaR = 0.f; ///< |R1-R2| + float deltaZ = 0.f; ///< z1-z2 + float deltaRxy = 0.f; ///< sqrt(dx^2+dy^2) + float deltaR3D = 0.f; ///< sqrt(dx^2+dy^2+dz^2) + // opening angle of conversion-point vectors + float opa = 0.f; + float cosOA = 0.f; + float drOverCosOA = 0.f; + float deta = 0.f, dphi = 0.f; + float pairEta = 0.f, pairPhi = 0.f; + float kt = 0.f, qinv = 0.f; + float cosTheta = 0.f; + float openingAngle = 0.f; + // validity flag + bool valid = true; + }; + + void addSinglePhotonQAHistogramsForStep(const std::string& path) + { + fRegistry.add((path + "hPt").c_str(), "p_{T};p_{T} (GeV/c);counts", kTH1D, {axisPt}, true); + fRegistry.add((path + "hEta").c_str(), "#eta;#eta;counts", kTH1D, {axisEta}, true); + fRegistry.add((path + "hPhi").c_str(), "#phi;#phi (rad);counts", kTH1D, {axisPhi}, true); + fRegistry.add((path + "hEtaVsPhi").c_str(), "acceptance;#phi (rad);#eta", kTH2D, {axisPhi, axisEta}, true); + fRegistry.add((path + "hR").c_str(), "R_{conv};R_{conv} (cm);counts", kTH1D, {axisR}, true); + fRegistry.add((path + "hZConv").c_str(), "z_{conv};z_{conv} (cm);counts", kTH1D, {axisZConv}, true); + fRegistry.add((path + "hRVsZConv").c_str(), "R_{conv} vs z_{conv};z_{conv} (cm);R_{conv} (cm)", kTH2D, {axisZConv, axisR}, true); + } + + void addFullRangeHistograms(const std::string& path) + { + fRegistry.add((path + "hDeltaRVsQinv").c_str(), "|R_{1}-R_{2}| vs q_{inv} (full range);q_{inv} (GeV/c);|R_{1}-R_{2}| (cm)", kTH2D, {axisQinv, axisDeltaR}, true); + fRegistry.add((path + "hDeltaZVsQinv").c_str(), "#Delta z vs q_{inv} (full range);q_{inv} (GeV/c);#Delta z (cm)", kTH2D, {axisQinv, axisDeltaZ}, true); + fRegistry.add((path + "hDeltaR3DVsQinv").c_str(), "#Delta r_{3D} vs q_{inv} (full range);q_{inv} (GeV/c);#Delta r_{3D} (cm)", kTH2D, {axisQinv, axisDeltaR3D}, true); + fRegistry.add((path + "hQinvVsCent").c_str(), "q_{inv} vs centrality (full range);centrality (%);q_{inv} (GeV/c)", kTH2D, {axisCentQA, axisQinv}, true); + fRegistry.add((path + "hQinvVsOccupancy").c_str(), "q_{inv} vs occupancy (full range);occupancy;q_{inv} (GeV/c)", kTH2D, {axisOccupancy, axisQinv}, true); + fRegistry.add((path + "hSparseDeltaRDeltaZQinv").c_str(), "|R_{1}-R_{2}|,#Delta z,q_{inv} (full range)", kTHnSparseD, {axisDeltaR, axisDeltaZ, axisQinv}, true); + fRegistry.add((path + "hDeltaRCosOAVsQinv").c_str(), "#Delta r/cos(#theta_{op}/2) vs q_{inv};q_{inv} (GeV/c);#Delta r/cos(#theta_{op}/2) (cm)", kTH2D, {axisQinv, {100, 0, 100}}, true); + } + + /// ev_id : 0 = same-event, 1 = mixed-event + template + inline void fillFullRangeQA(PairQAObservables const& obs, float cent, float occupancy) + { + constexpr auto base = fullRangePrefix(); + fRegistry.fill(HIST(base) + HIST("hDeltaRVsQinv"), obs.qinv, obs.deltaR); + fRegistry.fill(HIST(base) + HIST("hDeltaZVsQinv"), obs.qinv, obs.deltaZ); + fRegistry.fill(HIST(base) + HIST("hDeltaR3DVsQinv"), obs.qinv, obs.deltaR3D); + fRegistry.fill(HIST(base) + HIST("hQinvVsCent"), cent, obs.qinv); + fRegistry.fill(HIST(base) + HIST("hQinvVsOccupancy"), occupancy, obs.qinv); + fRegistry.fill(HIST(base) + HIST("hSparseDeltaRDeltaZQinv"), + obs.deltaR, obs.deltaZ, obs.qinv); + } + + template + inline void fillFullRangeDeltaRCosOA(float qinv, float drOverCosOA) + { + constexpr auto base = fullRangePrefix(); + fRegistry.fill(HIST(base) + HIST("hDeltaRCosOAVsQinv"), qinv, drOverCosOA); + } + void addQAHistogramsForStep(const std::string& path) { - fRegistry.add((path + "hDeltaEta").c_str(), "#Delta#eta;#Delta#eta;counts", kTH1F, {axisDeltaEta}, true); - fRegistry.add((path + "hDeltaPhi").c_str(), "#Delta#phi;#Delta#phi (rad);counts", kTH1F, {axisDeltaPhi}, true); - fRegistry.add((path + "hDEtaDPhi").c_str(), "#Delta#eta vs #Delta#phi;#Delta#eta;#Delta#phi (rad)", kTH2F, {axisDeltaEta, axisDeltaPhi}, true); - fRegistry.add((path + "hDeltaEtaVsPairEta").c_str(), "#Delta#eta vs #LT#eta#GT_{pair};#LT#eta#GT_{pair};#Delta#eta", kTH2F, {axisEta, axisDeltaEta}, true); - fRegistry.add((path + "hCosTheta").c_str(), "cos(#theta*) in pair rest frame;cos(#theta*);counts", kTH1F, {axisCosTheta}, true); - fRegistry.add((path + "hOpeningAngle").c_str(), "Opening angle between conversion points;#alpha (rad);counts", kTH1F, {axisOpeningAngle}, true); + // ── 1D: photon kinematics ──────────────────────────────────────────────── + fRegistry.add((path + "hPairEta").c_str(), "pair #eta;#eta_{pair};counts", kTH1D, {axisEta}, true); + fRegistry.add((path + "hPairPhi").c_str(), "pair #phi;#phi_{pair} (rad);counts", kTH1D, {axisPhi}, true); + fRegistry.add((path + "hPairKt").c_str(), "pair k_{T};k_{T} (GeV/c);counts", kTH1D, {axisKt}, true); + fRegistry.add((path + "hQinv").c_str(), "q_{inv};q_{inv} (GeV/c);counts", kTH1D, {axisQinv}, true); + + // ── 1D: angular ───────────────────────────────────────────────────────── + fRegistry.add((path + "hDeltaEta").c_str(), "#Delta#eta;#Delta#eta;counts", kTH1D, {axisDeltaEta}, true); + fRegistry.add((path + "hDeltaPhi").c_str(), "#Delta#phi;#Delta#phi (rad);counts", kTH1D, {axisDeltaPhi}, true); + fRegistry.add((path + "hCosTheta").c_str(), "cos(#theta*) in pair rest frame;cos(#theta*);counts", kTH1D, {axisCosTheta}, true); + fRegistry.add((path + "hOpeningAngle").c_str(), "Opening angle;#alpha (rad);counts", kTH1D, {axisOpeningAngle}, true); fRegistry.add((path + "hEllipseVal").c_str(), "(#Delta#eta/#sigma_{#eta})^{2}+(#Delta#phi/#sigma_{#phi})^{2};value;counts", kTH1D, {axisEllipseVal}, true); + + // ── 1D: geometry ──────────────────────────────────────────────────────── + fRegistry.add((path + "hR1").c_str(), "R_{conv,1};R_{1} (cm);counts", kTH1D, {axisR}, true); + fRegistry.add((path + "hR2").c_str(), "R_{conv,2};R_{2} (cm);counts", kTH1D, {axisR}, true); + fRegistry.add((path + "hDeltaR").c_str(), "|R_{1}-R_{2}|;|R_{1}-R_{2}| (cm);counts", kTH1D, {axisDeltaR}, true); + fRegistry.add((path + "hDeltaZ").c_str(), "#Delta z;#Delta z (cm);counts", kTH1D, {axisDeltaZ}, true); + fRegistry.add((path + "hDeltaRxy").c_str(), "#Delta r_{xy};#Delta r_{xy} (cm);counts", kTH1D, {axisDeltaRxy}, true); + fRegistry.add((path + "hDeltaR3D").c_str(), "|#vec{r}_{1}-#vec{r}_{2}|;#Delta r_{3D} (cm);counts", kTH1D, {axisDeltaR3D}, true); + + // ── 1D: event-level ───────────────────────────────────────────────────── + fRegistry.add((path + "hCent").c_str(), "centrality;centrality (%);counts", kTH1D, {axisCentQA}, true); + fRegistry.add((path + "hOccupancy").c_str(), "occupancy;occupancy;counts", kTH1D, {axisOccupancy}, true); + + // ── 2D: angular ───────────────────────────────────────────────────────── + fRegistry.add((path + "hDEtaDPhi").c_str(), "#Delta#eta vs #Delta#phi;#Delta#eta;#Delta#phi (rad)", kTH2D, {axisDeltaEta, axisDeltaPhi}, true); + fRegistry.add((path + "hDeltaEtaVsPairEta").c_str(), "#Delta#eta vs #LT#eta#GT_{pair};#LT#eta#GT_{pair};#Delta#eta", kTH2D, {axisEta, axisDeltaEta}, true); + + // ── 2D: geometry ──────────────────────────────────────────────────────── + fRegistry.add((path + "hR1VsR2").c_str(), "R_{1} vs R_{2};R_{1} (cm);R_{2} (cm)", kTH2D, {axisR, axisR}, true); + fRegistry.add((path + "hDeltaRVsDeltaZ").c_str(), "|R_{1}-R_{2}| vs #Delta z;|R_{1}-R_{2}| (cm);#Delta z (cm)", kTH2D, {axisDeltaR, axisDeltaZ}, true); + + // ── 2D: geometry vs kT ────────────────────────────────────────────────── + // Note: hDeltaRVsQinv, hDeltaZVsQinv live in FullRange/ (always filled, full q range) + fRegistry.add((path + "hDeltaRVsKt").c_str(), "|R_{1}-R_{2}| vs k_{T};k_{T} (GeV/c);|R_{1}-R_{2}| (cm)", kTH2D, {axisKt, axisDeltaR}, true); + fRegistry.add((path + "hDeltaZVsKt").c_str(), "#Delta z vs k_{T};k_{T} (GeV/c);#Delta z (cm)", kTH2D, {axisKt, axisDeltaZ}, true); + + // ── 2D: angular vs geometry ───────────────────────────────────────────── + fRegistry.add((path + "hDeltaPhiVsDeltaR").c_str(), "#Delta#phi vs |R_{1}-R_{2}|;|R_{1}-R_{2}| (cm);#Delta#phi (rad)", kTH2D, {axisDeltaR, axisDeltaPhi}, true); + fRegistry.add((path + "hDeltaEtaVsDeltaR").c_str(), "#Delta#eta vs |R_{1}-R_{2}|;|R_{1}-R_{2}| (cm);#Delta#eta", kTH2D, {axisDeltaR, axisDeltaEta}, true); + fRegistry.add((path + "hDeltaPhiVsDeltaZ").c_str(), "#Delta#phi vs #Delta z;#Delta z (cm);#Delta#phi (rad)", kTH2D, {axisDeltaZ, axisDeltaPhi}, true); + fRegistry.add((path + "hDeltaEtaVsDeltaZ").c_str(), "#Delta#eta vs #Delta z;#Delta z (cm);#Delta#eta", kTH2D, {axisDeltaZ, axisDeltaEta}, true); + + // ── 2D: vs event properties ───────────────────────────────────────────── + fRegistry.add((path + "hDeltaRVsCent").c_str(), "|R_{1}-R_{2}| vs centrality;centrality (%);|R_{1}-R_{2}| (cm)", kTH2D, {axisCentQA, axisDeltaR}, true); + fRegistry.add((path + "hDeltaRVsOccupancy").c_str(), "|R_{1}-R_{2}| vs occupancy;occupancy;|R_{1}-R_{2}| (cm)", kTH2D, {axisOccupancy, axisDeltaR}, true); + + fRegistry.add((path + "hSparseDEtaDPhiCent").c_str(), + "#Delta#eta,#Delta#phi,centrality", + kTHnSparseD, {axisDeltaEta, axisDeltaPhi, axisCentQA}, true); + fRegistry.add((path + "hSparseDEtaDPhiOcc").c_str(), + "#Delta#eta,#Delta#phi,occupancy", + kTHnSparseD, {axisDeltaEta, axisDeltaPhi, axisOccupancy}, true); + fRegistry.add((path + "hSparseDEtaDPhiKt").c_str(), + "#Delta#eta,#Delta#phi,k_{T}", + kTHnSparseD, {axisDeltaEta, axisDeltaPhi, axisKt}, true); + fRegistry.add((path + "hSparseDeltaRDeltaZKt").c_str(), + "|R_{1}-R_{2}|,#Delta z,k_{T}", + kTHnSparseD, {axisDeltaR, axisDeltaZ, axisKt}, true); } void addhistograms() { static constexpr std::string_view det[6] = {"FT0M", "FT0A", "FT0C", "BTot", "BPos", "BNeg"}; - fRegistry.add("Event/before/hEP2_CentFT0C_forMix", Form("2nd harmonics EP for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", det[cfgEP2EstimatorForMix].data()), kTH2F, {{110, 0, 110}, {180, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}}, false); - fRegistry.add("Event/after/hEP2_CentFT0C_forMix", Form("2nd harmonics EP for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", det[cfgEP2EstimatorForMix].data()), kTH2F, {{110, 0, 110}, {180, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}}, false); - - // ── Single-photon QA ───────────────────────────────────────────────────── - fRegistry.add("SinglePhoton/hPt", "V0 photon p_{T};p_{T} (GeV/c);counts", kTH1F, {axisPt}, true); - fRegistry.add("SinglePhoton/hEta", "V0 photon #eta;#eta;counts", kTH1F, {axisEta}, true); - fRegistry.add("SinglePhoton/hPhi", "V0 photon #phi;#phi (rad);counts", kTH1F, {axisPhi}, true); - fRegistry.add("SinglePhoton/hEtaVsPhi", "V0 photon acceptance;#phi (rad);#eta", kTH2F, {axisPhi, axisEta}, true); - - // ── HBT physics ────────────────────────────────────────────────────────── + fRegistry.add("Event/before/hEP2_CentFT0C_forMix", + Form("2nd harmonics EP for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", det[cfgEP2EstimatorForMix].data()), + kTH2D, {{110, 0, 110}, {180, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}}, false); + fRegistry.add("Event/after/hEP2_CentFT0C_forMix", + Form("2nd harmonics EP for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", det[cfgEP2EstimatorForMix].data()), + kTH2D, {{110, 0, 110}, {180, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}}, false); + + addSinglePhotonQAHistogramsForStep("SinglePhoton/Before/"); + addSinglePhotonQAHistogramsForStep("SinglePhoton/AfterDRCosOA/"); + addSinglePhotonQAHistogramsForStep("SinglePhoton/AfterRZ/"); + addSinglePhotonQAHistogramsForStep("SinglePhoton/AfterEllipse/"); + + // ── HBT correlation functions ───────────────────────────────────────────── if (cfgDo3D) { - fRegistry.add("Pair/same/hs_3d", "diphoton correlation 3D LCMS", kTHnSparseD, {axisQout, axisQside, axisQlong, axisKt}, true); + fRegistry.add("Pair/same/CF_3D", "diphoton correlation 3D LCMS", + kTHnSparseD, {axisQout, axisQside, axisQlong, axisKt}, true); + if (cfgDo2D) { + fRegistry.add("Pair/same/CF_2D", "diphoton correlation 2D (qout,qinv)", + kTHnSparseD, {axisQout, axisQinv, axisKt}, true); + } } else { if (cfgUseLCMS) { - fRegistry.add("Pair/same/hs_1d", "diphoton correlation 1D LCMS", kTHnSparseD, {axisQabsLcms, axisKt}, true); + fRegistry.add("Pair/same/CF_1D", "diphoton correlation 1D LCMS", kTH2D, {axisQabsLcms, axisKt}, true); } else { - fRegistry.add("Pair/same/hs_1d", "diphoton correlation 1D", kTHnSparseD, {axisQinv, axisKt}, true); + fRegistry.add("Pair/same/CF_1D", "diphoton correlation 1D (qinv)", kTH2D, {axisQinv, axisKt}, true); } } - fRegistry.add("Pair/same/hDeltaRCosOA", "distance between 2 conversion points;#Deltar/cos(#theta_{op}/2) (cm)", kTH1D, {{100, 0, 100}}, true); + fRegistry.add("Pair/same/hDeltaRCosOA", + "distance between 2 conversion points / cos(#theta_{op}/2);#Delta r / cos(#theta_{op}/2) (cm);counts", + kTH1D, {{100, 0, 100}}, true); + // ── QA steps (same-event; mix-event cloned below) ───────────────────────── addQAHistogramsForStep("Pair/same/QA/Before/"); - addQAHistogramsForStep("Pair/same/QA/After/"); + addQAHistogramsForStep("Pair/same/QA/AfterDRCosOA/"); + addQAHistogramsForStep("Pair/same/QA/AfterRZ/"); + addQAHistogramsForStep("Pair/same/QA/AfterEllipse/"); + + // ── MC truth histograms (same-event; mix-event cloned below) ───────────── + addMCHistograms(); + + // ── Full-range histograms: always filled, qinv as axis ─────────────────── + addFullRangeHistograms("Pair/same/FullRange/"); + // Clone all Pair/same/ histograms to Pair/mix/ fRegistry.addClone("Pair/same/", "Pair/mix/"); } - // --------------------------------------------------------------------------- - // DefineEMEventCut / DefinePCMCut - // --------------------------------------------------------------------------- + // ─── DefineEMEventCut ────────────────────────────────────────────────────── void DefineEMEventCut() { @@ -451,6 +704,8 @@ struct photonhbt { fEMEventCut.SetRequireGoodITSLayersAll(eventcuts.cfgRequireGoodITSLayersAll); } + // ─── DefinePCMCut ────────────────────────────────────────────────────────── + void DefinePCMCut() { fV0PhotonCut = V0PhotonCut("fV0PhotonCut", "fV0PhotonCut"); @@ -478,8 +733,37 @@ struct photonhbt { fV0PhotonCut.SetRequireTPConly(pcmcuts.cfgRequireV0WithTPCOnly); } + /// step_id: 0 = Before, 1 = AfterDRCosOA, 2 = AfterRZ, 3 = AfterEllipse + template + static constexpr const char* singlePhotonQAPrefix() + { + if constexpr (step_id == 0) + return "SinglePhoton/Before/"; + if constexpr (step_id == 1) + return "SinglePhoton/AfterDRCosOA/"; + if constexpr (step_id == 2) + return "SinglePhoton/AfterRZ/"; + return "SinglePhoton/AfterEllipse/"; + } + + template + inline void fillSinglePhotonQAStep(TPhoton const& g) + { + if (!qaflags.doSinglePhotonQa) + return; + constexpr auto base = singlePhotonQAPrefix(); + const float r = std::sqrt(g.vx() * g.vx() + g.vy() * g.vy()); + fRegistry.fill(HIST(base) + HIST("hPt"), g.pt()); + fRegistry.fill(HIST(base) + HIST("hEta"), g.eta()); + fRegistry.fill(HIST(base) + HIST("hPhi"), g.phi()); + fRegistry.fill(HIST(base) + HIST("hEtaVsPhi"), g.phi(), g.eta()); + fRegistry.fill(HIST(base) + HIST("hR"), r); + fRegistry.fill(HIST(base) + HIST("hZConv"), g.vz()); + fRegistry.fill(HIST(base) + HIST("hRVsZConv"), g.vz(), r); + } + template - void fillPairHistogram(TCollision const&, + void fillPairHistogram(TCollision const& /*collision*/, ROOT::Math::PtEtaPhiMVector v1, ROOT::Math::PtEtaPhiMVector v2, float weight = 1.f) @@ -504,60 +788,497 @@ struct photonhbt { float qlong_lcms = q3_lcms.Dot(uv_long); if (cfgDo3D) { - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_3d"), + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("CF_3D"), std::fabs(qout_lcms), std::fabs(qside_lcms), std::fabs(qlong_lcms), kt, weight); + if (cfgDo2D) { + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("CF_2D"), + std::fabs(qout_lcms), std::fabs(qinv), kt, weight); + } } else { - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_1d"), + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("CF_1D"), cfgUseLCMS ? qabs_lcms : qinv, kt, weight); } } - template - inline void fillPairQAStep(float deta, float dphi, float pairEta, - float cosTheta, float openingAngle) + template + void fillPairHistogramMC(TCollision const& /*collision*/, + ROOT::Math::PtEtaPhiMVector v1, + ROOT::Math::PtEtaPhiMVector v2, + float weight = 1.f) + { + + float rndm = std::pow(-1, dist01(engine) % 2); + auto k12 = 0.5 * (v1 + v2); + float kt = k12.Pt(); + float qinv = -(((v1 - v2) * rndm).M()); + + ROOT::Math::XYZVector uv_out(k12.Px() / k12.Pt(), k12.Py() / k12.Pt(), 0); + ROOT::Math::XYZVector uv_long(0, 0, 1); + ROOT::Math::XYZVector uv_side = uv_out.Cross(uv_long); + + ROOT::Math::PxPyPzEVector v1c(v1), v2c(v2); + float beta_z = (v1 + v2).Beta() * std::cos((v1 + v2).Theta()); + ROOT::Math::Boost bst_z(0, 0, -beta_z); + auto q12_lcms = bst_z((v1c - v2c) * rndm); + auto q3_lcms = q12_lcms.Vect(); + float qabs_lcms = q3_lcms.R(); + float qout_lcms = q3_lcms.Dot(uv_out); + float qside_lcms = q3_lcms.Dot(uv_side); + float qlong_lcms = q3_lcms.Dot(uv_long); + + constexpr auto mcDir = []() constexpr -> const char* { + if constexpr (ev_id == 0) { + if constexpr (TruthT == PairTruthType::TrueTrueDistinct) + return "Pair/same/MC/TrueTrueDistinct/"; + if constexpr (TruthT == PairTruthType::TrueTrueSamePhoton) + return "Pair/same/MC/TrueTrueSamePhoton/"; + if constexpr (TruthT == PairTruthType::SharedMcLeg) + return "Pair/same/MC/SharedMcLeg/"; + if constexpr (TruthT == PairTruthType::TrueFake) + return "Pair/same/MC/TrueFake/"; + if constexpr (TruthT == PairTruthType::FakeFake) + return "Pair/same/MC/FakeFake/"; + return "Pair/same/MC/Pi0Daughters/"; + } else { + if constexpr (TruthT == PairTruthType::TrueTrueDistinct) + return "Pair/mix/MC/TrueTrueDistinct/"; + if constexpr (TruthT == PairTruthType::TrueTrueSamePhoton) + return "Pair/mix/MC/TrueTrueSamePhoton/"; + if constexpr (TruthT == PairTruthType::SharedMcLeg) + return "Pair/mix/MC/SharedMcLeg/"; + if constexpr (TruthT == PairTruthType::TrueFake) + return "Pair/mix/MC/TrueFake/"; + if constexpr (TruthT == PairTruthType::FakeFake) + return "Pair/mix/MC/FakeFake/"; + return "Pair/mix/MC/Pi0Daughters/"; + } + }(); + + if (cfgDo3D) { + fRegistry.fill(HIST(mcDir) + HIST("CF_3D"), + std::fabs(qout_lcms), std::fabs(qside_lcms), std::fabs(qlong_lcms), kt, weight); + if (cfgDo2D) { + fRegistry.fill(HIST(mcDir) + HIST("CF_2D"), + std::fabs(qout_lcms), std::fabs(qinv), kt, weight); + } + } else { + fRegistry.fill(HIST(mcDir) + HIST("CF_1D"), + cfgUseLCMS ? qabs_lcms : qinv, kt, weight); + } + } + + template + PairQAObservables buildPairQAObservables(TG1 const& g1, TG2 const& g2) + { + PairQAObservables o{}; + + o.x1 = g1.vx(); + o.y1 = g1.vy(); + o.z1 = g1.vz(); + o.x2 = g2.vx(); + o.y2 = g2.vy(); + o.z2 = g2.vz(); + + o.r1 = std::sqrt(o.x1 * o.x1 + o.y1 * o.y1); + o.r2 = std::sqrt(o.x2 * o.x2 + o.y2 * o.y2); + + o.dx = o.x1 - o.x2; + o.dy = o.y1 - o.y2; + o.dz = o.z1 - o.z2; + + o.deltaR = std::fabs(o.r1 - o.r2); + o.deltaZ = o.dz; + o.deltaRxy = std::sqrt(o.dx * o.dx + o.dy * o.dy); + o.deltaR3D = std::sqrt(o.dx * o.dx + o.dy * o.dy + o.dz * o.dz); + + ROOT::Math::XYZVector cp1(o.x1, o.y1, o.z1); + ROOT::Math::XYZVector cp2(o.x2, o.y2, o.z2); + const float mag1 = std::sqrt(cp1.Mag2()); + const float mag2 = std::sqrt(cp2.Mag2()); + if (mag1 < 1e-12f || mag2 < 1e-12f) { + o.valid = false; + return o; + } + float cosPA = static_cast(cp1.Dot(cp2) / (mag1 * mag2)); + cosPA = std::clamp(cosPA, -1.f, 1.f); + o.opa = std::acos(cosPA); + o2::math_utils::bringTo02Pi(o.opa); + if (o.opa > o2::constants::math::PI) + o.opa -= o2::constants::math::PI; + o.cosOA = std::cos(o.opa / 2.f); + o.drOverCosOA = (std::fabs(o.cosOA) < 1e-12f) ? 1e12f : (o.deltaR3D / o.cosOA); + + o.v1 = ROOT::Math::PtEtaPhiMVector(g1.pt(), g1.eta(), g1.phi(), 0.f); + o.v2 = ROOT::Math::PtEtaPhiMVector(g2.pt(), g2.eta(), g2.phi(), 0.f); + o.k12 = 0.5f * (o.v1 + o.v2); + + o.deta = g1.eta() - g2.eta(); + o.dphi = RecoDecay::constrainAngle(g1.phi() - g2.phi(), -o2::constants::math::PI); // dphi in [-pi, pi] + o.pairEta = 0.5f * (g1.eta() + g2.eta()); + o.pairPhi = RecoDecay::constrainAngle(o.k12.Phi(), 0.f); // pair phi in [0, 2pi] — matches axisPhi + o.kt = o.k12.Pt(); + o.qinv = std::fabs((o.v1 - o.v2).M()); + o.cosTheta = std::fabs(computeCosTheta(o.v1, o.v2)); + o.openingAngle = o.opa; + + return o; + } + + template + inline void fillPairQAStep(PairQAObservables const& o, float cent, float occupancy) { if (!qaflags.doPairQa) return; + constexpr auto base = qaPrefix(); + + // 1D: kinematics + fRegistry.fill(HIST(base) + HIST("hPairEta"), o.pairEta); + fRegistry.fill(HIST(base) + HIST("hPairPhi"), o.pairPhi); + fRegistry.fill(HIST(base) + HIST("hPairKt"), o.kt); + fRegistry.fill(HIST(base) + HIST("hQinv"), o.qinv); + + // 1D: angular + fRegistry.fill(HIST(base) + HIST("hDeltaEta"), o.deta); + fRegistry.fill(HIST(base) + HIST("hDeltaPhi"), o.dphi); + fRegistry.fill(HIST(base) + HIST("hCosTheta"), o.cosTheta); + fRegistry.fill(HIST(base) + HIST("hOpeningAngle"), o.openingAngle); + + // 1D: geometry + fRegistry.fill(HIST(base) + HIST("hR1"), o.r1); + fRegistry.fill(HIST(base) + HIST("hR2"), o.r2); + fRegistry.fill(HIST(base) + HIST("hDeltaR"), o.deltaR); + fRegistry.fill(HIST(base) + HIST("hDeltaZ"), o.deltaZ); + fRegistry.fill(HIST(base) + HIST("hDeltaRxy"), o.deltaRxy); + fRegistry.fill(HIST(base) + HIST("hDeltaR3D"), o.deltaR3D); + + // 1D: event + fRegistry.fill(HIST(base) + HIST("hCent"), cent); + fRegistry.fill(HIST(base) + HIST("hOccupancy"), occupancy); + + // 1D: ellipse value (diagnostic, conditional on cut being configured) const float sE = ggpaircuts.cfgEllipseSigEta.value; const float sP = ggpaircuts.cfgEllipseSigPhi.value; + if (sE > 1e-9f && sP > 1e-9f) { + const float ellipseVal = (o.deta / sE) * (o.deta / sE) + (o.dphi / sP) * (o.dphi / sP); + fRegistry.fill(HIST(base) + HIST("hEllipseVal"), ellipseVal); + } - if constexpr (ev_id == 0 && IsBefore) { - fRegistry.fill(HIST("Pair/same/QA/Before/hDeltaEta"), deta); - fRegistry.fill(HIST("Pair/same/QA/Before/hDeltaPhi"), dphi); - fRegistry.fill(HIST("Pair/same/QA/Before/hDEtaDPhi"), deta, dphi); - fRegistry.fill(HIST("Pair/same/QA/Before/hDeltaEtaVsPairEta"), pairEta, deta); - fRegistry.fill(HIST("Pair/same/QA/Before/hCosTheta"), cosTheta); - fRegistry.fill(HIST("Pair/same/QA/Before/hOpeningAngle"), openingAngle); - if (sE > 1e-9f && sP > 1e-9f) - fRegistry.fill(HIST("Pair/same/QA/Before/hEllipseVal"), (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP)); - } else if constexpr (ev_id == 0 && !IsBefore) { - fRegistry.fill(HIST("Pair/same/QA/After/hDeltaEta"), deta); - fRegistry.fill(HIST("Pair/same/QA/After/hDeltaPhi"), dphi); - fRegistry.fill(HIST("Pair/same/QA/After/hDEtaDPhi"), deta, dphi); - fRegistry.fill(HIST("Pair/same/QA/After/hDeltaEtaVsPairEta"), pairEta, deta); - fRegistry.fill(HIST("Pair/same/QA/After/hCosTheta"), cosTheta); - fRegistry.fill(HIST("Pair/same/QA/After/hOpeningAngle"), openingAngle); - if (sE > 1e-9f && sP > 1e-9f) - fRegistry.fill(HIST("Pair/same/QA/After/hEllipseVal"), (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP)); - } else if constexpr (ev_id == 1 && IsBefore) { - fRegistry.fill(HIST("Pair/mix/QA/Before/hDeltaEta"), deta); - fRegistry.fill(HIST("Pair/mix/QA/Before/hDeltaPhi"), dphi); - fRegistry.fill(HIST("Pair/mix/QA/Before/hDEtaDPhi"), deta, dphi); - fRegistry.fill(HIST("Pair/mix/QA/Before/hDeltaEtaVsPairEta"), pairEta, deta); - fRegistry.fill(HIST("Pair/mix/QA/Before/hCosTheta"), cosTheta); - fRegistry.fill(HIST("Pair/mix/QA/Before/hOpeningAngle"), openingAngle); - if (sE > 1e-9f && sP > 1e-9f) - fRegistry.fill(HIST("Pair/mix/QA/Before/hEllipseVal"), (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP)); - } else { - fRegistry.fill(HIST("Pair/mix/QA/After/hDeltaEta"), deta); - fRegistry.fill(HIST("Pair/mix/QA/After/hDeltaPhi"), dphi); - fRegistry.fill(HIST("Pair/mix/QA/After/hDEtaDPhi"), deta, dphi); - fRegistry.fill(HIST("Pair/mix/QA/After/hDeltaEtaVsPairEta"), pairEta, deta); - fRegistry.fill(HIST("Pair/mix/QA/After/hCosTheta"), cosTheta); - fRegistry.fill(HIST("Pair/mix/QA/After/hOpeningAngle"), openingAngle); - if (sE > 1e-9f && sP > 1e-9f) - fRegistry.fill(HIST("Pair/mix/QA/After/hEllipseVal"), (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP)); + // 2D: angular + fRegistry.fill(HIST(base) + HIST("hDEtaDPhi"), o.deta, o.dphi); + fRegistry.fill(HIST(base) + HIST("hDeltaEtaVsPairEta"), o.pairEta, o.deta); + + // 2D: geometry + fRegistry.fill(HIST(base) + HIST("hR1VsR2"), o.r1, o.r2); + fRegistry.fill(HIST(base) + HIST("hDeltaRVsDeltaZ"), o.deltaR, o.deltaZ); + + // 2D: geometry vs kT (qinv variants live in FullRange/) + fRegistry.fill(HIST(base) + HIST("hDeltaRVsKt"), o.kt, o.deltaR); + fRegistry.fill(HIST(base) + HIST("hDeltaZVsKt"), o.kt, o.deltaZ); + + // 2D: angular vs geometry + fRegistry.fill(HIST(base) + HIST("hDeltaPhiVsDeltaR"), o.deltaR, o.dphi); + fRegistry.fill(HIST(base) + HIST("hDeltaEtaVsDeltaR"), o.deltaR, o.deta); + fRegistry.fill(HIST(base) + HIST("hDeltaPhiVsDeltaZ"), o.deltaZ, o.dphi); + fRegistry.fill(HIST(base) + HIST("hDeltaEtaVsDeltaZ"), o.deltaZ, o.deta); + + // 2D: vs event properties (qinv variants live in FullRange/) + fRegistry.fill(HIST(base) + HIST("hDeltaRVsCent"), cent, o.deltaR); + fRegistry.fill(HIST(base) + HIST("hDeltaRVsOccupancy"), occupancy, o.deltaR); + + // THnSparse (hSparseDeltaRDeltaZQinv lives in FullRange/) + fRegistry.fill(HIST(base) + HIST("hSparseDEtaDPhiCent"), o.deta, o.dphi, cent); + fRegistry.fill(HIST(base) + HIST("hSparseDEtaDPhiOcc"), o.deta, o.dphi, occupancy); + fRegistry.fill(HIST(base) + HIST("hSparseDEtaDPhiKt"), o.deta, o.dphi, o.kt); + fRegistry.fill(HIST(base) + HIST("hSparseDeltaRDeltaZKt"), o.deltaR, o.deltaZ, o.kt); + } + + template + static PhotonMCInfo buildPhotonMCInfo(TPhoton const& g, + TMCParticles const& mcParticles) + { + PhotonMCInfo info{}; + + const auto pos = g.template posTrack_as(); + const auto neg = g.template negTrack_as(); + + // PWGEM uses emmcparticle, not the standard mcParticle accessor + if (!pos.has_emmcparticle() || !neg.has_emmcparticle()) + return info; + + info.hasMC = true; + info.mcPosId = pos.emmcparticleId(); + info.mcNegId = neg.emmcparticleId(); + + const auto mcPos = pos.template emmcparticle_as(); + const auto mcNeg = neg.template emmcparticle_as(); + + if (!mcPos.has_mothers() || !mcNeg.has_mothers()) + return info; + + const int mothIdPos = mcPos.mothersIds()[0]; + const int mothIdNeg = mcNeg.mothersIds()[0]; + if (mothIdPos != mothIdNeg) + return info; + + info.sameMother = true; + info.motherId = mothIdPos; + + const auto mother = mcParticles.iteratorAt(mothIdPos); + info.motherPdg = mother.pdgCode(); + info.isTruePhoton = (info.motherPdg == 22); + info.isPhysicalPrimary = mother.isPhysicalPrimary(); + + return info; + } + + static PairTruthType classifyPairTruth(PhotonMCInfo const& m1, + PhotonMCInfo const& m2) + { + const bool t1 = m1.hasMC && m1.sameMother && m1.isTruePhoton; + const bool t2 = m2.hasMC && m2.sameMother && m2.isTruePhoton; + + if (m1.hasMC && m2.hasMC) { + if ((m1.mcPosId >= 0 && (m1.mcPosId == m2.mcPosId || m1.mcPosId == m2.mcNegId)) || + (m1.mcNegId >= 0 && (m1.mcNegId == m2.mcPosId || m1.mcNegId == m2.mcNegId))) + return PairTruthType::SharedMcLeg; + } + + if (!t1 && !t2) + return PairTruthType::FakeFake; + if (t1 != t2) + return PairTruthType::TrueFake; + + // Both are true photons — same or different MC photon? + if (m1.motherId >= 0 && m1.motherId == m2.motherId) + return PairTruthType::TrueTrueSamePhoton; + + return PairTruthType::TrueTrueDistinct; + } + + template + static bool isPi0DaughterPair(PhotonMCInfo const& m1, + PhotonMCInfo const& m2, + TMCParticles const& mcParticles) + { + if (!m1.isTruePhoton || !m2.isTruePhoton) + return false; + if (m1.motherId < 0 || m2.motherId < 0) + return false; + // The photons themselves must have the same grandmother = pi0 + const auto ph1 = mcParticles.iteratorAt(m1.motherId); + const auto ph2 = mcParticles.iteratorAt(m2.motherId); + if (!ph1.has_mothers() || !ph2.has_mothers()) + return false; + const int gm1 = ph1.mothersIds()[0]; + const int gm2 = ph2.mothersIds()[0]; + if (gm1 != gm2) + return false; + return (std::abs(mcParticles.iteratorAt(gm1).pdgCode()) == 111); + } + + static constexpr std::string_view pairTruthLabel(PairTruthType t) + { + switch (t) { + case PairTruthType::TrueTrueDistinct: + return "TrueTrueDistinct/"; + case PairTruthType::TrueTrueSamePhoton: + return "TrueTrueSamePhoton/"; + case PairTruthType::SharedMcLeg: + return "SharedMcLeg/"; + case PairTruthType::TrueFake: + return "TrueFake/"; + case PairTruthType::FakeFake: + return "FakeFake/"; + case PairTruthType::Pi0Daughters: + return "Pi0Daughters/"; + default: + return "Unknown/"; + } + } + + void addMCHistograms() + { + const AxisSpec axisTruthType{{0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5}, "truth type (1=TrueTrueDistinct,2=TrueTrueSamePhoton,3=SharedMcLeg,4=TrueFake,5=FakeFake,6=Pi0Daughters)"}; + + static constexpr std::array kTypes = { + "TrueTrueDistinct/", + "TrueTrueSamePhoton/", + "SharedMcLeg/", + "TrueFake/", + "FakeFake/", + "Pi0Daughters/"}; + + for (const auto& label : kTypes) { + const std::string base = std::string("Pair/same/MC/") + std::string(label); + + if (cfgDo3D) { + fRegistry.add((base + "CF_3D").c_str(), "MC CF 3D LCMS", kTHnSparseD, {axisQout, axisQside, axisQlong, axisKt}, true); + if (cfgDo2D) { + fRegistry.add((base + "CF_2D").c_str(), "MC CF 2D", kTHnSparseD, {axisQout, axisQinv, axisKt}, true); + } + } else { + if (cfgUseLCMS) { + fRegistry.add((base + "CF_1D").c_str(), "MC CF 1D LCMS", kTH2D, {axisQabsLcms, axisKt}, true); + } else { + fRegistry.add((base + "CF_1D").c_str(), "MC CF 1D (qinv)", kTH2D, {axisQinv, axisKt}, true); + } + } + + fRegistry.add((base + "hQinv").c_str(), "q_{inv};q_{inv} (GeV/c);counts", kTH1D, {axisQinv}, true); + fRegistry.add((base + "hDeltaEta").c_str(), "#Delta#eta;#Delta#eta;counts", kTH1D, {axisDeltaEta}, true); + fRegistry.add((base + "hDeltaPhi").c_str(), "#Delta#phi;#Delta#phi (rad);counts", kTH1D, {axisDeltaPhi}, true); + fRegistry.add((base + "hDEtaDPhi").c_str(), "#Delta#eta vs #Delta#phi;#Delta#eta;#Delta#phi", kTH2D, {axisDeltaEta, axisDeltaPhi}, true); + fRegistry.add((base + "hDeltaR").c_str(), "|R_{1}-R_{2}|;|R_{1}-R_{2}| (cm);counts", kTH1D, {axisDeltaR}, true); + fRegistry.add((base + "hDeltaZ").c_str(), "#Delta z;#Delta z (cm);counts", kTH1D, {axisDeltaZ}, true); + fRegistry.add((base + "hDeltaR3D").c_str(), "#Delta r_{3D};#Delta r_{3D} (cm);counts", kTH1D, {axisDeltaR3D}, true); + fRegistry.add((base + "hKt").c_str(), "k_{T};k_{T} (GeV/c);counts", kTH1D, {axisKt}, true); + fRegistry.add((base + "hDeltaRVsQinv").c_str(), "|R_{1}-R_{2}| vs q_{inv};q_{inv} (GeV/c);|R_{1}-R_{2}| (cm)", kTH2D, {axisQinv, axisDeltaR}, true); + fRegistry.add((base + "hDeltaZVsQinv").c_str(), "#Delta z vs q_{inv};q_{inv} (GeV/c);#Delta z (cm)", kTH2D, {axisQinv, axisDeltaZ}, true); + fRegistry.add((base + "hDeltaR3DVsQinv").c_str(), "#Delta r_{3D} vs q_{inv};q_{inv} (GeV/c);#Delta r_{3D} (cm)", kTH2D, {axisQinv, axisDeltaR3D}, true); + fRegistry.add((base + "hDEtaDPhiVsQinv").c_str(), "#Delta#eta vs #Delta#phi vs q_{inv};#Delta#eta;#Delta#phi;q_{inv}", kTHnSparseD, {axisDeltaEta, axisDeltaPhi, axisQinv}, true); + fRegistry.add((base + "hSparseDeltaRDeltaZQinv").c_str(), "|R_{1}-R_{2}|,#Delta z,q_{inv};|R_{1}-R_{2}| (cm);#Delta z (cm);q_{inv} (GeV/c)", kTHnSparseD, {axisDeltaR, axisDeltaZ, axisQinv}, true); + } + + fRegistry.add("Pair/same/MC/hTruthTypeVsQinv", "truth type vs q_{inv};q_{inv} (GeV/c);truth type", kTH2D, {axisQinv, axisTruthType}, true); + fRegistry.add("Pair/same/MC/hTruthTypeVsKt", "truth type vs k_{T};k_{T} (GeV/c);truth type", kTH2D, {axisKt, axisTruthType}, true); + } + + template + inline void fillMCPairQATyped(PairQAObservables const& obs) + { + constexpr auto base = []() constexpr -> const char* { + if constexpr (!IsMix) { + if constexpr (TruthT == PairTruthType::TrueTrueDistinct) + return "Pair/same/MC/TrueTrueDistinct/"; + if constexpr (TruthT == PairTruthType::TrueTrueSamePhoton) + return "Pair/same/MC/TrueTrueSamePhoton/"; + if constexpr (TruthT == PairTruthType::SharedMcLeg) + return "Pair/same/MC/SharedMcLeg/"; + if constexpr (TruthT == PairTruthType::TrueFake) + return "Pair/same/MC/TrueFake/"; + if constexpr (TruthT == PairTruthType::FakeFake) + return "Pair/same/MC/FakeFake/"; + return "Pair/same/MC/Pi0Daughters/"; + } else { + if constexpr (TruthT == PairTruthType::TrueTrueDistinct) + return "Pair/mix/MC/TrueTrueDistinct/"; + if constexpr (TruthT == PairTruthType::TrueTrueSamePhoton) + return "Pair/mix/MC/TrueTrueSamePhoton/"; + if constexpr (TruthT == PairTruthType::SharedMcLeg) + return "Pair/mix/MC/SharedMcLeg/"; + if constexpr (TruthT == PairTruthType::TrueFake) + return "Pair/mix/MC/TrueFake/"; + if constexpr (TruthT == PairTruthType::FakeFake) + return "Pair/mix/MC/FakeFake/"; + return "Pair/mix/MC/Pi0Daughters/"; + } + }(); + + fRegistry.fill(HIST(base) + HIST("hQinv"), obs.qinv); + fRegistry.fill(HIST(base) + HIST("hDeltaEta"), obs.deta); + fRegistry.fill(HIST(base) + HIST("hDeltaPhi"), obs.dphi); + fRegistry.fill(HIST(base) + HIST("hDEtaDPhi"), obs.deta, obs.dphi); + fRegistry.fill(HIST(base) + HIST("hDeltaR"), obs.deltaR); + fRegistry.fill(HIST(base) + HIST("hDeltaZ"), obs.deltaZ); + fRegistry.fill(HIST(base) + HIST("hDeltaR3D"), obs.deltaR3D); + fRegistry.fill(HIST(base) + HIST("hKt"), obs.kt); + + constexpr auto summaryDir = IsMix ? "Pair/mix/MC/" : "Pair/same/MC/"; + const int typeIdx = static_cast(TruthT); + fRegistry.fill(HIST(summaryDir) + HIST("hTruthTypeVsQinv"), obs.qinv, typeIdx); + fRegistry.fill(HIST(summaryDir) + HIST("hTruthTypeVsKt"), obs.kt, typeIdx); + } + + template + inline void fillMCPairQA(PairTruthType truthType, PairQAObservables const& obs) + { + switch (truthType) { + case PairTruthType::TrueTrueDistinct: + fillMCPairQATyped(obs); + break; + case PairTruthType::TrueTrueSamePhoton: + fillMCPairQATyped(obs); + break; + case PairTruthType::SharedMcLeg: + fillMCPairQATyped(obs); + break; + case PairTruthType::TrueFake: + fillMCPairQATyped(obs); + break; + case PairTruthType::FakeFake: + fillMCPairQATyped(obs); + break; + case PairTruthType::Pi0Daughters: + fillMCPairQATyped(obs); + break; + default: + break; + } + } + + template + inline void fillMCPairQAFullRangeTyped(PairQAObservables const& obs) + { + constexpr auto base = []() constexpr -> const char* { + if constexpr (!IsMix) { + if constexpr (TruthT == PairTruthType::TrueTrueDistinct) + return "Pair/same/MC/TrueTrueDistinct/"; + if constexpr (TruthT == PairTruthType::TrueTrueSamePhoton) + return "Pair/same/MC/TrueTrueSamePhoton/"; + if constexpr (TruthT == PairTruthType::SharedMcLeg) + return "Pair/same/MC/SharedMcLeg/"; + if constexpr (TruthT == PairTruthType::TrueFake) + return "Pair/same/MC/TrueFake/"; + if constexpr (TruthT == PairTruthType::FakeFake) + return "Pair/same/MC/FakeFake/"; + return "Pair/same/MC/Pi0Daughters/"; + } else { + if constexpr (TruthT == PairTruthType::TrueTrueDistinct) + return "Pair/mix/MC/TrueTrueDistinct/"; + if constexpr (TruthT == PairTruthType::TrueTrueSamePhoton) + return "Pair/mix/MC/TrueTrueSamePhoton/"; + if constexpr (TruthT == PairTruthType::SharedMcLeg) + return "Pair/mix/MC/SharedMcLeg/"; + if constexpr (TruthT == PairTruthType::TrueFake) + return "Pair/mix/MC/TrueFake/"; + if constexpr (TruthT == PairTruthType::FakeFake) + return "Pair/mix/MC/FakeFake/"; + return "Pair/mix/MC/Pi0Daughters/"; + } + }(); + + fRegistry.fill(HIST(base) + HIST("hDeltaRVsQinv"), obs.qinv, obs.deltaR); + fRegistry.fill(HIST(base) + HIST("hDeltaZVsQinv"), obs.qinv, obs.deltaZ); + fRegistry.fill(HIST(base) + HIST("hDeltaR3DVsQinv"), obs.qinv, obs.deltaR3D); + fRegistry.fill(HIST(base) + HIST("hDEtaDPhiVsQinv"), obs.deta, obs.dphi, obs.qinv); + fRegistry.fill(HIST(base) + HIST("hSparseDeltaRDeltaZQinv"), obs.deltaR, obs.deltaZ, obs.qinv); + } + + template + inline void fillMCPairQAFullRange(PairTruthType truthType, PairQAObservables const& obs) + { + switch (truthType) { + case PairTruthType::TrueTrueDistinct: + fillMCPairQAFullRangeTyped(obs); + break; + case PairTruthType::TrueTrueSamePhoton: + fillMCPairQAFullRangeTyped(obs); + break; + case PairTruthType::SharedMcLeg: + fillMCPairQAFullRangeTyped(obs); + break; + case PairTruthType::TrueFake: + fillMCPairQAFullRangeTyped(obs); + break; + case PairTruthType::FakeFake: + fillMCPairQAFullRangeTyped(obs); + break; + case PairTruthType::Pi0Daughters: + fillMCPairQAFullRangeTyped(obs); + break; + default: + break; } } @@ -576,127 +1297,167 @@ struct photonhbt { initCCDB(collision); int ndiphoton = 0; + // ── Centrality selection ────────────────────────────────────────────── const float cent[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; if (cent[cfgCentEstimator] < cfgCentMin || cfgCentMax < cent[cfgCentEstimator]) continue; - const std::array epArr = {collision.ep2ft0m(), collision.ep2ft0a(), collision.ep2ft0c(), - collision.ep2fv0a(), collision.ep2btot(), collision.ep2bpos(), collision.ep2bneg()}; - float ep2 = epArr[cfgEP2EstimatorForMix]; + const std::array epArr = { + collision.ep2ft0m(), collision.ep2ft0a(), collision.ep2ft0c(), + collision.ep2fv0a(), collision.ep2btot(), collision.ep2bpos(), collision.ep2bneg()}; + const float ep2 = epArr[cfgEP2EstimatorForMix]; + // ── Event QA and event cut ──────────────────────────────────────────── fRegistry.fill(HIST("Event/before/hEP2_CentFT0C_forMix"), collision.centFT0C(), ep2); o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(&fRegistry, collision, 1.f); if (!fEMEventCut.IsSelected(collision)) continue; o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(&fRegistry, collision, 1.f); + fRegistry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted + fRegistry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted fRegistry.fill(HIST("Event/after/hEP2_CentFT0C_forMix"), collision.centFT0C(), ep2); - // Mixing-bin indices — uses static binOf helper: - int zbin = binOf(ztxBinEdges, collision.posZ()); - int centbin = binOf(centBinEdges, cent[cfgCentEstimator]); - int epbin = binOf(epBinEgdes, ep2); - int occbin = binOf(occBinEdges, - cfgOccupancyEstimator == 1 - ? static_cast(collision.trackOccupancyInTimeRange()) - : collision.ft0cOccupancyInTimeRange()); + // ── Event mixing bins ───────────────────────────────────────────────── + const float occupancy = (cfgOccupancyEstimator == 1) + ? static_cast(collision.trackOccupancyInTimeRange()) + : collision.ft0cOccupancyInTimeRange(); + const float centForQA = cent[cfgCentEstimator]; + + const int zbin = binOf(ztxBinEdges, collision.posZ()); + const int centbin = binOf(centBinEdges, centForQA); + const int epbin = binOf(epBinEgdes, ep2); + const int occbin = binOf(occBinEdges, occupancy); auto keyBin = std::make_tuple(zbin, centbin, epbin, occbin); auto keyDFCollision = std::make_pair(ndf, collision.globalIndex()); + // ── Slice photons for this collision ────────────────────────────────── auto photons1Coll = photons1.sliceBy(perCollision1, collision.globalIndex()); auto photons2Coll = photons2.sliceBy(perCollision2, collision.globalIndex()); - // ── Single-photon QA ───────────────────────────────────────────────── + // ── Single-photon QA if (qaflags.doSinglePhotonQa) { for (const auto& g : photons1Coll) { if (!cut1.template IsSelected(g)) continue; - fRegistry.fill(HIST("SinglePhoton/hPt"), g.pt()); - fRegistry.fill(HIST("SinglePhoton/hEta"), g.eta()); - fRegistry.fill(HIST("SinglePhoton/hPhi"), g.phi()); - fRegistry.fill(HIST("SinglePhoton/hEtaVsPhi"), g.phi(), g.eta()); + fillSinglePhotonQAStep<0>(g); } } + std::unordered_set photonIdsAfterDRCosOA; + std::unordered_set photonIdsAfterRZ; + std::unordered_set photonIdsAfterEllipse; + // ── Same-event pair loop ────────────────────────────────────────────── for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photons1Coll, photons2Coll))) { if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) continue; - auto pos1 = g1.template posTrack_as(); - auto ele1 = g1.template negTrack_as(); - auto pos2 = g2.template posTrack_as(); - auto ele2 = g2.template negTrack_as(); - if (pos1.trackId() == pos2.trackId() || ele1.trackId() == ele2.trackId()) + const auto pos1 = g1.template posTrack_as(); + const auto ele1 = g1.template negTrack_as(); + const auto pos2 = g2.template posTrack_as(); + const auto ele2 = g2.template negTrack_as(); + if (pos1.trackId() == pos2.trackId() || + pos1.trackId() == ele2.trackId() || + ele1.trackId() == pos2.trackId() || + ele1.trackId() == ele2.trackId()) + continue; + + auto obs = buildPairQAObservables(g1, g2); + if (!obs.valid) + continue; + + const bool doQA = passQinvQAGate(obs.qinv); + const bool doFullRange = passQinvFullRangeGate(obs.qinv); + + // ── QA: Before any pair cut ─────────────────────────────────────── + if (doQA) + fillPairQAStep<0, 0>(obs, centForQA, occupancy); + + // ── Cut 1: dr/cosOA ─────────────────────────────────────────────── + if (doFullRange) + fillFullRangeDeltaRCosOA<0>(obs.qinv, obs.drOverCosOA); + fRegistry.fill(HIST("Pair/same/hDeltaRCosOA"), obs.drOverCosOA); + if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA) continue; - ROOT::Math::XYZVector cp1(g1.vx(), g1.vy(), g1.vz()); - ROOT::Math::XYZVector cp2(g2.vx(), g2.vy(), g2.vz()); - float dr = std::sqrt(std::pow(g1.vx() - g2.vx(), 2) + std::pow(g1.vy() - g2.vy(), 2) + std::pow(g1.vz() - g2.vz(), 2)); - float opa = std::acos(std::clamp(static_cast(cp1.Dot(cp2) / (std::sqrt(cp1.Mag2()) * std::sqrt(cp2.Mag2()))), -1.f, 1.f)); - o2::math_utils::bringTo02Pi(opa); - if (opa > o2::constants::math::PI) - opa -= o2::constants::math::PI; - float cosOA = std::cos(opa / 2.f); - if (dr / cosOA < ggpaircuts.cfgMinDR_CosOA) + photonIdsAfterDRCosOA.insert(g1.globalIndex()); + photonIdsAfterDRCosOA.insert(g2.globalIndex()); + + // ── QA: After dr/cosOA cut ──────────────────────────────────────── + if (doQA) + fillPairQAStep<0, 1>(obs, centForQA, occupancy); + + // ── Cut 2: R/Z geometry ─────────────────────────────────────────── + if (!passRZCut(obs.deltaR, obs.deltaZ)) continue; - fRegistry.fill(HIST("Pair/same/hDeltaRCosOA"), dr / cosOA); - - // Kinematic variables for QA - ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); - ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); - float deta = g1.eta() - g2.eta(); - float dphi = RecoDecay::constrainAngle(g1.phi() - g2.phi(), -o2::constants::math::PI); - float pairEta = 0.5f * (g1.eta() + g2.eta()); - float cosTheta = std::fabs(computeCosTheta(v1, v2)); - float openingAngle = opa; - - // ── QA: Before ellipse cut ────────────────────────────────────── - fillPairQAStep<0, true>(deta, dphi, pairEta, cosTheta, openingAngle); - - // ── Ellipse cut ───────────────────────────────────────────────── - if (isInsideEllipse(deta, dphi)) + + photonIdsAfterRZ.insert(g1.globalIndex()); + photonIdsAfterRZ.insert(g2.globalIndex()); + + // ── QA: After R/Z cut ───────────────────────────────────────────── + if (doQA) + fillPairQAStep<0, 2>(obs, centForQA, occupancy); + + // ── Cut 3: Ellipse in (DeltaEta, DeltaPhi) ──────────────────────── + if (isInsideEllipse(obs.deta, obs.dphi)) continue; - // ── QA: After ellipse cut ─────────────────────────────────────── - fillPairQAStep<0, false>(deta, dphi, pairEta, cosTheta, openingAngle); + photonIdsAfterEllipse.insert(g1.globalIndex()); + photonIdsAfterEllipse.insert(g2.globalIndex()); + + // ── QA: After ellipse cut = final accepted pairs ────────────────── + if (doQA) + fillPairQAStep<0, 3>(obs, centForQA, occupancy); + + if (doFullRange) + fillFullRangeQA<0>(obs, centForQA, occupancy); - // ── Physics ───────────────────────────────────────────────── - fillPairHistogram<0>(collision, v1, v2, 1.f); + fillPairHistogram<0>(collision, obs.v1, obs.v2, 1.f); ndiphoton++; auto addToPool = [&](auto const& g) { - if (std::find(usedPhotonIdsPerCol.begin(), usedPhotonIdsPerCol.end(), - g.globalIndex()) == usedPhotonIdsPerCol.end()) { - EMPair gtmp(g.pt(), g.eta(), g.phi(), 0); + if (usedPhotonIdsPerCol.insert(g.globalIndex()).second) { + EMPair gtmp(g.pt(), g.eta(), g.phi(), 0.f); gtmp.setConversionPointXYZ(g.vx(), g.vy(), g.vz()); emh1->AddTrackToEventPool(keyDFCollision, gtmp); - usedPhotonIdsPerCol.emplace_back(g.globalIndex()); } }; addToPool(g1); addToPool(g2); + } // end same-event pair loop - // end same-event pair loop + if (qaflags.doSinglePhotonQa) { + for (const auto& g : photons1Coll) { + if (!cut1.template IsSelected(g)) + continue; + const int gid = g.globalIndex(); + if (photonIdsAfterDRCosOA.count(gid)) + fillSinglePhotonQAStep<1>(g); + if (photonIdsAfterRZ.count(gid)) + fillSinglePhotonQAStep<2>(g); + if (photonIdsAfterEllipse.count(gid)) + fillSinglePhotonQAStep<3>(g); + } } usedPhotonIdsPerCol.clear(); - usedPhotonIdsPerCol.shrink_to_fit(); - // ── Mixed-event loop ──────────────────────────────────────────────────── - if (!cfgDoMix || !(ndiphoton > 0)) + if (!cfgDoMix || ndiphoton == 0) continue; auto selectedPhotons = emh1->GetTracksPerCollision(keyDFCollision); auto poolIDs = emh1->GetCollisionIdsFromEventPool(keyBin); for (const auto& mixID : poolIDs) { + // skip same event if (mixID.second == collision.globalIndex() && mixID.first == ndf) continue; - uint64_t bcMix = mapMixedEventIdToGlobalBC[mixID]; - uint64_t diffBC = std::max(collision.globalBC(), bcMix) - std::min(collision.globalBC(), bcMix); + const uint64_t bcMix = mapMixedEventIdToGlobalBC[mixID]; + const uint64_t diffBC = std::max(collision.globalBC(), bcMix) - + std::min(collision.globalBC(), bcMix); fRegistry.fill(HIST("Pair/mix/hDiffBC"), diffBC); if (diffBC < ndiffBCMix) continue; @@ -706,40 +1467,53 @@ struct photonhbt { for (const auto& g1 : selectedPhotons) { for (const auto& g2 : poolPhotons) { - ROOT::Math::XYZVector cp1(g1.vx(), g1.vy(), g1.vz()); - ROOT::Math::XYZVector cp2(g2.vx(), g2.vy(), g2.vz()); - float dr = std::sqrt(std::pow(g1.vx() - g2.vx(), 2) + std::pow(g1.vy() - g2.vy(), 2) + std::pow(g1.vz() - g2.vz(), 2)); - float opa = std::acos(std::clamp(static_cast(cp1.Dot(cp2) / (std::sqrt(cp1.Mag2()) * std::sqrt(cp2.Mag2()))), -1.f, 1.f)); - o2::math_utils::bringTo02Pi(opa); - if (opa > o2::constants::math::PI) - opa -= o2::constants::math::PI; - float cosOA = std::cos(opa / 2.f); - if (dr / cosOA < ggpaircuts.cfgMinDR_CosOA) + auto obs = buildPairQAObservables(g1, g2); + if (!obs.valid) continue; - fRegistry.fill(HIST("Pair/mix/hDeltaRCosOA"), dr / cosOA); - ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); - ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); - float deta = g1.eta() - g2.eta(); - float dphi = RecoDecay::constrainAngle(g1.phi() - g2.phi(), -o2::constants::math::PI); - float pairEta = 0.5f * (g1.eta() + g2.eta()); - float cosTheta = std::fabs(computeCosTheta(v1, v2)); - float openingAngle = opa; + const bool doQA = passQinvQAGate(obs.qinv); + const bool doFullRange = passQinvFullRangeGate(obs.qinv); + + // ── QA: Before any pair cut ───────────────────────────────── + if (doQA) + fillPairQAStep<1, 0>(obs, centForQA, occupancy); + + // ── Cut 1: dr/cosOA ───────────────────────────────────────── + if (doFullRange) + fillFullRangeDeltaRCosOA<1>(obs.qinv, obs.drOverCosOA); + fRegistry.fill(HIST("Pair/mix/hDeltaRCosOA"), obs.drOverCosOA); + if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA) + continue; + + // ── QA: After dr/cosOA cut ────────────────────────────────── + if (doQA) + fillPairQAStep<1, 1>(obs, centForQA, occupancy); + + // ── Cut 2: R/Z geometry ───────────────────────────────────── + if (!passRZCut(obs.deltaR, obs.deltaZ)) + continue; - // QA Before cut — mix/QA/Before/ histograms. - fillPairQAStep<1, true>(deta, dphi, pairEta, cosTheta, openingAngle); + // ── QA: After R/Z cut ─────────────────────────────────────── + if (doQA) + fillPairQAStep<1, 2>(obs, centForQA, occupancy); - // Apply ellipse cut - if (isInsideEllipse(deta, dphi)) + // ── Cut 3: Ellipse ────────────────────────────────────────── + if (isInsideEllipse(obs.deta, obs.dphi)) continue; - // QA After cut - fillPairQAStep<1, false>(deta, dphi, pairEta, cosTheta, openingAngle); + // ── QA: After ellipse cut ─────────────────────────────────── + if (doQA) + fillPairQAStep<1, 3>(obs, centForQA, occupancy); - fillPairHistogram<1>(collision, v1, v2, 1.f); + // ── Full-range fills ──────────────────────────────────────── + if (doFullRange) + fillFullRangeQA<1>(obs, centForQA, occupancy); + + // ── Fill CF histogram — always ────────────────────────────── + fillPairHistogram<1>(collision, obs.v1, obs.v2, 1.f); } } - } + } // end mixed-event loop if (ndiphoton > 0) { emh1->AddCollisionIdAtLast(keyBin, keyDFCollision); @@ -749,36 +1523,311 @@ struct photonhbt { } // end collision loop } - using MyEMH = o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, EMPair>; + using MyEMH = o2::aod::pwgem::dilepton::utils::EventMixingHandler< + std::tuple, + std::pair, + EMPair>; + MyEMH* emh1 = nullptr; MyEMH* emh2 = nullptr; - std::vector usedPhotonIdsPerCol; + + std::unordered_set usedPhotonIdsPerCol; std::map, uint64_t> mapMixedEventIdToGlobalBC; SliceCache cache; Preslice perCollisionPCM = aod::v0photonkf::pmeventId; - Filter collisionFilterCentrality = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || - (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || - (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); - Filter collisionFilterOccupancyTrack = eventcuts.cfgTrackOccupancyMin <= o2::aod::evsel::trackOccupancyInTimeRange && - o2::aod::evsel::trackOccupancyInTimeRange < eventcuts.cfgTrackOccupancyMax; - Filter collisionFilterOccupancyFT0c = eventcuts.cfgFT0COccupancyMin <= o2::aod::evsel::ft0cOccupancyInTimeRange && - o2::aod::evsel::ft0cOccupancyInTimeRange < eventcuts.cfgFT0COccupancyMax; + Filter collisionFilterCentrality = + (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || + (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || + (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); + Filter collisionFilterOccupancyTrack = + eventcuts.cfgTrackOccupancyMin <= o2::aod::evsel::trackOccupancyInTimeRange && + o2::aod::evsel::trackOccupancyInTimeRange < eventcuts.cfgTrackOccupancyMax; + Filter collisionFilterOccupancyFT0c = + eventcuts.cfgFT0COccupancyMin <= o2::aod::evsel::ft0cOccupancyInTimeRange && + o2::aod::evsel::ft0cOccupancyInTimeRange < eventcuts.cfgFT0COccupancyMax; + using FilteredMyCollisions = soa::Filtered; int ndf = 0; - void processAnalysis(FilteredMyCollisions const& collisions, MyV0Photons const& v0photons, aod::V0Legs const& v0legs) + + void processAnalysis(FilteredMyCollisions const& collisions, + MyV0Photons const& v0photons, + aod::V0Legs const& v0legs) { - runPairing(collisions, v0photons, v0photons, v0legs, v0legs, - perCollisionPCM, perCollisionPCM, fV0PhotonCut, fV0PhotonCut); + runPairing(collisions, + v0photons, v0photons, + v0legs, v0legs, + perCollisionPCM, perCollisionPCM, + fV0PhotonCut, fV0PhotonCut); ndf++; } + PROCESS_SWITCH(photonhbt, processAnalysis, "pairing for analysis", true); + + template + void runPairingMC(TCollisions const& collisions, + TPhotons const& photons, + TLegs const& /*legs*/, + TMCParticles const& mcParticles, + TPreslice const& perCollision, + TCut const& cut) + { + for (const auto& collision : collisions) { + initCCDB(collision); + int ndiphoton = 0; + + const float cent[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; + if (cent[cfgCentEstimator] < cfgCentMin || cfgCentMax < cent[cfgCentEstimator]) + continue; + + const std::array epArr = { + collision.ep2ft0m(), collision.ep2ft0a(), collision.ep2ft0c(), + collision.ep2fv0a(), collision.ep2btot(), collision.ep2bpos(), collision.ep2bneg()}; + const float ep2 = epArr[cfgEP2EstimatorForMix]; + + fRegistry.fill(HIST("Event/before/hEP2_CentFT0C_forMix"), collision.centFT0C(), ep2); + o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(&fRegistry, collision, 1.f); + if (!fEMEventCut.IsSelected(collision)) + continue; + + o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(&fRegistry, collision, 1.f); + fRegistry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted + fRegistry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted + fRegistry.fill(HIST("Event/after/hEP2_CentFT0C_forMix"), collision.centFT0C(), ep2); + + const float occupancy = (cfgOccupancyEstimator == 1) + ? static_cast(collision.trackOccupancyInTimeRange()) + : collision.ft0cOccupancyInTimeRange(); + const float centForQA = cent[cfgCentEstimator]; + + const int zbin = binOf(ztxBinEdges, collision.posZ()); + const int centbin = binOf(centBinEdges, centForQA); + const int epbin = binOf(epBinEgdes, ep2); + const int occbin = binOf(occBinEdges, occupancy); + + auto keyBin = std::make_tuple(zbin, centbin, epbin, occbin); + auto keyDFCollision = std::make_pair(ndf, collision.globalIndex()); + + auto photonsColl = photons.sliceBy(perCollision, collision.globalIndex()); + + if (qaflags.doSinglePhotonQa) { + for (const auto& g : photonsColl) { + if (!cut.template IsSelected(g)) + continue; + fillSinglePhotonQAStep<0>(g); + } + } + + std::unordered_set photonIdsAfterDRCosOA; + std::unordered_set photonIdsAfterRZ; + std::unordered_set photonIdsAfterEllipse; + + // ── Same-event pair loop ────────────────────────────────────────────── + for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photonsColl, photonsColl))) { + if (!cut.template IsSelected(g1) || + !cut.template IsSelected(g2)) + continue; + + const auto pos1 = g1.template posTrack_as(); + const auto ele1 = g1.template negTrack_as(); + const auto pos2 = g2.template posTrack_as(); + const auto ele2 = g2.template negTrack_as(); + if (pos1.trackId() == pos2.trackId() || + pos1.trackId() == ele2.trackId() || + ele1.trackId() == pos2.trackId() || + ele1.trackId() == ele2.trackId()) + continue; + + // ── MC truth classification ─────────────────────────────────────── + const auto mc1 = buildPhotonMCInfo(g1, mcParticles); + const auto mc2 = buildPhotonMCInfo(g2, mcParticles); + auto truthType = classifyPairTruth(mc1, mc2); + if (truthType == PairTruthType::TrueTrueDistinct && + isPi0DaughterPair(mc1, mc2, mcParticles)) + truthType = PairTruthType::Pi0Daughters; + + auto obs = buildPairQAObservables(g1, g2); + if (!obs.valid) + continue; + + const bool doQA = passQinvQAGate(obs.qinv); + const bool doFullRange = passQinvFullRangeGate(obs.qinv); + + // ── Pair QA: Before ─────────────────────────────────────────────── + if (doQA) + fillPairQAStep<0, 0>(obs, centForQA, occupancy); + + // ── Cut 1: dr/cosOA ─────────────────────────────────────────────── + if (doFullRange) + fillFullRangeDeltaRCosOA<0>(obs.qinv, obs.drOverCosOA); + fRegistry.fill(HIST("Pair/same/hDeltaRCosOA"), obs.drOverCosOA); + if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA) + continue; + + photonIdsAfterDRCosOA.insert(g1.globalIndex()); + photonIdsAfterDRCosOA.insert(g2.globalIndex()); + if (doQA) + fillPairQAStep<0, 1>(obs, centForQA, occupancy); + + // ── Cut 2: R/Z geometry ─────────────────────────────────────────── + if (!passRZCut(obs.deltaR, obs.deltaZ)) + continue; + + photonIdsAfterRZ.insert(g1.globalIndex()); + photonIdsAfterRZ.insert(g2.globalIndex()); + if (doQA) + fillPairQAStep<0, 2>(obs, centForQA, occupancy); + + // ── Cut 3: Ellipse ──────────────────────────────────────────────── + if (isInsideEllipse(obs.deta, obs.dphi)) + continue; + + photonIdsAfterEllipse.insert(g1.globalIndex()); + photonIdsAfterEllipse.insert(g2.globalIndex()); + if (doQA) + fillPairQAStep<0, 3>(obs, centForQA, occupancy); + + // ── Full-range fills ────────────────────────────────────────────── + if (doFullRange) + fillFullRangeQA<0>(obs, centForQA, occupancy); + + // ── Fill inclusive CF — always ──────────────────────────────────── + fillPairHistogram<0>(collision, obs.v1, obs.v2, 1.f); + ndiphoton++; + + if (doQA) + fillMCPairQA(truthType, obs); + if (doFullRange) + fillMCPairQAFullRange(truthType, obs); + switch (truthType) { + case PairTruthType::TrueTrueDistinct: + fillPairHistogramMC<0, PairTruthType::TrueTrueDistinct>(collision, obs.v1, obs.v2); + break; + case PairTruthType::TrueTrueSamePhoton: + fillPairHistogramMC<0, PairTruthType::TrueTrueSamePhoton>(collision, obs.v1, obs.v2); + break; + case PairTruthType::SharedMcLeg: + fillPairHistogramMC<0, PairTruthType::SharedMcLeg>(collision, obs.v1, obs.v2); + break; + case PairTruthType::TrueFake: + fillPairHistogramMC<0, PairTruthType::TrueFake>(collision, obs.v1, obs.v2); + break; + case PairTruthType::FakeFake: + fillPairHistogramMC<0, PairTruthType::FakeFake>(collision, obs.v1, obs.v2); + break; + case PairTruthType::Pi0Daughters: + fillPairHistogramMC<0, PairTruthType::Pi0Daughters>(collision, obs.v1, obs.v2); + break; + default: + break; + } + + auto addToPool = [&](auto const& g) { + if (usedPhotonIdsPerCol.insert(g.globalIndex()).second) { + EMPair gtmp(g.pt(), g.eta(), g.phi(), 0.f); + gtmp.setConversionPointXYZ(g.vx(), g.vy(), g.vz()); + emh1->AddTrackToEventPool(keyDFCollision, gtmp); + } + }; + addToPool(g1); + addToPool(g2); + } // end same-event pair loop + + if (qaflags.doSinglePhotonQa) { + for (const auto& g : photonsColl) { + if (!cut.template IsSelected(g)) + continue; + const int gid = g.globalIndex(); + if (photonIdsAfterDRCosOA.count(gid)) + fillSinglePhotonQAStep<1>(g); + if (photonIdsAfterRZ.count(gid)) + fillSinglePhotonQAStep<2>(g); + if (photonIdsAfterEllipse.count(gid)) + fillSinglePhotonQAStep<3>(g); + } + } + + usedPhotonIdsPerCol.clear(); + + if (!cfgDoMix || ndiphoton == 0) + continue; + + auto selectedPhotons = emh1->GetTracksPerCollision(keyDFCollision); + auto poolIDs = emh1->GetCollisionIdsFromEventPool(keyBin); + + for (const auto& mixID : poolIDs) { + if (mixID.second == collision.globalIndex() && mixID.first == ndf) + continue; + const uint64_t bcMix = mapMixedEventIdToGlobalBC[mixID]; + const uint64_t diffBC = std::max(collision.globalBC(), bcMix) - + std::min(collision.globalBC(), bcMix); + fRegistry.fill(HIST("Pair/mix/hDiffBC"), diffBC); + if (diffBC < ndiffBCMix) + continue; + + auto poolPhotons = emh1->GetTracksPerCollision(mixID); + for (const auto& g1 : selectedPhotons) { + for (const auto& g2 : poolPhotons) { + auto obs = buildPairQAObservables(g1, g2); + if (!obs.valid) + continue; + const bool doQA = passQinvQAGate(obs.qinv); + const bool doFullRange = passQinvFullRangeGate(obs.qinv); + if (doQA) + fillPairQAStep<1, 0>(obs, centForQA, occupancy); + if (doFullRange) + fillFullRangeDeltaRCosOA<1>(obs.qinv, obs.drOverCosOA); + fRegistry.fill(HIST("Pair/mix/hDeltaRCosOA"), obs.drOverCosOA); + if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA) + continue; + if (doQA) + fillPairQAStep<1, 1>(obs, centForQA, occupancy); + if (!passRZCut(obs.deltaR, obs.deltaZ)) + continue; + if (doQA) + fillPairQAStep<1, 2>(obs, centForQA, occupancy); + if (isInsideEllipse(obs.deta, obs.dphi)) + continue; + if (doQA) + fillPairQAStep<1, 3>(obs, centForQA, occupancy); + if (doFullRange) + fillFullRangeQA<1>(obs, centForQA, occupancy); + fillPairHistogram<1>(collision, obs.v1, obs.v2, 1.f); + } + } + } + + if (ndiphoton > 0) { + emh1->AddCollisionIdAtLast(keyBin, keyDFCollision); + emh2->AddCollisionIdAtLast(keyBin, keyDFCollision); + mapMixedEventIdToGlobalBC[keyDFCollision] = collision.globalBC(); + } + } // end collision loop + } + + void processMC(FilteredMyCollisions const& collisions, + MyV0Photons const& v0photons, + MyMCV0Legs const& v0legs, + aod::EMMCParticles const& mcParticles, + aod::EMMCEvents const& /*mcEvents*/) + { + runPairingMC(collisions, v0photons, v0legs, mcParticles, + perCollisionPCM, fV0PhotonCut); + ndf++; + } + + PROCESS_SWITCH(photonhbt, processMC, "pairing with MC truth classification", false); }; +// ============================================================================ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{ - adaptAnalysisTask(cfgc, TaskName{"photonhbt"})}; + return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"photonhbt"})}; }