From 00f5e3f9020fbf3dfbc956ea1a9fe5cd3df5a7f0 Mon Sep 17 00:00:00 2001 From: Stefanie Mrozinski Date: Mon, 16 Mar 2026 14:42:26 +0100 Subject: [PATCH 1/3] Add new QC histograms for photon pairs --- PWGEM/PhotonMeson/Core/PhotonHBT.h | 854 +++++++++++++++++------------ 1 file changed, 494 insertions(+), 360 deletions(-) diff --git a/PWGEM/PhotonMeson/Core/PhotonHBT.h b/PWGEM/PhotonMeson/Core/PhotonHBT.h index fcd9d9d9ce7..364189d0fe7 100644 --- a/PWGEM/PhotonMeson/Core/PhotonHBT.h +++ b/PWGEM/PhotonMeson/Core/PhotonHBT.h @@ -64,32 +64,74 @@ enum class ggHBTPairType : int { }; } // namespace o2::aod::pwgem::photon::core::photonhbt -using namespace o2; -using namespace o2::aod; -using namespace o2::framework; -using namespace o2::framework::expressions; -using namespace o2::soa; -using namespace o2::aod::pwgem::dilepton::utils; -using namespace o2::aod::pwgem::photon::core::photonhbt; +using namespace o2; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) +using namespace o2::aod; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) +using namespace o2::framework; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) +using namespace o2::framework::expressions; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) +using namespace o2::soa; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) +using namespace o2::aod::pwgem::dilepton::utils; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) +using namespace o2::aod::pwgem::photon::core::photonhbt; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) using MyCollisions = soa::Join; using MyCollision = MyCollisions::iterator; -using MyV0Photons = soa::Join; +using MyV0Photons = soa::Join; using MyV0Photon = MyV0Photons::iterator; template struct PhotonHBT { - // Configurables - // Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - // Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; - // Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - // Configurable skipGRPOquery{"skipGRPOquery", true, "skip grpo query"}; - // Configurable d_bz_input{"d_bz_input", -999, "bz field in kG, -999 is automatic"}; + + // Single-photon: + // 0 = Inclusive (all photons) + // 1 = ITSTPC_ITSTPC — both legs have ITS+TPC + // 2 = ITSTPC_TPCOnly — one ITS+TPC, one TPC-only + // 3 = TPCOnly_TPCOnly — both legs TPC-only + // + // Pair combo. + // 0 = Inclusive (all recognised pairs) + // 1 = ITSTPC_ITSTPC × ITSTPC_ITSTPC + // 2 = ITSTPC_ITSTPC × ITSTPC_TPCOnly + // 3 = ITSTPC_ITSTPC × TPCOnly_TPCOnly + // 4 = ITSTPC_TPCOnly × ITSTPC_TPCOnly + // 5 = ITSTPC_TPCOnly × TPCOnly_TPCOnly + // 6 = TPCOnly_TPCOnly × TPCOnly_TPCOnly + + template + static inline int classifyV0ComboIdx(TGamma const& g) + { + auto pos = g.template posTrack_as(); + auto neg = g.template negTrack_as(); + const bool posII = pos.hasITS() && pos.hasTPC(); + const bool posTPC = !pos.hasITS() && pos.hasTPC(); + const bool negII = neg.hasITS() && neg.hasTPC(); + const bool negTPC = !neg.hasITS() && neg.hasTPC(); + if (posII && negII) + return 1; + if ((posII && negTPC) || (posTPC && negII)) + return 2; + if (posTPC && negTPC) + return 3; + return 0; + } + + static inline int pairComboBin(int c1, int c2) + { + if (c1 <= 0 || c2 <= 0) + return 0; + if (c1 > c2) + std::swap(c1, c2); + static constexpr int kTable[4][4] = { + {0, 0, 0, 0}, + {0, 1, 2, 3}, + {0, 2, 4, 5}, + {0, 3, 5, 6}}; + return kTable[c1][c2]; + } Configurable cfgDo3D{"cfgDo3D", false, "enable 3D analysis"}; - Configurable cfgEP2Estimator_for_Mix{"cfgEP2Estimator_for_Mix", 3, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5"}; + Configurable cfgEp2EstimatorForMix{"cfgEp2EstimatorForMix", 3, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5"}; + Configurable cfgEp3EstimatorForMix{"cfgEp3EstimatorForMix", 3, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5 (for psi3 mixing)"}; Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; Configurable cfgOccupancyEstimator{"cfgOccupancyEstimator", 0, "FT0C:0, Track:1"}; Configurable cfgCentMin{"cfgCentMin", -1, "min. centrality"}; @@ -97,183 +139,194 @@ struct PhotonHBT { Configurable maxY{"maxY", 0.8, "maximum rapidity for reconstructed particles"}; Configurable cfgDoMix{"cfgDoMix", true, "flag for event mixing"}; Configurable ndepth{"ndepth", 100, "depth for event mixing"}; - Configurable ndiff_bc_mix{"ndiff_bc_mix", 594, "difference in global BC required in mixed events"}; - ConfigurableAxis ConfVtxBins{"ConfVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; - ConfigurableAxis ConfCentBins{"ConfCentBins", {VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.f, 999.f}, "Mixing bins - centrality"}; - ConfigurableAxis ConfEPBins{"ConfEPBins", {16, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}, "Mixing bins - event plane angle"}; - ConfigurableAxis ConfOccupancyBins{"ConfOccupancyBins", {VARIABLE_WIDTH, -1, 1e+10}, "Mixing bins - occupancy"}; - Configurable cfgUseLCMS{"cfgUseLCMS", true, "measure relative momentum in LCMS for 1D"}; // always in LCMS for 3D + Configurable ndiffBcMix{"ndiffBcMix", 594, "difference in global BC required in mixed events"}; + Configurable cfgUseLcms{"cfgUseLcms", true, "measure relative momentum in LCMS for 1D"}; + + ConfigurableAxis confVtxBins{"confVtxBins", {VARIABLE_WIDTH, -10.f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; + ConfigurableAxis confCentBins{"confCentBins", {VARIABLE_WIDTH, 0.f, 5.f, 10.f, 20.f, 30.f, 40.f, 50.f, 60.f, 70.f, 80.f, 90.f, 100.f, 999.f}, "Mixing bins - centrality"}; + ConfigurableAxis confEpBins{"confEpBins", {8, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}, "Mixing bins - psi2 event plane"}; + + ConfigurableAxis confEp3Bins{"confEp3Bins", {1, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}, "Mixing bins - psi3 event plane (set to 1 bin to disable)"}; + ConfigurableAxis confOccupancyBins{"confOccupancyBins", {VARIABLE_WIDTH, -1, 1e+10}, "Mixing bins - occupancy"}; - ConfigurableAxis ConfQBins{"ConfQBins", {60, 0, +0.3f}, "q bins for output histograms"}; - ConfigurableAxis ConfKtBins{"ConfKtBins", {VARIABLE_WIDTH, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0}, "kT bins for output histograms"}; + ConfigurableAxis confQBins{"confQBins", {60, 0, +0.3f}, "q bins for output histograms"}; // o2-linter: disable=name/configurable (Q is a physics symbol for momentum transfer; cannot be lowercased without losing meaning) + ConfigurableAxis confKtBins{"confKtBins", {VARIABLE_WIDTH, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0}, "kT bins for output histograms"}; + + // QA axis configurables + 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 confDeltaEtaBins{"confDeltaEtaBins", {100, -0.5f, +0.5f}, "Delta-eta bins"}; + ConfigurableAxis confDeltaPhiBins{"confDeltaPhiBins", {100, -0.5f, +0.5f}, "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)"}; + + ConfigurableAxis confRxyBins{"confRxyBins", {90, 0.f, 90.f}, "Rxy bins (cm)"}; + + ConfigurableAxis confPsiPairBins{"confPsiPairBins", {100, -0.2f, 0.2f}, "psi_pair bins (rad)"}; + + // =========================================================================== + // QA flags + // =========================================================================== + + Configurable doPairQa{"doPairQa", true, "fill pair QA histograms (Before/After ellipse cut, with combo axis)"}; + Configurable doSinglePhotonQa{"doSinglePhotonQa", true, "fill single-photon QA histograms (pT, eta, phi with V0 combo axis)"}; + + // =========================================================================== + // Pair cuts + // =========================================================================== + + struct : ConfigurableGroup { + std::string prefix = "ggpaircut_group"; + Configurable cfgMinDrCosOa{"cfgMinDrCosOa", -1, "min. dr/cosOA for kPCMPCM"}; + Configurable cfgApplyEllipseCut{"cfgApplyEllipseCut", 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"}; + } ggpaircuts; EMPhotonEventCut fEMEventCut; struct : ConfigurableGroup { std::string prefix = "eventcut_group"; Configurable cfgZvtxMin{"cfgZvtxMin", -10.f, "min. Zvtx"}; Configurable cfgZvtxMax{"cfgZvtxMax", +10.f, "max. Zvtx"}; - Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8 in event cut"}; - Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND in event cut"}; - Configurable cfgRequireNoTFB{"cfgRequireNoTFB", true, "require No time frame border in event cut"}; - Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", true, "require no ITS readout frame border in event cut"}; - Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup in event cut"}; - Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC in event cut"}; // ITS-TPC matched track contributes PV. - Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx between FT0 vs. PV in event cut"}; + Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8"}; + Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgRequireNoTFB{"cfgRequireNoTFB", true, "require no TF border"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", true, "require no ITS ROF border"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup"}; + Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx FT0 vs PV"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) Configurable cfgTrackOccupancyMin{"cfgTrackOccupancyMin", -2, "min. track occupancy"}; Configurable cfgTrackOccupancyMax{"cfgTrackOccupancyMax", 1000000000, "max. track occupancy"}; - Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -2, "min. FT0C occupancy"}; - Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000, "max. FT0C occupancy"}; + Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -2.f, "min. FT0C occupancy"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000.f, "max. FT0C occupancy"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) Configurable cfgRequireNoCollInTimeRangeStandard{"cfgRequireNoCollInTimeRangeStandard", false, "require no collision in time range standard"}; Configurable cfgRequireNoCollInTimeRangeStrict{"cfgRequireNoCollInTimeRangeStrict", false, "require no collision in time range strict"}; - Configurable cfgRequireNoCollInITSROFStandard{"cfgRequireNoCollInITSROFStandard", false, "require no collision in time range standard"}; - Configurable cfgRequireNoCollInITSROFStrict{"cfgRequireNoCollInITSROFStrict", false, "require no collision in time range strict"}; - Configurable cfgRequireNoHighMultCollInPrevRof{"cfgRequireNoHighMultCollInPrevRof", false, "require no HM collision in previous ITS ROF"}; - Configurable cfgRequireGoodITSLayer3{"cfgRequireGoodITSLayer3", false, "number of inactive chips on ITS layer 3 are below threshold "}; - Configurable cfgRequireGoodITSLayer0123{"cfgRequireGoodITSLayer0123", false, "number of inactive chips on ITS layers 0-3 are below threshold "}; - Configurable cfgRequireGoodITSLayersAll{"cfgRequireGoodITSLayersAll", false, "number of inactive chips on all ITS layers are below threshold "}; + Configurable cfgRequireNoCollInITSROFStandard{"cfgRequireNoCollInITSROFStandard", false, "require no collision in time range standard"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgRequireNoCollInITSROFStrict{"cfgRequireNoCollInITSROFStrict", false, "require no collision in time range strict"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgRequireNoHighMultCollInPrevRof{"cfgRequireNoHighMultCollInPrevRof", false, "require no HM collision in previous ITS ROF"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgRequireGoodITSLayer3{"cfgRequireGoodITSLayer3", false, "ITS layer 3 chips OK"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgRequireGoodITSLayer0123{"cfgRequireGoodITSLayer0123", false, "ITS layers 0-3 chips OK"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgRequireGoodITSLayersAll{"cfgRequireGoodITSLayersAll", false, "all ITS layers chips OK"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) } eventcuts; V0PhotonCut fV0PhotonCut; struct : ConfigurableGroup { std::string prefix = "pcmcut_group"; - Configurable cfg_require_v0_with_itstpc{"cfg_require_v0_with_itstpc", false, "flag to select V0s with ITS-TPC matched tracks"}; - Configurable cfg_require_v0_with_itsonly{"cfg_require_v0_with_itsonly", false, "flag to select V0s with ITSonly tracks"}; - Configurable cfg_require_v0_with_tpconly{"cfg_require_v0_with_tpconly", false, "flag to select V0s with TPConly tracks"}; - Configurable cfg_min_pt_v0{"cfg_min_pt_v0", 0.1, "min pT for v0 photons at PV"}; - Configurable cfg_max_eta_v0{"cfg_max_eta_v0", 0.8, "max eta for v0 photons at PV"}; - Configurable cfg_min_v0radius{"cfg_min_v0radius", 16.0, "min v0 radius"}; - Configurable cfg_max_v0radius{"cfg_max_v0radius", 90.0, "max v0 radius"}; - Configurable cfg_max_alpha_ap{"cfg_max_alpha_ap", 0.95, "max alpha for AP cut"}; - Configurable cfg_max_qt_ap{"cfg_max_qt_ap", 0.01, "max qT for AP cut"}; - Configurable cfg_min_cospa{"cfg_min_cospa", 0.997, "min V0 CosPA"}; - Configurable cfg_max_pca{"cfg_max_pca", 3.0, "max distance btween 2 legs"}; - Configurable cfg_max_chi2kf{"cfg_max_chi2kf", 1e+10, "max chi2/ndf with KF"}; - Configurable cfg_reject_v0_on_itsib{"cfg_reject_v0_on_itsib", true, "flag to reject V0s on ITSib"}; - - Configurable cfg_disable_itsonly_track{"cfg_disable_itsonly_track", false, "flag to disable ITSonly tracks"}; - Configurable cfg_disable_tpconly_track{"cfg_disable_tpconly_track", false, "flag to disable TPConly tracks"}; - Configurable cfg_min_ncluster_tpc{"cfg_min_ncluster_tpc", 0, "min ncluster tpc"}; - Configurable cfg_min_ncrossedrows{"cfg_min_ncrossedrows", 40, "min ncrossed rows"}; - Configurable cfg_max_frac_shared_clusters_tpc{"cfg_max_frac_shared_clusters_tpc", 999.f, "max fraction of shared clusters in TPC"}; - Configurable cfg_max_chi2tpc{"cfg_max_chi2tpc", 4.0, "max chi2/NclsTPC"}; - Configurable cfg_max_chi2its{"cfg_max_chi2its", 36.0, "max chi2/NclsITS"}; - Configurable cfg_min_TPCNsigmaEl{"cfg_min_TPCNsigmaEl", -3.0, "min. TPC n sigma for electron"}; - Configurable cfg_max_TPCNsigmaEl{"cfg_max_TPCNsigmaEl", +3.0, "max. TPC n sigma for electron"}; + Configurable cfgRequireV0WithItstpc{"cfgRequireV0WithItstpc", false, "flag to select V0s with ITS-TPC matched tracks"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgRequireV0WithItsonly{"cfgRequireV0WithItsonly", false, "flag to select V0s with ITSonly tracks"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgRequireV0WithTpconly{"cfgRequireV0WithTpconly", false, "flag to select V0s with TPConly tracks"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMinPtV0{"cfgMinPtV0", 0.1, "min pT for v0 photons at PV"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMaxEtaV0{"cfgMaxEtaV0", 0.8, "max eta for v0 photons at PV"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMinV0Radius{"cfgMinV0Radius", 16.0, "min v0 radius"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMaxV0Radius{"cfgMaxV0Radius", 90.0, "max v0 radius"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMaxAlphaAp{"cfgMaxAlphaAp", 0.95, "max alpha for AP cut"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMaxQtAp{"cfgMaxQtAp", 0.01, "max qT for AP cut"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMinCospa{"cfgMinCospa", 0.997, "min V0 CosPA"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMaxPca{"cfgMaxPca", 3.0, "max distance btween 2 legs"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMaxChi2Kf{"cfgMaxChi2Kf", 1e+10, "max chi2/ndf with KF"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgRejectV0OnItsib{"cfgRejectV0OnItsib", true, "flag to reject V0s on ITSib"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgDisableItsonlyTrack{"cfgDisableItsonlyTrack", false, "flag to disable ITSonly tracks"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgDisableTpconlyTrack{"cfgDisableTpconlyTrack", false, "flag to disable TPConly tracks"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMinNclusterTpc{"cfgMinNclusterTpc", 0, "min ncluster tpc"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMinNcrossedrows{"cfgMinNcrossedrows", 40, "min ncrossed rows"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMaxFracSharedClustersTpc{"cfgMaxFracSharedClustersTpc", 999.f, "max fraction of shared clusters in TPC"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMaxChi2Tpc{"cfgMaxChi2Tpc", 4.0, "max chi2/NclsTPC"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMaxChi2Its{"cfgMaxChi2Its", 36.0, "max chi2/NclsITS"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMinTpcNsigmaEl{"cfgMinTpcNsigmaEl", -3.0, "min. TPC n sigma for electron"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) + Configurable cfgMaxTpcNsigmaEl{"cfgMaxTpcNsigmaEl", +3.0, "max. TPC n sigma for electron"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) } pcmcuts; - struct : ConfigurableGroup { - std::string prefix = "ggpaircut_group"; - Configurable cfgMinDR_CosOA{"cfgMinDR_CosOA", -1, "min. dr/cosOA for kPCMPCM"}; - } ggpaircuts; - ~PhotonHBT() { delete emh1; emh1 = 0x0; delete emh2; emh2 = 0x0; - map_mixed_eventId_to_globalBC.clear(); - used_photonIds_per_col.clear(); used_photonIds_per_col.shrink_to_fit(); } HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; - // static constexpr std::string_view event_types[2] = {"before", "after"}; static constexpr std::string_view event_pair_types[2] = {"same/", "mix/"}; std::mt19937 engine; std::uniform_int_distribution dist01; - - // o2::ccdb::CcdbApi ccdbApi; - // Service ccdb; int mRunNumber; - // float d_bz; std::vector zvtx_bin_edges; std::vector cent_bin_edges; - std::vector ep_bin_edges; + std::vector ep2_bin_edges; + std::vector ep3_bin_edges; std::vector occ_bin_edges; + static constexpr float MinFloatTol = 1e-9f; + static constexpr float CosineMin = -1.f; + static constexpr float CosineMax = 1.f; + static constexpr float HalfAngle = 2.f; + static constexpr float PairEtaWeight = 0.5f; + + inline bool isInsideEllipse(float deta, float dphi) const + { + if (!ggpaircuts.cfgApplyEllipseCut.value) + return false; + const float sE = ggpaircuts.cfgEllipseSigEta.value; + const float sP = ggpaircuts.cfgEllipseSigPhi.value; + if (sE < MinFloatTol || sP < MinFloatTol) + return false; + return (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP) < ggpaircuts.cfgEllipseR2.value; + } + + static inline float normDphi(float dphi) + { + while (dphi > o2::constants::math::PI) + dphi -= o2::constants::math::TwoPI; // o2-linter: disable=two-pi-add-subtract (PWGEM has no PWGHF dependency; manual wrap is established pattern) + while (dphi < -o2::constants::math::PI) + dphi += o2::constants::math::TwoPI; // o2-linter: disable=two-pi-add-subtract (PWGEM has no PWGHF dependency; manual wrap is established pattern) + return dphi; + } + + static inline float computeCosTheta(const ROOT::Math::PtEtaPhiMVector& v1, + const ROOT::Math::PtEtaPhiMVector& v2) + { + ROOT::Math::PxPyPzEVector p1(v1), p2(v2); + ROOT::Math::PxPyPzEVector pair = p1 + p2; + ROOT::Math::Boost boost(-pair.BoostToCM()); + ROOT::Math::PxPyPzEVector p1cm = boost(p1); + ROOT::Math::XYZVector pairDir(pair.Px(), pair.Py(), pair.Pz()); + ROOT::Math::XYZVector p1cmDir(p1cm.Px(), p1cm.Py(), p1cm.Pz()); + if (pairDir.R() < MinFloatTol || p1cmDir.R() < MinFloatTol) + return -1.f; + return static_cast(pairDir.Unit().Dot(p1cmDir.Unit())); + } + void init(InitContext& /*context*/) { mRunNumber = 0; - // d_bz = 0; - - // ccdb->setURL(ccdburl); - // ccdb->setCaching(true); - // ccdb->setLocalObjectValidityChecking(); - // ccdb->setFatalWhenNull(false); - - if (ConfVtxBins.value[0] == VARIABLE_WIDTH) { - zvtx_bin_edges = std::vector(ConfVtxBins.value.begin(), ConfVtxBins.value.end()); - zvtx_bin_edges.erase(zvtx_bin_edges.begin()); - for (const auto& edge : zvtx_bin_edges) { - LOGF(info, "VARIABLE_WIDTH: zvtx_bin_edges = %f", edge); - } - } else { - int nbins = static_cast(ConfVtxBins.value[0]); - float xmin = static_cast(ConfVtxBins.value[1]); - float xmax = static_cast(ConfVtxBins.value[2]); - zvtx_bin_edges.resize(nbins + 1); - for (int i = 0; i < nbins + 1; i++) { - zvtx_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; - LOGF(info, "FIXED_WIDTH: zvtx_bin_edges[%d] = %f", i, zvtx_bin_edges[i]); - } - } - if (ConfCentBins.value[0] == VARIABLE_WIDTH) { - cent_bin_edges = std::vector(ConfCentBins.value.begin(), ConfCentBins.value.end()); - cent_bin_edges.erase(cent_bin_edges.begin()); - for (const auto& edge : cent_bin_edges) { - LOGF(info, "VARIABLE_WIDTH: cent_bin_edges = %f", edge); - } - } else { - int nbins = static_cast(ConfCentBins.value[0]); - float xmin = static_cast(ConfCentBins.value[1]); - float xmax = static_cast(ConfCentBins.value[2]); - cent_bin_edges.resize(nbins + 1); - for (int i = 0; i < nbins + 1; i++) { - cent_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; - LOGF(info, "FIXED_WIDTH: cent_bin_edges[%d] = %f", i, cent_bin_edges[i]); + auto parseBinsVerbatim = [](const ConfigurableAxis& cfg, std::vector& edges) { + if (cfg.value[0] == VARIABLE_WIDTH) { + edges = std::vector(cfg.value.begin(), cfg.value.end()); + edges.erase(edges.begin()); + } else { + int nbins = static_cast(cfg.value[0]); + float xmin = static_cast(cfg.value[1]); + float xmax = static_cast(cfg.value[2]); + edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) + edges[i] = (xmax - xmin) / nbins * i + xmin; } - } - - if (ConfEPBins.value[0] == VARIABLE_WIDTH) { - ep_bin_edges = std::vector(ConfEPBins.value.begin(), ConfEPBins.value.end()); - ep_bin_edges.erase(ep_bin_edges.begin()); - for (const auto& edge : ep_bin_edges) { - LOGF(info, "VARIABLE_WIDTH: ep_bin_edges = %f", edge); - } - } else { - int nbins = static_cast(ConfEPBins.value[0]); - float xmin = static_cast(ConfEPBins.value[1]); - float xmax = static_cast(ConfEPBins.value[2]); - ep_bin_edges.resize(nbins + 1); - for (int i = 0; i < nbins + 1; i++) { - ep_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; - LOGF(info, "FIXED_WIDTH: ep_bin_edges[%d] = %f", i, ep_bin_edges[i]); - } - } + }; + parseBinsVerbatim(confVtxBins, zvtx_bin_edges); + parseBinsVerbatim(confCentBins, cent_bin_edges); + parseBinsVerbatim(confEpBins, ep2_bin_edges); + parseBinsVerbatim(confEp3Bins, ep3_bin_edges); LOGF(info, "cfgOccupancyEstimator = %d", cfgOccupancyEstimator.value); - if (ConfOccupancyBins.value[0] == VARIABLE_WIDTH) { - occ_bin_edges = std::vector(ConfOccupancyBins.value.begin(), ConfOccupancyBins.value.end()); - occ_bin_edges.erase(occ_bin_edges.begin()); - for (const auto& edge : occ_bin_edges) { - LOGF(info, "VARIABLE_WIDTH: occ_bin_edges = %f", edge); - } - } else { - int nbins = static_cast(ConfOccupancyBins.value[0]); - float xmin = static_cast(ConfOccupancyBins.value[1]); - float xmax = static_cast(ConfOccupancyBins.value[2]); - occ_bin_edges.resize(nbins + 1); - for (int i = 0; i < nbins + 1; i++) { - occ_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; - LOGF(info, "FIXED_WIDTH: occ_bin_edges[%d] = %f", i, occ_bin_edges[i]); - } - } + parseBinsVerbatim(confOccupancyBins, occ_bin_edges); emh1 = new MyEMH(ndepth); emh2 = new MyEMH(ndepth); @@ -281,82 +334,88 @@ struct PhotonHBT { o2::aod::pwgem::photonmeson::utils::eventhistogram::addEventHistograms(&fRegistry); DefineEMEventCut(); DefinePCMCut(); - addhistograms(); std::random_device seed_gen; engine = std::mt19937(seed_gen()); dist01 = std::uniform_int_distribution(0, 1); - fRegistry.add("Pair/mix/hDiffBC", "diff. global BC in mixed event;|BC_{current} - BC_{mixed}|", kTH1D, {{10001, -0.5, 10000.5}}, true); + fRegistry.add("Pair/mix/hDiffBC", + "diff. global BC in mixed event;|BC_{current} - BC_{mixed}|", + kTH1D, {{10001, -0.5, 10000.5}}, true); } template void initCCDB(TCollision const& collision) { - if (mRunNumber == collision.runNumber()) { + if (mRunNumber == collision.runNumber()) return; - } - - // // In case override, don't proceed, please - no CCDB access required - // if (d_bz_input > -990) { - // d_bz = d_bz_input; - // o2::parameters::GRPMagField grpmag; - // if (std::fabs(d_bz) > 1e-5) { - // grpmag.setL3Current(30000.f / (d_bz / 5.0f)); - // } - // mRunNumber = collision.runNumber(); - // return; - // } - - // auto run3grp_timestamp = collision.timestamp(); - // o2::parameters::GRPObject* grpo = 0x0; - // o2::parameters::GRPMagField* grpmag = 0x0; - // if (!skipGRPOquery) - // grpo = ccdb->getForTimeStamp(grpPath, run3grp_timestamp); - // if (grpo) { - // // Fetch magnetic field from ccdb for current collision - // d_bz = grpo->getNominalL3Field(); - // LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG"; - // } else { - // grpmag = ccdb->getForTimeStamp(grpmagPath, run3grp_timestamp); - // if (!grpmag) { - // LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grp_timestamp; - // } - // // Fetch magnetic field from ccdb for current collision - // d_bz = std::lround(5.f * grpmag->getL3Current() / 30000.f); - // LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG"; - // } mRunNumber = collision.runNumber(); } void addhistograms() { - // o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms<-1>(&fRegistry); static constexpr std::string_view qvec_det_names[6] = {"FT0M", "FT0A", "FT0C", "BTot", "BPos", "BNeg"}; - fRegistry.add("Event/before/hEP2_CentFT0C_forMix", Form("2nd harmonics event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[cfgEP2Estimator_for_Mix].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 event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[cfgEP2Estimator_for_Mix].data()), kTH2F, {{110, 0, 110}, {180, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}}, false); - - // pair info - const AxisSpec axis_kt{ConfKtBins, "k_{T} (GeV/c)"}; - const AxisSpec axis_qinv{ConfQBins, "q_{inv} (GeV/c)"}; - const AxisSpec axis_qabs_lcms{ConfQBins, "|#bf{q}|^{LCMS} (GeV/c)"}; - const AxisSpec axis_qout{ConfQBins, "q_{out} (GeV/c)"}; // qout does not change between LAB and LCMS frame - const AxisSpec axis_qside{ConfQBins, "q_{side} (GeV/c)"}; // qside does not change between LAB and LCMS frame - const AxisSpec axis_qlong{ConfQBins, "q_{long} (GeV/c)"}; - - if (cfgDo3D) { // 3D + fRegistry.add("Event/before/hEP2_CentFT0C_forMix", + Form("2nd harmonics event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[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 event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[cfgEp2EstimatorForMix].data()), + kTH2F, {{110, 0, 110}, {180, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}}, false); + + const AxisSpec axis_kt{confKtBins, "k_{T} (GeV/c)"}; + const AxisSpec axis_qinv{confQBins, "q_{inv} (GeV/c)"}; + const AxisSpec axis_qabs_lcms{confQBins, "|#bf{q}|^{LCMS} (GeV/c)"}; + const AxisSpec axis_qout{confQBins, "q_{out} (GeV/c)"}; + const AxisSpec axis_qside{confQBins, "q_{side} (GeV/c)"}; + const AxisSpec axis_qlong{confQBins, "q_{long} (GeV/c)"}; + + const AxisSpec axisPt{confPtBins, "p_{T} (GeV/c)"}; + const AxisSpec axisEta{confEtaBins, "#eta"}; + const AxisSpec axisPhi{confPhiBins, "#phi (rad)"}; + const AxisSpec axisDeltaEta{confDeltaEtaBins, "#Delta#eta"}; + const AxisSpec axisDeltaPhi{confDeltaPhiBins, "#Delta#phi (rad)"}; + const AxisSpec axisEllipseVal{confEllipseValBins, "Ellipse val"}; + const AxisSpec axisCosTheta{confCosThetaBins, "cos(#theta*)"}; + const AxisSpec axisOpeningAngle{confOpeningAngleBins, "#alpha (rad)"}; + + const AxisSpec axisV0Combo{4, -0.5f, 3.5f, "V0 combo (0=Incl,1=II,2=IT,3=TT)"}; + const AxisSpec axisPairCombo{7, -0.5f, 6.5f, "Pair combo (0=Incl,1=II-II,2=II-IT,3=II-TT,4=IT-IT,5=IT-TT,6=TT-TT)"}; + const AxisSpec axisRxy{confRxyBins, "R_{xy} (cm)"}; + const AxisSpec axisPsiPair{confPsiPairBins, "#psi_{pair} (rad)"}; + + // ── Single-photon QA ───────────────────────────────────────────────────── + + fRegistry.add("SinglePhoton/hPt", "V0 photon p_{T};p_{T} (GeV/c);V0 combo", kTH2F, {axisPt, axisV0Combo}, true); + fRegistry.add("SinglePhoton/hEta", "V0 photon #eta;#eta;V0 combo", kTH2F, {axisEta, axisV0Combo}, true); + fRegistry.add("SinglePhoton/hPhi", "V0 photon #phi;#phi (rad);V0 combo", kTH2F, {axisPhi, axisV0Combo}, true); + fRegistry.add("SinglePhoton/hEtaVsPhi", "V0 photon acceptance;#phi (rad);#eta;V0 combo", kTH3F, {axisPhi, axisEta, axisV0Combo}, true); + fRegistry.add("SinglePhoton/hRxy", "Conversion R_{xy};R_{xy} (cm);V0 combo", kTH2F, {axisRxy, axisV0Combo}, true); + fRegistry.add("SinglePhoton/hPsiPair", "#psi_{pair};#psi_{pair} (rad);V0 combo", kTH2F, {axisPsiPair, axisV0Combo}, true); + + if (cfgDo3D) { fRegistry.add("Pair/same/hs_3d", "diphoton correlation 3D LCMS", kTHnSparseD, {axis_qout, axis_qside, axis_qlong, axis_kt}, true); - } else { // 1D - if (cfgUseLCMS) { + } else { + if (cfgUseLcms) fRegistry.add("Pair/same/hs_1d", "diphoton correlation 1D LCMS", kTHnSparseD, {axis_qabs_lcms, axis_kt}, true); - } else { + else fRegistry.add("Pair/same/hs_1d", "diphoton correlation 1D", kTHnSparseD, {axis_qinv, axis_kt}, true); - } } + if constexpr (pairtype == ggHBTPairType::kPCMPCM) + fRegistry.add("Pair/same/hDeltaRCosOA", "distance between 2 conversion points;#Deltar/cos(#theta_{op}/2) (cm)", kTH1D, {{100, 0, 100}}, true); - if constexpr (pairtype == ggHBTPairType::kPCMPCM) { - fRegistry.add("Pair/same/hDeltaRCosOA", "distance between 2 conversion points;#Deltar/cos(#theta_{op}/2) (cm)", kTH1D, {{100, 0, 100}}, true); // dr/cosOA of conversion points + fRegistry.add("Pair/same/hKt", "k_{T};k_{T} (GeV/c);counts", kTH1F, {axis_kt}, true); + + for (const auto& step : {"Before", "After"}) { + const std::string s = std::string("Pair/same/QA/") + step + "/"; + + fRegistry.add((s + "hDeltaEta").c_str(), "#Delta#eta;#Delta#eta;Pair combo", kTH2F, {axisDeltaEta, axisPairCombo}, true); + fRegistry.add((s + "hDeltaPhi").c_str(), "#Delta#phi;#Delta#phi (rad);Pair combo", kTH2F, {axisDeltaPhi, axisPairCombo}, true); + fRegistry.add((s + "hDEtaDPhi").c_str(), "#Delta#eta vs #Delta#phi;#Delta#eta;#Delta#phi (rad);Pair combo", kTH3F, {axisDeltaEta, axisDeltaPhi, axisPairCombo}, true); + fRegistry.add((s + "hDeltaEtaVsPairEta").c_str(), "#Delta#eta vs #LT#eta#GT_{pair};#LT#eta#GT_{pair};#Delta#eta;Pair combo", kTH3F, {axisEta, axisDeltaEta, axisPairCombo}, true); + fRegistry.add((s + "hCosTheta").c_str(), "cos(#theta*);cos(#theta*);Pair combo", kTH2F, {axisCosTheta, axisPairCombo}, true); + fRegistry.add((s + "hOpeningAngle").c_str(), "Opening angle;#alpha (rad);Pair combo", kTH2F, {axisOpeningAngle, axisPairCombo}, true); + fRegistry.add((s + "hEllipseVal").c_str(), "Ellipse value;value;Pair combo", kTH2F, {axisEllipseVal, axisPairCombo}, true); } fRegistry.addClone("Pair/same/", "Pair/mix/"); @@ -386,179 +445,223 @@ struct PhotonHBT { void DefinePCMCut() { fV0PhotonCut = V0PhotonCut("fV0PhotonCut", "fV0PhotonCut"); - - // for v0 - fV0PhotonCut.SetV0PtRange(pcmcuts.cfg_min_pt_v0, 1e10f); - fV0PhotonCut.SetV0EtaRange(-pcmcuts.cfg_max_eta_v0, +pcmcuts.cfg_max_eta_v0); - fV0PhotonCut.SetMinCosPA(pcmcuts.cfg_min_cospa); - fV0PhotonCut.SetMaxPCA(pcmcuts.cfg_max_pca); - fV0PhotonCut.SetMaxChi2KF(pcmcuts.cfg_max_chi2kf); - fV0PhotonCut.SetRxyRange(pcmcuts.cfg_min_v0radius, pcmcuts.cfg_max_v0radius); - fV0PhotonCut.SetAPRange(pcmcuts.cfg_max_alpha_ap, pcmcuts.cfg_max_qt_ap); - fV0PhotonCut.RejectITSib(pcmcuts.cfg_reject_v0_on_itsib); - - // for track - fV0PhotonCut.SetMinNClustersTPC(pcmcuts.cfg_min_ncluster_tpc); - fV0PhotonCut.SetMinNCrossedRowsTPC(pcmcuts.cfg_min_ncrossedrows); + fV0PhotonCut.SetV0PtRange(pcmcuts.cfgMinPtV0, 1e10f); + fV0PhotonCut.SetV0EtaRange(-pcmcuts.cfgMaxEtaV0, +pcmcuts.cfgMaxEtaV0); + fV0PhotonCut.SetMinCosPA(pcmcuts.cfgMinCospa); + fV0PhotonCut.SetMaxPCA(pcmcuts.cfgMaxPca); + fV0PhotonCut.SetMaxChi2KF(pcmcuts.cfgMaxChi2Kf); + fV0PhotonCut.SetRxyRange(pcmcuts.cfgMinV0Radius, pcmcuts.cfgMaxV0Radius); + fV0PhotonCut.SetAPRange(pcmcuts.cfgMaxAlphaAp, pcmcuts.cfgMaxQtAp); + fV0PhotonCut.RejectITSib(pcmcuts.cfgRejectV0OnItsib); + fV0PhotonCut.SetMinNClustersTPC(pcmcuts.cfgMinNclusterTpc); + fV0PhotonCut.SetMinNCrossedRowsTPC(pcmcuts.cfgMinNcrossedrows); fV0PhotonCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); - fV0PhotonCut.SetMaxFracSharedClustersTPC(pcmcuts.cfg_max_frac_shared_clusters_tpc); - fV0PhotonCut.SetChi2PerClusterTPC(0.0, pcmcuts.cfg_max_chi2tpc); - fV0PhotonCut.SetTPCNsigmaElRange(pcmcuts.cfg_min_TPCNsigmaEl, pcmcuts.cfg_max_TPCNsigmaEl); - fV0PhotonCut.SetChi2PerClusterITS(-1e+10, pcmcuts.cfg_max_chi2its); - fV0PhotonCut.SetDisableITSonly(pcmcuts.cfg_disable_itsonly_track); - fV0PhotonCut.SetDisableTPConly(pcmcuts.cfg_disable_tpconly_track); + fV0PhotonCut.SetMaxFracSharedClustersTPC(pcmcuts.cfgMaxFracSharedClustersTpc); + fV0PhotonCut.SetChi2PerClusterTPC(0.0, pcmcuts.cfgMaxChi2Tpc); + fV0PhotonCut.SetTPCNsigmaElRange(pcmcuts.cfgMinTpcNsigmaEl, pcmcuts.cfgMaxTpcNsigmaEl); + fV0PhotonCut.SetChi2PerClusterITS(-1e+10, pcmcuts.cfgMaxChi2Its); + fV0PhotonCut.SetDisableITSonly(pcmcuts.cfgDisableItsonlyTrack); + fV0PhotonCut.SetDisableTPConly(pcmcuts.cfgDisableTpconlyTrack); fV0PhotonCut.SetNClustersITS(0, 7); fV0PhotonCut.SetMeanClusterSizeITSob(0.0, 16.0); - fV0PhotonCut.SetRequireITSTPC(pcmcuts.cfg_require_v0_with_itstpc); - fV0PhotonCut.SetRequireITSonly(pcmcuts.cfg_require_v0_with_itsonly); - fV0PhotonCut.SetRequireTPConly(pcmcuts.cfg_require_v0_with_tpconly); + fV0PhotonCut.SetRequireITSTPC(pcmcuts.cfgRequireV0WithItstpc); + fV0PhotonCut.SetRequireITSonly(pcmcuts.cfgRequireV0WithItsonly); + fV0PhotonCut.SetRequireTPConly(pcmcuts.cfgRequireV0WithTpconly); } template - void fillPairHistogram(TCollision const&, const ROOT::Math::PtEtaPhiMVector v1, const ROOT::Math::PtEtaPhiMVector v2, const float weight = 1.f) + void fillPairHistogram(TCollision const&, + const ROOT::Math::PtEtaPhiMVector v1, + const ROOT::Math::PtEtaPhiMVector v2, + const float weight = 1.f) { - float rndm = std::pow(-1, dist01(engine) % 2); // +1 or -1 to randomize order between 1 and 2. - // Lab. frame + float rndm = std::pow(-1, dist01(engine) % 2); ROOT::Math::PtEtaPhiMVector q12 = (v1 - v2) * rndm; ROOT::Math::PtEtaPhiMVector k12 = 0.5 * (v1 + v2); - float qinv = -q12.M(); // for identical particles -> qinv = 2 x kstar + float qinv = -q12.M(); float kt = k12.Pt(); - - ROOT::Math::XYZVector uv_out(k12.Px() / k12.Pt(), k12.Py() / k12.Pt(), 0); // unit vector for out. i.e. parallel to kt - ROOT::Math::XYZVector uv_long(0, 0, 1); // unit vector for long, beam axis - ROOT::Math::XYZVector uv_side = uv_out.Cross(uv_long); // unit vector for side - - ROOT::Math::PxPyPzEVector v1_cartesian(v1); - ROOT::Math::PxPyPzEVector v2_cartesian(v2); - ROOT::Math::PxPyPzEVector q12_cartesian = (v1_cartesian - v2_cartesian) * rndm; - float beta = (v1 + v2).Beta(); - // float beta_x = beta * std::cos((v1 + v2).Phi()) * std::sin((v1 + v2).Theta()); - // float beta_y = beta * std::sin((v1 + v2).Phi()) * std::sin((v1 + v2).Theta()); - float beta_z = beta * std::cos((v1 + v2).Theta()); - - // longitudinally co-moving system (LCMS) - ROOT::Math::Boost bst_z(0, 0, -beta_z); // Boost supports only PxPyPzEVector - ROOT::Math::PxPyPzEVector q12_lcms = bst_z(q12_cartesian); - ROOT::Math::XYZVector q_3d_lcms = q12_lcms.Vect(); // 3D q vector in LCMS + 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); + ROOT::Math::PxPyPzEVector q12c = (v1c - v2c) * rndm; + float beta_z = (v1 + v2).Beta() * std::cos((v1 + v2).Theta()); + ROOT::Math::Boost bst_z(0, 0, -beta_z); + ROOT::Math::PxPyPzEVector q12_lcms = bst_z(q12c); + ROOT::Math::XYZVector q_3d_lcms = q12_lcms.Vect(); + float qabs_lcms = q_3d_lcms.R(); float qout_lcms = q_3d_lcms.Dot(uv_out); float qside_lcms = q_3d_lcms.Dot(uv_side); float qlong_lcms = q_3d_lcms.Dot(uv_long); - float qabs_lcms = q_3d_lcms.R(); + if (cfgDo3D) + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_3d"), + std::fabs(qout_lcms), std::fabs(qside_lcms), std::fabs(qlong_lcms), kt, weight); + else { + if (cfgUseLcms) + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_1d"), qabs_lcms, kt, weight); + else + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_1d"), qinv, kt, weight); + } - // float qabs_lcms_tmp = std::sqrt(std::pow(qout_lcms, 2) + std::pow(qside_lcms, 2) + std::pow(qlong_lcms, 2)); - // LOGF(info, "qabs_lcms = %f, qabs_lcms_tmp = %f", qabs_lcms, qabs_lcms_tmp); - - // // pair rest frame (PRF) - // ROOT::Math::Boost boostPRF = ROOT::Math::Boost(-beta_x, -beta_y, -beta_z); - // ROOT::Math::PxPyPzEVector v1_prf = boostPRF(v1_cartesian); - // ROOT::Math::PxPyPzEVector v2_prf = boostPRF(v2_cartesian); - // ROOT::Math::PxPyPzEVector rel_k = (v1_prf - v2_prf) * rndm; - // float kstar = 0.5 * rel_k.P(); - // // LOGF(info, "qabs_lcms = %f, qinv = %f, kstar = %f", qabs_lcms, qinv, kstar); - - // ROOT::Math::PxPyPzEVector v1_lcms_cartesian = bst_z(v1_cartesian); - // ROOT::Math::PxPyPzEVector v2_lcms_cartesian = bst_z(v2_cartesian); - // ROOT::Math::PxPyPzEVector q12_lcms_cartesian = bst_z(q12_cartesian); - // LOGF(info, "q12.Pz() = %f, q12_cartesian.Pz() = %f", q12.Pz(), q12_cartesian.Pz()); - // LOGF(info, "v1.Pz() = %f, v2.Pz() = %f", v1.Pz(), v2.Pz()); - // LOGF(info, "v1_lcms_cartesian.Pz() = %f, v2_lcms_cartesian.Pz() = %f", v1_lcms_cartesian.Pz(), v2_lcms_cartesian.Pz()); - // LOGF(info, "q12_lcms_cartesian.Pz() = %f", q12_lcms_cartesian.Pz()); - // LOGF(info, "q_3d_lcms.Dot(uv_out) = %f, q_3d_lcms.Dot(uv_side) = %f, q_3d.Dot(uv_out) = %f, q_3d.Dot(uv_side) = %f", q_3d_lcms.Dot(uv_out), q_3d_lcms.Dot(uv_side), q_3d.Dot(uv_out), q_3d.Dot(uv_side)); - // LOGF(info, "q12_lcms.Pz() = %f, q_3d_lcms.Dot(uv_long) = %f", q12_lcms.Pz(), q_3d_lcms.Dot(uv_long)); - // ROOT::Math::PxPyPzEVector q12_lcms_tmp = bst_z(v1_cartesian) - bst_z(v2_cartesian); - // LOGF(info, "q12_lcms.Px() = %f, q12_lcms.Py() = %f, q12_lcms.Pz() = %f, q12_lcms_tmp.Px() = %f, q12_lcms_tmp.Py() = %f, q12_lcms_tmp.Pz() = %f", q12_lcms.Px(), q12_lcms.Py(), q12_lcms.Pz(), q12_lcms_tmp.Px(), q12_lcms_tmp.Py(), q12_lcms_tmp.Pz()); - // float qabs_lcms_tmp = q12_lcms.P(); - // LOGF(info, "qabs_lcms = %f, qabs_lcms_tmp = %f", qabs_lcms, qabs_lcms_tmp); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hKt"), kt, weight); + } - if (cfgDo3D) { - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_3d"), std::fabs(qout_lcms), std::fabs(qside_lcms), std::fabs(qlong_lcms), kt, weight); // qosl can be [-inf, +inf] and CF is symmetric for pos and neg qosl. To reduce stat. unc. absolute value is taken here. + template + inline void fillPairQAStep(float deta, float dphi, float pairEta, + float cosTheta, float openingAngle, float ellVal, + int pairComboIdx) + { + if (!doPairQa) + return; + + if constexpr (IsBefore) { + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDeltaEta"), deta, 0.f); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDeltaPhi"), dphi, 0.f); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDEtaDPhi"), deta, dphi, 0.f); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDeltaEtaVsPairEta"), pairEta, deta, 0.f); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hCosTheta"), cosTheta, 0.f); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hOpeningAngle"), openingAngle, 0.f); + if (ellVal >= 0.f) + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hEllipseVal"), ellVal, 0.f); + if (pairComboIdx > 0) { + const float pcb = float(pairComboIdx); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDeltaEta"), deta, pcb); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDeltaPhi"), dphi, pcb); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDEtaDPhi"), deta, dphi, pcb); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDeltaEtaVsPairEta"), pairEta, deta, pcb); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hCosTheta"), cosTheta, pcb); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hOpeningAngle"), openingAngle, pcb); + if (ellVal >= 0.f) + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hEllipseVal"), ellVal, pcb); + } } else { - if (cfgUseLCMS) { - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_1d"), qabs_lcms, kt, weight); - } else { - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_1d"), qinv, kt, weight); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDeltaEta"), deta, 0.f); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDeltaPhi"), dphi, 0.f); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDEtaDPhi"), deta, dphi, 0.f); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDeltaEtaVsPairEta"), pairEta, deta, 0.f); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hCosTheta"), cosTheta, 0.f); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hOpeningAngle"), openingAngle, 0.f); + if (ellVal >= 0.f) + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hEllipseVal"), ellVal, 0.f); + if (pairComboIdx > 0) { + const float pcb = float(pairComboIdx); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDeltaEta"), deta, pcb); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDeltaPhi"), dphi, pcb); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDEtaDPhi"), deta, dphi, pcb); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDeltaEtaVsPairEta"), pairEta, deta, pcb); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hCosTheta"), cosTheta, pcb); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hOpeningAngle"), openingAngle, pcb); + if (ellVal >= 0.f) + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hEllipseVal"), ellVal, pcb); } } } - template - void runPairing(TCollisions const& collisions, TPhotons1 const& photons1, TPhotons2 const& photons2, TSubInfos1 const&, TSubInfos2 const&, TPreslice1 const& perCollision1, TPreslice2 const& perCollision2, TCut1 const& cut1, TCut2 const& cut2) + // =========================================================================== + // runPairing + // =========================================================================== + + template + void runPairing(TCollisions const& collisions, + TPhotons1 const& photons1, TPhotons2 const& photons2, + TSubInfos1 const&, TSubInfos2 const&, + TPreslice1 const& perCollision1, TPreslice2 const& perCollision2, + TCut1 const& cut1, TCut2 const& cut2) { for (const auto& collision : collisions) { initCCDB(collision); int ndiphoton = 0; + const float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; - if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { + if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) continue; - } - const float eventplanes_2_for_mix[6] = {collision.ep2ft0m(), collision.ep2ft0a(), collision.ep2ft0c(), collision.ep2btot(), collision.ep2bpos(), collision.ep2bneg()}; - float ep2 = eventplanes_2_for_mix[cfgEP2Estimator_for_Mix]; - fRegistry.fill(HIST("Event/before/hEP2_CentFT0C_forMix"), collision.centFT0C(), ep2); + const float ep2_arr[6] = {collision.ep2ft0m(), collision.ep2ft0a(), collision.ep2ft0c(), + collision.ep2btot(), collision.ep2bpos(), collision.ep2bneg()}; + const float ep3_arr[6] = {collision.ep3ft0m(), collision.ep3ft0a(), collision.ep3ft0c(), + collision.ep3btot(), collision.ep3bpos(), collision.ep3bneg()}; + float ep2 = ep2_arr[cfgEp2EstimatorForMix]; + float ep3 = ep3_arr[cfgEp3EstimatorForMix]; + + 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)) { + if (!fEMEventCut.IsSelected(collision)) continue; - } o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(&fRegistry, collision, 1.f); fRegistry.fill(HIST("Event/after/hEP2_CentFT0C_forMix"), collision.centFT0C(), ep2); - int zbin = lower_bound(zvtx_bin_edges.begin(), zvtx_bin_edges.end(), collision.posZ()) - zvtx_bin_edges.begin() - 1; - if (zbin < 0) { - zbin = 0; - } else if (static_cast(zvtx_bin_edges.size()) - 2 < zbin) { - zbin = static_cast(zvtx_bin_edges.size()) - 2; - } - - float centrality = centralities[cfgCentEstimator]; - int centbin = lower_bound(cent_bin_edges.begin(), cent_bin_edges.end(), centrality) - cent_bin_edges.begin() - 1; - if (centbin < 0) { - centbin = 0; - } else if (static_cast(cent_bin_edges.size()) - 2 < centbin) { - centbin = static_cast(cent_bin_edges.size()) - 2; - } - - int epbin = lower_bound(ep_bin_edges.begin(), ep_bin_edges.end(), ep2) - ep_bin_edges.begin() - 1; - if (epbin < 0) { - epbin = 0; - } else if (static_cast(ep_bin_edges.size()) - 2 < epbin) { - epbin = static_cast(ep_bin_edges.size()) - 2; - } - - int occbin = -1; - if (cfgOccupancyEstimator == 0) { - occbin = lower_bound(occ_bin_edges.begin(), occ_bin_edges.end(), collision.ft0cOccupancyInTimeRange()) - occ_bin_edges.begin() - 1; - } else if (cfgOccupancyEstimator == 1) { - occbin = lower_bound(occ_bin_edges.begin(), occ_bin_edges.end(), collision.trackOccupancyInTimeRange()) - occ_bin_edges.begin() - 1; - } else { - occbin = lower_bound(occ_bin_edges.begin(), occ_bin_edges.end(), collision.ft0cOccupancyInTimeRange()) - occ_bin_edges.begin() - 1; - } - - if (occbin < 0) { - occbin = 0; - } else if (static_cast(occ_bin_edges.size()) - 2 < occbin) { - occbin = static_cast(occ_bin_edges.size()) - 2; - } - - // LOGF(info, "collision.globalIndex() = %d, collision.posZ() = %f, centrality = %f, ep2 = %f, collision.trackOccupancyInTimeRange() = %d, zbin = %d, centbin = %d, epbin = %d, occbin = %d", collision.globalIndex(), collision.posZ(), centrality, ep2, collision.trackOccupancyInTimeRange(), zbin, centbin, epbin, occbin); - - auto key_bin = std::make_tuple(zbin, centbin, epbin, occbin); + auto clampBin = [](int b, int nmax) { + return (b < 0) ? 0 : (b > nmax ? nmax : b); + }; + auto binOf = [&](const std::vector& edges, float val) { + int b = static_cast(std::lower_bound(edges.begin(), edges.end(), val) - edges.begin()) - 1; + return clampBin(b, static_cast(edges.size()) - 2); + }; + + int zbin = binOf(zvtx_bin_edges, collision.posZ()); + int centbin = binOf(cent_bin_edges, centralities[cfgCentEstimator]); + int ep2bin = binOf(ep2_bin_edges, ep2); + int ep3bin = binOf(ep3_bin_edges, ep3); + int occbin = binOf(occ_bin_edges, + cfgOccupancyEstimator == 1 + ? static_cast(collision.trackOccupancyInTimeRange()) + : collision.ft0cOccupancyInTimeRange()); + + auto key_bin = std::make_tuple(zbin, centbin, ep2bin, ep3bin, occbin); auto key_df_collision = std::make_pair(ndf, collision.globalIndex()); if constexpr (pairtype == ggHBTPairType::kPCMPCM) { + auto photons1_coll = photons1.sliceBy(perCollision1, collision.globalIndex()); auto photons2_coll = photons2.sliceBy(perCollision2, collision.globalIndex()); + + // ── Single-photon QA ───────────────────────────────────────────────── + if (doSinglePhotonQa) { + for (const auto& g : photons1_coll) { + if (!cut1.template IsSelected(g)) + continue; + const int ci = classifyV0ComboIdx(g); // 0=Other/1/2/3 + fRegistry.fill(HIST("SinglePhoton/hPt"), g.pt(), 0.f); + fRegistry.fill(HIST("SinglePhoton/hEta"), g.eta(), 0.f); + fRegistry.fill(HIST("SinglePhoton/hPhi"), g.phi(), 0.f); + fRegistry.fill(HIST("SinglePhoton/hEtaVsPhi"), g.phi(), g.eta(), 0.f); + fRegistry.fill(HIST("SinglePhoton/hRxy"), std::hypot(g.vx(), g.vy()), 0.f); + if constexpr (requires { g.psipair(); }) { + fRegistry.fill(HIST("SinglePhoton/hPsiPair"), g.psipair(), 0.f); + } + if (ci > 0) { + const float cb = float(ci); + fRegistry.fill(HIST("SinglePhoton/hPt"), g.pt(), cb); + fRegistry.fill(HIST("SinglePhoton/hEta"), g.eta(), cb); + fRegistry.fill(HIST("SinglePhoton/hPhi"), g.phi(), cb); + fRegistry.fill(HIST("SinglePhoton/hEtaVsPhi"), g.phi(), g.eta(), cb); + fRegistry.fill(HIST("SinglePhoton/hRxy"), std::hypot(g.vx(), g.vy()), cb); + if constexpr (requires { g.psipair(); }) { + fRegistry.fill(HIST("SinglePhoton/hPsiPair"), g.psipair(), cb); + } + } + } + } + + // ── Same-event pair loop ────────────────────────────────────────────── for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photons1_coll, photons2_coll))) { - if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { + 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()) { // never happens. only for protection. + if (pos1.trackId() == pos2.trackId() || ele1.trackId() == ele2.trackId()) continue; - } + + const int c1 = classifyV0ComboIdx(g1); + const int c2 = classifyV0ComboIdx(g2); + const int pcb = pairComboBin(c1, c2); // 0 if either is Other ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); @@ -566,17 +669,35 @@ struct PhotonHBT { 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)); ROOT::Math::XYZVector cp1(g1.vx(), g1.vy(), g1.vz()); ROOT::Math::XYZVector cp2(g2.vx(), g2.vy(), g2.vz()); - float opa = std::acos(cp1.Dot(cp2) / (std::sqrt(cp1.Mag2()) * std::sqrt(cp2.Mag2()))); // opening angle between 2 conversion points + float opa = std::acos(std::clamp(static_cast(cp1.Dot(cp2) / (std::sqrt(cp1.Mag2()) * std::sqrt(cp2.Mag2()))), CosineMin, CosineMax)); o2::math_utils::bringTo02Pi(opa); - if (opa > o2::constants::math::PI) { + if (opa > o2::constants::math::PI) opa -= o2::constants::math::PI; - } - float cosOA = std::cos(opa / 2.f); - if (dr / cosOA < ggpaircuts.cfgMinDR_CosOA) { + float cosOA = std::cos(opa / HalfAngle); + if (dr / cosOA < ggpaircuts.cfgMinDrCosOa) continue; - } fRegistry.fill(HIST("Pair/same/hDeltaRCosOA"), dr / cosOA); + float deta = g1.eta() - g2.eta(); + float dphi = normDphi(g1.phi() - g2.phi()); + float pairEta = PairEtaWeight * (g1.eta() + g2.eta()); + float cosTheta = std::fabs(computeCosTheta(v1, v2)); + const float sE = ggpaircuts.cfgEllipseSigEta.value; + const float sP = ggpaircuts.cfgEllipseSigPhi.value; + float ellVal = (sE > MinFloatTol && sP > MinFloatTol) + ? (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP) + : -1.f; + + // QA Before ellipse cut + fillPairQAStep<0, true>(deta, dphi, pairEta, cosTheta, opa, ellVal, pcb); + + // Ellipse cut + if (isInsideEllipse(deta, dphi)) + continue; + + // QA After ellipse cut + fillPairQAStep<0, false>(deta, dphi, pairEta, cosTheta, opa, ellVal, pcb); + fillPairHistogram<0>(collision, v1, v2, 1.f); ndiphoton++; @@ -592,42 +713,33 @@ struct PhotonHBT { emh1->AddTrackToEventPool(key_df_collision, g2tmp); used_photonIds_per_col.emplace_back(g2.globalIndex()); } - } // end of pairing loop + } // end same-event pair loop } used_photonIds_per_col.clear(); used_photonIds_per_col.shrink_to_fit(); - // event mixing - if (!cfgDoMix || !(ndiphoton > 0)) { + // ── Mixed-event loop ──────────────────────────────────────────────────── + if (!cfgDoMix || !(ndiphoton > 0)) continue; - } - // make a vector of selected photons in this collision. auto selected_photons1_in_this_event = emh1->GetTracksPerCollision(key_df_collision); - auto selected_photons2_in_this_event = emh2->GetTracksPerCollision(key_df_collision); - auto collisionIds1_in_mixing_pool = emh1->GetCollisionIdsFromEventPool(key_bin); - auto collisionIds2_in_mixing_pool = emh2->GetCollisionIdsFromEventPool(key_bin); if constexpr (pairtype == ggHBTPairType::kPCMPCM) { for (const auto& mix_dfId_collisionId : collisionIds1_in_mixing_pool) { int mix_dfId = mix_dfId_collisionId.first; int64_t mix_collisionId = mix_dfId_collisionId.second; - - if (collision.globalIndex() == mix_collisionId && ndf == mix_dfId) { // this never happens. only protection. + if (collision.globalIndex() == mix_collisionId && ndf == mix_dfId) continue; - } auto globalBC_mix = map_mixed_eventId_to_globalBC[mix_dfId_collisionId]; uint64_t diffBC = std::max(collision.globalBC(), globalBC_mix) - std::min(collision.globalBC(), globalBC_mix); fRegistry.fill(HIST("Pair/mix/hDiffBC"), diffBC); - if (diffBC < ndiff_bc_mix) { + if (diffBC < ndiffBcMix) continue; - } auto photons1_from_event_pool = emh1->GetTracksPerCollision(mix_dfId_collisionId); - // LOGF(info, "Do event mixing: current event (%d, %d), ngamma = %d | event pool (%d, %d), ngamma = %d", ndf, collision.globalIndex(), selected_photons1_in_this_event.size(), mix_dfId, mix_collisionId, photons1_from_event_pool.size()); for (const auto& g1 : selected_photons1_in_this_event) { for (const auto& g2 : photons1_from_event_pool) { @@ -637,21 +749,36 @@ struct PhotonHBT { 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)); ROOT::Math::XYZVector cp1(g1.vx(), g1.vy(), g1.vz()); ROOT::Math::XYZVector cp2(g2.vx(), g2.vy(), g2.vz()); - float opa = std::acos(cp1.Dot(cp2) / (std::sqrt(cp1.Mag2()) * std::sqrt(cp2.Mag2()))); // opening angle between 2 conversion points + float opa = std::acos(std::clamp(static_cast(cp1.Dot(cp2) / (std::sqrt(cp1.Mag2()) * std::sqrt(cp2.Mag2()))), CosineMin, CosineMax)); o2::math_utils::bringTo02Pi(opa); - if (opa > o2::constants::math::PI) { + if (opa > o2::constants::math::PI) opa -= o2::constants::math::PI; - } - float cosOA = std::cos(opa / 2.f); - if (dr / cosOA < ggpaircuts.cfgMinDR_CosOA) { + float cosOA = std::cos(opa / HalfAngle); + if (dr / cosOA < ggpaircuts.cfgMinDrCosOa) continue; - } fRegistry.fill(HIST("Pair/mix/hDeltaRCosOA"), dr / cosOA); + float deta = g1.eta() - g2.eta(); + float dphi = normDphi(g1.phi() - g2.phi()); + float pairEta = PairEtaWeight * (g1.eta() + g2.eta()); + float cosTheta = std::fabs(computeCosTheta(v1, v2)); + const float sE = ggpaircuts.cfgEllipseSigEta.value; + const float sP = ggpaircuts.cfgEllipseSigPhi.value; + float ellVal = (sE > MinFloatTol && sP > MinFloatTol) + ? (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP) + : -1.f; + + fillPairQAStep<1, true>(deta, dphi, pairEta, cosTheta, opa, ellVal, 0); + + if (isInsideEllipse(deta, dphi)) + continue; + + fillPairQAStep<1, false>(deta, dphi, pairEta, cosTheta, opa, ellVal, 0); + fillPairHistogram<1>(collision, v1, v2, 1.f); } } - } // end of loop over mixed event pool + } // end mixed event pool loop } if (ndiphoton > 0) { @@ -659,21 +786,27 @@ struct PhotonHBT { emh2->AddCollisionIdAtLast(key_bin, key_df_collision); map_mixed_eventId_to_globalBC[key_df_collision] = collision.globalBC(); } - } // end of collision loop + } // 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 used_photonIds_per_col; // + std::vector used_photonIds_per_col; std::map, uint64_t> map_mixed_eventId_to_globalBC; SliceCache cache; Preslice perCollision_pcm = aod::v0photonkf::pmeventId; - Filter collisionFilter_centrality = (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 collisionFilter_occupancy_track = eventcuts.cfgTrackOccupancyMin <= o2::aod::evsel::trackOccupancyInTimeRange && o2::aod::evsel::trackOccupancyInTimeRange < eventcuts.cfgTrackOccupancyMax; - Filter collisionFilter_occupancy_ft0c = eventcuts.cfgFT0COccupancyMin <= o2::aod::evsel::ft0cOccupancyInTimeRange && o2::aod::evsel::ft0cOccupancyInTimeRange < eventcuts.cfgFT0COccupancyMax; + Filter collisionFilter_centrality = (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 collisionFilter_occupancy_track = eventcuts.cfgTrackOccupancyMin <= o2::aod::evsel::trackOccupancyInTimeRange && + o2::aod::evsel::trackOccupancyInTimeRange < eventcuts.cfgTrackOccupancyMax; + Filter collisionFilter_occupancy_ft0c = eventcuts.cfgFT0COccupancyMin <= o2::aod::evsel::ft0cOccupancyInTimeRange && + o2::aod::evsel::ft0cOccupancyInTimeRange < eventcuts.cfgFT0COccupancyMax; using FilteredMyCollisions = soa::Filtered; int ndf = 0; @@ -682,7 +815,8 @@ struct PhotonHBT { if constexpr (pairtype == ggHBTPairType::kPCMPCM) { auto v0photons = std::get<0>(std::tie(args...)); auto v0legs = std::get<1>(std::tie(args...)); - runPairing(collisions, v0photons, v0photons, v0legs, v0legs, perCollision_pcm, perCollision_pcm, fV0PhotonCut, fV0PhotonCut); + runPairing(collisions, v0photons, v0photons, v0legs, v0legs, + perCollision_pcm, perCollision_pcm, fV0PhotonCut, fV0PhotonCut); } ndf++; } From 2fb63882830a2d648e1c34c0ad04f4e5a0537f7b Mon Sep 17 00:00:00 2001 From: Stefanie Mrozinski Date: Mon, 16 Mar 2026 14:47:56 +0100 Subject: [PATCH 2/3] Fix formatting --- PWGEM/PhotonMeson/Core/PhotonHBT.h | 1 - 1 file changed, 1 deletion(-) diff --git a/PWGEM/PhotonMeson/Core/PhotonHBT.h b/PWGEM/PhotonMeson/Core/PhotonHBT.h index 364189d0fe7..c283622dae8 100644 --- a/PWGEM/PhotonMeson/Core/PhotonHBT.h +++ b/PWGEM/PhotonMeson/Core/PhotonHBT.h @@ -81,7 +81,6 @@ using MyV0Photon = MyV0Photons::iterator; template struct PhotonHBT { - // Single-photon: // 0 = Inclusive (all photons) // 1 = ITSTPC_ITSTPC — both legs have ITS+TPC From 43e113bf24f14740a09c6d2c7a50a5b27b1c21c0 Mon Sep 17 00:00:00 2001 From: Stefanie Mrozinski Date: Tue, 17 Mar 2026 11:57:20 +0100 Subject: [PATCH 3/3] Rewrite Header file + Add QA histograms --- PWGEM/PhotonMeson/Core/PhotonHBT.h | 828 -------------------- PWGEM/PhotonMeson/Tasks/CMakeLists.txt | 10 +- PWGEM/PhotonMeson/Tasks/PhotonHBTPCMPCM.cxx | 27 - PWGEM/PhotonMeson/Tasks/photonhbt.cxx | 784 ++++++++++++++++++ 4 files changed, 789 insertions(+), 860 deletions(-) delete mode 100644 PWGEM/PhotonMeson/Core/PhotonHBT.h delete mode 100644 PWGEM/PhotonMeson/Tasks/PhotonHBTPCMPCM.cxx create mode 100644 PWGEM/PhotonMeson/Tasks/photonhbt.cxx diff --git a/PWGEM/PhotonMeson/Core/PhotonHBT.h b/PWGEM/PhotonMeson/Core/PhotonHBT.h deleted file mode 100644 index c283622dae8..00000000000 --- a/PWGEM/PhotonMeson/Core/PhotonHBT.h +++ /dev/null @@ -1,828 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file PhotonHBT.h -/// \brief This code loops over v0 photons and makes pairs for photon HBT analysis. -/// \author Daiki Sekihata, daiki.sekihata@cern.ch - -#ifndef PWGEM_PHOTONMESON_CORE_PHOTONHBT_H_ -#define PWGEM_PHOTONMESON_CORE_PHOTONHBT_H_ - -#include "PWGEM/Dilepton/Utils/EMTrack.h" -#include "PWGEM/Dilepton/Utils/EventMixingHandler.h" -#include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" -#include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" -#include "PWGEM/PhotonMeson/DataModel/EventTables.h" -#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" -#include "PWGEM/PhotonMeson/Utils/EventHistograms.h" -// -#include "Common/DataModel/Centrality.h" -#include "Common/DataModel/EventSelection.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include // IWYU pragma: keep -#include -#include // IWYU pragma: keep -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace o2::aod::pwgem::photon::core::photonhbt -{ -enum class ggHBTPairType : int { - kPCMPCM = 0, -}; -} // namespace o2::aod::pwgem::photon::core::photonhbt - -using namespace o2; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) -using namespace o2::aod; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) -using namespace o2::framework; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) -using namespace o2::framework::expressions; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) -using namespace o2::soa; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) -using namespace o2::aod::pwgem::dilepton::utils; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) -using namespace o2::aod::pwgem::photon::core::photonhbt; // o2-linter: disable=using-directive (required by O2 framework, inherited from upstream PWGEM) - -using MyCollisions = soa::Join; -using MyCollision = MyCollisions::iterator; - -using MyV0Photons = soa::Join; -using MyV0Photon = MyV0Photons::iterator; - -template -struct PhotonHBT { - - // Single-photon: - // 0 = Inclusive (all photons) - // 1 = ITSTPC_ITSTPC — both legs have ITS+TPC - // 2 = ITSTPC_TPCOnly — one ITS+TPC, one TPC-only - // 3 = TPCOnly_TPCOnly — both legs TPC-only - // - // Pair combo. - // 0 = Inclusive (all recognised pairs) - // 1 = ITSTPC_ITSTPC × ITSTPC_ITSTPC - // 2 = ITSTPC_ITSTPC × ITSTPC_TPCOnly - // 3 = ITSTPC_ITSTPC × TPCOnly_TPCOnly - // 4 = ITSTPC_TPCOnly × ITSTPC_TPCOnly - // 5 = ITSTPC_TPCOnly × TPCOnly_TPCOnly - // 6 = TPCOnly_TPCOnly × TPCOnly_TPCOnly - - template - static inline int classifyV0ComboIdx(TGamma const& g) - { - auto pos = g.template posTrack_as(); - auto neg = g.template negTrack_as(); - const bool posII = pos.hasITS() && pos.hasTPC(); - const bool posTPC = !pos.hasITS() && pos.hasTPC(); - const bool negII = neg.hasITS() && neg.hasTPC(); - const bool negTPC = !neg.hasITS() && neg.hasTPC(); - if (posII && negII) - return 1; - if ((posII && negTPC) || (posTPC && negII)) - return 2; - if (posTPC && negTPC) - return 3; - return 0; - } - - static inline int pairComboBin(int c1, int c2) - { - if (c1 <= 0 || c2 <= 0) - return 0; - if (c1 > c2) - std::swap(c1, c2); - static constexpr int kTable[4][4] = { - {0, 0, 0, 0}, - {0, 1, 2, 3}, - {0, 2, 4, 5}, - {0, 3, 5, 6}}; - return kTable[c1][c2]; - } - - Configurable cfgDo3D{"cfgDo3D", false, "enable 3D analysis"}; - Configurable cfgEp2EstimatorForMix{"cfgEp2EstimatorForMix", 3, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5"}; - Configurable cfgEp3EstimatorForMix{"cfgEp3EstimatorForMix", 3, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5 (for psi3 mixing)"}; - Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; - Configurable cfgOccupancyEstimator{"cfgOccupancyEstimator", 0, "FT0C:0, Track:1"}; - Configurable cfgCentMin{"cfgCentMin", -1, "min. centrality"}; - Configurable cfgCentMax{"cfgCentMax", 999, "max. centrality"}; - Configurable maxY{"maxY", 0.8, "maximum rapidity for reconstructed particles"}; - 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 cfgUseLcms{"cfgUseLcms", true, "measure relative momentum in LCMS for 1D"}; - - ConfigurableAxis confVtxBins{"confVtxBins", {VARIABLE_WIDTH, -10.f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; - ConfigurableAxis confCentBins{"confCentBins", {VARIABLE_WIDTH, 0.f, 5.f, 10.f, 20.f, 30.f, 40.f, 50.f, 60.f, 70.f, 80.f, 90.f, 100.f, 999.f}, "Mixing bins - centrality"}; - ConfigurableAxis confEpBins{"confEpBins", {8, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}, "Mixing bins - psi2 event plane"}; - - ConfigurableAxis confEp3Bins{"confEp3Bins", {1, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}, "Mixing bins - psi3 event plane (set to 1 bin to disable)"}; - ConfigurableAxis confOccupancyBins{"confOccupancyBins", {VARIABLE_WIDTH, -1, 1e+10}, "Mixing bins - occupancy"}; - - ConfigurableAxis confQBins{"confQBins", {60, 0, +0.3f}, "q bins for output histograms"}; // o2-linter: disable=name/configurable (Q is a physics symbol for momentum transfer; cannot be lowercased without losing meaning) - ConfigurableAxis confKtBins{"confKtBins", {VARIABLE_WIDTH, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0}, "kT bins for output histograms"}; - - // QA axis configurables - 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 confDeltaEtaBins{"confDeltaEtaBins", {100, -0.5f, +0.5f}, "Delta-eta bins"}; - ConfigurableAxis confDeltaPhiBins{"confDeltaPhiBins", {100, -0.5f, +0.5f}, "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)"}; - - ConfigurableAxis confRxyBins{"confRxyBins", {90, 0.f, 90.f}, "Rxy bins (cm)"}; - - ConfigurableAxis confPsiPairBins{"confPsiPairBins", {100, -0.2f, 0.2f}, "psi_pair bins (rad)"}; - - // =========================================================================== - // QA flags - // =========================================================================== - - Configurable doPairQa{"doPairQa", true, "fill pair QA histograms (Before/After ellipse cut, with combo axis)"}; - Configurable doSinglePhotonQa{"doSinglePhotonQa", true, "fill single-photon QA histograms (pT, eta, phi with V0 combo axis)"}; - - // =========================================================================== - // Pair cuts - // =========================================================================== - - struct : ConfigurableGroup { - std::string prefix = "ggpaircut_group"; - Configurable cfgMinDrCosOa{"cfgMinDrCosOa", -1, "min. dr/cosOA for kPCMPCM"}; - Configurable cfgApplyEllipseCut{"cfgApplyEllipseCut", 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"}; - } ggpaircuts; - - EMPhotonEventCut fEMEventCut; - struct : ConfigurableGroup { - std::string prefix = "eventcut_group"; - Configurable cfgZvtxMin{"cfgZvtxMin", -10.f, "min. Zvtx"}; - Configurable cfgZvtxMax{"cfgZvtxMax", +10.f, "max. Zvtx"}; - Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8"}; - Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRequireNoTFB{"cfgRequireNoTFB", true, "require no TF border"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", true, "require no ITS ROF border"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup"}; - Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx FT0 vs PV"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgTrackOccupancyMin{"cfgTrackOccupancyMin", -2, "min. track occupancy"}; - Configurable cfgTrackOccupancyMax{"cfgTrackOccupancyMax", 1000000000, "max. track occupancy"}; - Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -2.f, "min. FT0C occupancy"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000.f, "max. FT0C occupancy"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRequireNoCollInTimeRangeStandard{"cfgRequireNoCollInTimeRangeStandard", false, "require no collision in time range standard"}; - Configurable cfgRequireNoCollInTimeRangeStrict{"cfgRequireNoCollInTimeRangeStrict", false, "require no collision in time range strict"}; - Configurable cfgRequireNoCollInITSROFStandard{"cfgRequireNoCollInITSROFStandard", false, "require no collision in time range standard"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRequireNoCollInITSROFStrict{"cfgRequireNoCollInITSROFStrict", false, "require no collision in time range strict"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRequireNoHighMultCollInPrevRof{"cfgRequireNoHighMultCollInPrevRof", false, "require no HM collision in previous ITS ROF"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRequireGoodITSLayer3{"cfgRequireGoodITSLayer3", false, "ITS layer 3 chips OK"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRequireGoodITSLayer0123{"cfgRequireGoodITSLayer0123", false, "ITS layers 0-3 chips OK"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRequireGoodITSLayersAll{"cfgRequireGoodITSLayersAll", false, "all ITS layers chips OK"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - } eventcuts; - - V0PhotonCut fV0PhotonCut; - struct : ConfigurableGroup { - std::string prefix = "pcmcut_group"; - Configurable cfgRequireV0WithItstpc{"cfgRequireV0WithItstpc", false, "flag to select V0s with ITS-TPC matched tracks"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRequireV0WithItsonly{"cfgRequireV0WithItsonly", false, "flag to select V0s with ITSonly tracks"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRequireV0WithTpconly{"cfgRequireV0WithTpconly", false, "flag to select V0s with TPConly tracks"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMinPtV0{"cfgMinPtV0", 0.1, "min pT for v0 photons at PV"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMaxEtaV0{"cfgMaxEtaV0", 0.8, "max eta for v0 photons at PV"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMinV0Radius{"cfgMinV0Radius", 16.0, "min v0 radius"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMaxV0Radius{"cfgMaxV0Radius", 90.0, "max v0 radius"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMaxAlphaAp{"cfgMaxAlphaAp", 0.95, "max alpha for AP cut"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMaxQtAp{"cfgMaxQtAp", 0.01, "max qT for AP cut"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMinCospa{"cfgMinCospa", 0.997, "min V0 CosPA"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMaxPca{"cfgMaxPca", 3.0, "max distance btween 2 legs"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMaxChi2Kf{"cfgMaxChi2Kf", 1e+10, "max chi2/ndf with KF"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgRejectV0OnItsib{"cfgRejectV0OnItsib", true, "flag to reject V0s on ITSib"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgDisableItsonlyTrack{"cfgDisableItsonlyTrack", false, "flag to disable ITSonly tracks"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgDisableTpconlyTrack{"cfgDisableTpconlyTrack", false, "flag to disable TPConly tracks"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMinNclusterTpc{"cfgMinNclusterTpc", 0, "min ncluster tpc"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMinNcrossedrows{"cfgMinNcrossedrows", 40, "min ncrossed rows"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMaxFracSharedClustersTpc{"cfgMaxFracSharedClustersTpc", 999.f, "max fraction of shared clusters in TPC"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMaxChi2Tpc{"cfgMaxChi2Tpc", 4.0, "max chi2/NclsTPC"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMaxChi2Its{"cfgMaxChi2Its", 36.0, "max chi2/NclsITS"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMinTpcNsigmaEl{"cfgMinTpcNsigmaEl", -3.0, "min. TPC n sigma for electron"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - Configurable cfgMaxTpcNsigmaEl{"cfgMaxTpcNsigmaEl", +3.0, "max. TPC n sigma for electron"}; // o2-linter: disable=name/configurable (upstream PWGEM name; JSON string cannot be changed without breaking existing configurations) - } pcmcuts; - - ~PhotonHBT() - { - delete emh1; - emh1 = 0x0; - delete emh2; - emh2 = 0x0; - map_mixed_eventId_to_globalBC.clear(); - used_photonIds_per_col.clear(); - used_photonIds_per_col.shrink_to_fit(); - } - - HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; - static constexpr std::string_view event_pair_types[2] = {"same/", "mix/"}; - - std::mt19937 engine; - std::uniform_int_distribution dist01; - int mRunNumber; - - std::vector zvtx_bin_edges; - std::vector cent_bin_edges; - std::vector ep2_bin_edges; - std::vector ep3_bin_edges; - std::vector occ_bin_edges; - - static constexpr float MinFloatTol = 1e-9f; - static constexpr float CosineMin = -1.f; - static constexpr float CosineMax = 1.f; - static constexpr float HalfAngle = 2.f; - static constexpr float PairEtaWeight = 0.5f; - - inline bool isInsideEllipse(float deta, float dphi) const - { - if (!ggpaircuts.cfgApplyEllipseCut.value) - return false; - const float sE = ggpaircuts.cfgEllipseSigEta.value; - const float sP = ggpaircuts.cfgEllipseSigPhi.value; - if (sE < MinFloatTol || sP < MinFloatTol) - return false; - return (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP) < ggpaircuts.cfgEllipseR2.value; - } - - static inline float normDphi(float dphi) - { - while (dphi > o2::constants::math::PI) - dphi -= o2::constants::math::TwoPI; // o2-linter: disable=two-pi-add-subtract (PWGEM has no PWGHF dependency; manual wrap is established pattern) - while (dphi < -o2::constants::math::PI) - dphi += o2::constants::math::TwoPI; // o2-linter: disable=two-pi-add-subtract (PWGEM has no PWGHF dependency; manual wrap is established pattern) - return dphi; - } - - static inline float computeCosTheta(const ROOT::Math::PtEtaPhiMVector& v1, - const ROOT::Math::PtEtaPhiMVector& v2) - { - ROOT::Math::PxPyPzEVector p1(v1), p2(v2); - ROOT::Math::PxPyPzEVector pair = p1 + p2; - ROOT::Math::Boost boost(-pair.BoostToCM()); - ROOT::Math::PxPyPzEVector p1cm = boost(p1); - ROOT::Math::XYZVector pairDir(pair.Px(), pair.Py(), pair.Pz()); - ROOT::Math::XYZVector p1cmDir(p1cm.Px(), p1cm.Py(), p1cm.Pz()); - if (pairDir.R() < MinFloatTol || p1cmDir.R() < MinFloatTol) - return -1.f; - return static_cast(pairDir.Unit().Dot(p1cmDir.Unit())); - } - - void init(InitContext& /*context*/) - { - mRunNumber = 0; - - auto parseBinsVerbatim = [](const ConfigurableAxis& cfg, std::vector& edges) { - if (cfg.value[0] == VARIABLE_WIDTH) { - edges = std::vector(cfg.value.begin(), cfg.value.end()); - edges.erase(edges.begin()); - } else { - int nbins = static_cast(cfg.value[0]); - float xmin = static_cast(cfg.value[1]); - float xmax = static_cast(cfg.value[2]); - edges.resize(nbins + 1); - for (int i = 0; i < nbins + 1; i++) - edges[i] = (xmax - xmin) / nbins * i + xmin; - } - }; - - parseBinsVerbatim(confVtxBins, zvtx_bin_edges); - parseBinsVerbatim(confCentBins, cent_bin_edges); - parseBinsVerbatim(confEpBins, ep2_bin_edges); - parseBinsVerbatim(confEp3Bins, ep3_bin_edges); - LOGF(info, "cfgOccupancyEstimator = %d", cfgOccupancyEstimator.value); - parseBinsVerbatim(confOccupancyBins, occ_bin_edges); - - emh1 = new MyEMH(ndepth); - emh2 = new MyEMH(ndepth); - - o2::aod::pwgem::photonmeson::utils::eventhistogram::addEventHistograms(&fRegistry); - DefineEMEventCut(); - DefinePCMCut(); - addhistograms(); - - std::random_device seed_gen; - engine = std::mt19937(seed_gen()); - dist01 = std::uniform_int_distribution(0, 1); - - fRegistry.add("Pair/mix/hDiffBC", - "diff. global BC in mixed event;|BC_{current} - BC_{mixed}|", - kTH1D, {{10001, -0.5, 10000.5}}, true); - } - - template - void initCCDB(TCollision const& collision) - { - if (mRunNumber == collision.runNumber()) - return; - mRunNumber = collision.runNumber(); - } - - void addhistograms() - { - static constexpr std::string_view qvec_det_names[6] = {"FT0M", "FT0A", "FT0C", "BTot", "BPos", "BNeg"}; - fRegistry.add("Event/before/hEP2_CentFT0C_forMix", - Form("2nd harmonics event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[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 event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[cfgEp2EstimatorForMix].data()), - kTH2F, {{110, 0, 110}, {180, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}}, false); - - const AxisSpec axis_kt{confKtBins, "k_{T} (GeV/c)"}; - const AxisSpec axis_qinv{confQBins, "q_{inv} (GeV/c)"}; - const AxisSpec axis_qabs_lcms{confQBins, "|#bf{q}|^{LCMS} (GeV/c)"}; - const AxisSpec axis_qout{confQBins, "q_{out} (GeV/c)"}; - const AxisSpec axis_qside{confQBins, "q_{side} (GeV/c)"}; - const AxisSpec axis_qlong{confQBins, "q_{long} (GeV/c)"}; - - const AxisSpec axisPt{confPtBins, "p_{T} (GeV/c)"}; - const AxisSpec axisEta{confEtaBins, "#eta"}; - const AxisSpec axisPhi{confPhiBins, "#phi (rad)"}; - const AxisSpec axisDeltaEta{confDeltaEtaBins, "#Delta#eta"}; - const AxisSpec axisDeltaPhi{confDeltaPhiBins, "#Delta#phi (rad)"}; - const AxisSpec axisEllipseVal{confEllipseValBins, "Ellipse val"}; - const AxisSpec axisCosTheta{confCosThetaBins, "cos(#theta*)"}; - const AxisSpec axisOpeningAngle{confOpeningAngleBins, "#alpha (rad)"}; - - const AxisSpec axisV0Combo{4, -0.5f, 3.5f, "V0 combo (0=Incl,1=II,2=IT,3=TT)"}; - const AxisSpec axisPairCombo{7, -0.5f, 6.5f, "Pair combo (0=Incl,1=II-II,2=II-IT,3=II-TT,4=IT-IT,5=IT-TT,6=TT-TT)"}; - const AxisSpec axisRxy{confRxyBins, "R_{xy} (cm)"}; - const AxisSpec axisPsiPair{confPsiPairBins, "#psi_{pair} (rad)"}; - - // ── Single-photon QA ───────────────────────────────────────────────────── - - fRegistry.add("SinglePhoton/hPt", "V0 photon p_{T};p_{T} (GeV/c);V0 combo", kTH2F, {axisPt, axisV0Combo}, true); - fRegistry.add("SinglePhoton/hEta", "V0 photon #eta;#eta;V0 combo", kTH2F, {axisEta, axisV0Combo}, true); - fRegistry.add("SinglePhoton/hPhi", "V0 photon #phi;#phi (rad);V0 combo", kTH2F, {axisPhi, axisV0Combo}, true); - fRegistry.add("SinglePhoton/hEtaVsPhi", "V0 photon acceptance;#phi (rad);#eta;V0 combo", kTH3F, {axisPhi, axisEta, axisV0Combo}, true); - fRegistry.add("SinglePhoton/hRxy", "Conversion R_{xy};R_{xy} (cm);V0 combo", kTH2F, {axisRxy, axisV0Combo}, true); - fRegistry.add("SinglePhoton/hPsiPair", "#psi_{pair};#psi_{pair} (rad);V0 combo", kTH2F, {axisPsiPair, axisV0Combo}, true); - - if (cfgDo3D) { - fRegistry.add("Pair/same/hs_3d", "diphoton correlation 3D LCMS", kTHnSparseD, {axis_qout, axis_qside, axis_qlong, axis_kt}, true); - } else { - if (cfgUseLcms) - fRegistry.add("Pair/same/hs_1d", "diphoton correlation 1D LCMS", kTHnSparseD, {axis_qabs_lcms, axis_kt}, true); - else - fRegistry.add("Pair/same/hs_1d", "diphoton correlation 1D", kTHnSparseD, {axis_qinv, axis_kt}, true); - } - if constexpr (pairtype == ggHBTPairType::kPCMPCM) - 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/hKt", "k_{T};k_{T} (GeV/c);counts", kTH1F, {axis_kt}, true); - - for (const auto& step : {"Before", "After"}) { - const std::string s = std::string("Pair/same/QA/") + step + "/"; - - fRegistry.add((s + "hDeltaEta").c_str(), "#Delta#eta;#Delta#eta;Pair combo", kTH2F, {axisDeltaEta, axisPairCombo}, true); - fRegistry.add((s + "hDeltaPhi").c_str(), "#Delta#phi;#Delta#phi (rad);Pair combo", kTH2F, {axisDeltaPhi, axisPairCombo}, true); - fRegistry.add((s + "hDEtaDPhi").c_str(), "#Delta#eta vs #Delta#phi;#Delta#eta;#Delta#phi (rad);Pair combo", kTH3F, {axisDeltaEta, axisDeltaPhi, axisPairCombo}, true); - fRegistry.add((s + "hDeltaEtaVsPairEta").c_str(), "#Delta#eta vs #LT#eta#GT_{pair};#LT#eta#GT_{pair};#Delta#eta;Pair combo", kTH3F, {axisEta, axisDeltaEta, axisPairCombo}, true); - fRegistry.add((s + "hCosTheta").c_str(), "cos(#theta*);cos(#theta*);Pair combo", kTH2F, {axisCosTheta, axisPairCombo}, true); - fRegistry.add((s + "hOpeningAngle").c_str(), "Opening angle;#alpha (rad);Pair combo", kTH2F, {axisOpeningAngle, axisPairCombo}, true); - fRegistry.add((s + "hEllipseVal").c_str(), "Ellipse value;value;Pair combo", kTH2F, {axisEllipseVal, axisPairCombo}, true); - } - - fRegistry.addClone("Pair/same/", "Pair/mix/"); - } - - void DefineEMEventCut() - { - fEMEventCut = EMPhotonEventCut("fEMEventCut", "fEMEventCut"); - fEMEventCut.SetRequireSel8(eventcuts.cfgRequireSel8); - fEMEventCut.SetRequireFT0AND(eventcuts.cfgRequireFT0AND); - fEMEventCut.SetZvtxRange(eventcuts.cfgZvtxMin, eventcuts.cfgZvtxMax); - fEMEventCut.SetRequireNoTFB(eventcuts.cfgRequireNoTFB); - fEMEventCut.SetRequireNoITSROFB(eventcuts.cfgRequireNoITSROFB); - fEMEventCut.SetRequireNoSameBunchPileup(eventcuts.cfgRequireNoSameBunchPileup); - fEMEventCut.SetRequireVertexITSTPC(eventcuts.cfgRequireVertexITSTPC); - fEMEventCut.SetRequireGoodZvtxFT0vsPV(eventcuts.cfgRequireGoodZvtxFT0vsPV); - fEMEventCut.SetRequireNoCollInTimeRangeStandard(eventcuts.cfgRequireNoCollInTimeRangeStandard); - fEMEventCut.SetRequireNoCollInTimeRangeStrict(eventcuts.cfgRequireNoCollInTimeRangeStrict); - fEMEventCut.SetRequireNoCollInITSROFStandard(eventcuts.cfgRequireNoCollInITSROFStandard); - fEMEventCut.SetRequireNoCollInITSROFStrict(eventcuts.cfgRequireNoCollInITSROFStrict); - fEMEventCut.SetRequireNoHighMultCollInPrevRof(eventcuts.cfgRequireNoHighMultCollInPrevRof); - fEMEventCut.SetRequireGoodITSLayer3(eventcuts.cfgRequireGoodITSLayer3); - fEMEventCut.SetRequireGoodITSLayer0123(eventcuts.cfgRequireGoodITSLayer0123); - fEMEventCut.SetRequireGoodITSLayersAll(eventcuts.cfgRequireGoodITSLayersAll); - } - - void DefinePCMCut() - { - fV0PhotonCut = V0PhotonCut("fV0PhotonCut", "fV0PhotonCut"); - fV0PhotonCut.SetV0PtRange(pcmcuts.cfgMinPtV0, 1e10f); - fV0PhotonCut.SetV0EtaRange(-pcmcuts.cfgMaxEtaV0, +pcmcuts.cfgMaxEtaV0); - fV0PhotonCut.SetMinCosPA(pcmcuts.cfgMinCospa); - fV0PhotonCut.SetMaxPCA(pcmcuts.cfgMaxPca); - fV0PhotonCut.SetMaxChi2KF(pcmcuts.cfgMaxChi2Kf); - fV0PhotonCut.SetRxyRange(pcmcuts.cfgMinV0Radius, pcmcuts.cfgMaxV0Radius); - fV0PhotonCut.SetAPRange(pcmcuts.cfgMaxAlphaAp, pcmcuts.cfgMaxQtAp); - fV0PhotonCut.RejectITSib(pcmcuts.cfgRejectV0OnItsib); - fV0PhotonCut.SetMinNClustersTPC(pcmcuts.cfgMinNclusterTpc); - fV0PhotonCut.SetMinNCrossedRowsTPC(pcmcuts.cfgMinNcrossedrows); - fV0PhotonCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); - fV0PhotonCut.SetMaxFracSharedClustersTPC(pcmcuts.cfgMaxFracSharedClustersTpc); - fV0PhotonCut.SetChi2PerClusterTPC(0.0, pcmcuts.cfgMaxChi2Tpc); - fV0PhotonCut.SetTPCNsigmaElRange(pcmcuts.cfgMinTpcNsigmaEl, pcmcuts.cfgMaxTpcNsigmaEl); - fV0PhotonCut.SetChi2PerClusterITS(-1e+10, pcmcuts.cfgMaxChi2Its); - fV0PhotonCut.SetDisableITSonly(pcmcuts.cfgDisableItsonlyTrack); - fV0PhotonCut.SetDisableTPConly(pcmcuts.cfgDisableTpconlyTrack); - fV0PhotonCut.SetNClustersITS(0, 7); - fV0PhotonCut.SetMeanClusterSizeITSob(0.0, 16.0); - fV0PhotonCut.SetRequireITSTPC(pcmcuts.cfgRequireV0WithItstpc); - fV0PhotonCut.SetRequireITSonly(pcmcuts.cfgRequireV0WithItsonly); - fV0PhotonCut.SetRequireTPConly(pcmcuts.cfgRequireV0WithTpconly); - } - - template - void fillPairHistogram(TCollision const&, - const ROOT::Math::PtEtaPhiMVector v1, - const ROOT::Math::PtEtaPhiMVector v2, - const float weight = 1.f) - { - float rndm = std::pow(-1, dist01(engine) % 2); - ROOT::Math::PtEtaPhiMVector q12 = (v1 - v2) * rndm; - ROOT::Math::PtEtaPhiMVector k12 = 0.5 * (v1 + v2); - float qinv = -q12.M(); - float kt = k12.Pt(); - 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); - ROOT::Math::PxPyPzEVector q12c = (v1c - v2c) * rndm; - float beta_z = (v1 + v2).Beta() * std::cos((v1 + v2).Theta()); - ROOT::Math::Boost bst_z(0, 0, -beta_z); - ROOT::Math::PxPyPzEVector q12_lcms = bst_z(q12c); - ROOT::Math::XYZVector q_3d_lcms = q12_lcms.Vect(); - float qabs_lcms = q_3d_lcms.R(); - float qout_lcms = q_3d_lcms.Dot(uv_out); - float qside_lcms = q_3d_lcms.Dot(uv_side); - float qlong_lcms = q_3d_lcms.Dot(uv_long); - if (cfgDo3D) - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_3d"), - std::fabs(qout_lcms), std::fabs(qside_lcms), std::fabs(qlong_lcms), kt, weight); - else { - if (cfgUseLcms) - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_1d"), qabs_lcms, kt, weight); - else - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_1d"), qinv, kt, weight); - } - - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hKt"), kt, weight); - } - - template - inline void fillPairQAStep(float deta, float dphi, float pairEta, - float cosTheta, float openingAngle, float ellVal, - int pairComboIdx) - { - if (!doPairQa) - return; - - if constexpr (IsBefore) { - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDeltaEta"), deta, 0.f); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDeltaPhi"), dphi, 0.f); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDEtaDPhi"), deta, dphi, 0.f); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDeltaEtaVsPairEta"), pairEta, deta, 0.f); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hCosTheta"), cosTheta, 0.f); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hOpeningAngle"), openingAngle, 0.f); - if (ellVal >= 0.f) - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hEllipseVal"), ellVal, 0.f); - if (pairComboIdx > 0) { - const float pcb = float(pairComboIdx); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDeltaEta"), deta, pcb); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDeltaPhi"), dphi, pcb); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDEtaDPhi"), deta, dphi, pcb); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hDeltaEtaVsPairEta"), pairEta, deta, pcb); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hCosTheta"), cosTheta, pcb); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hOpeningAngle"), openingAngle, pcb); - if (ellVal >= 0.f) - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/Before/hEllipseVal"), ellVal, pcb); - } - } else { - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDeltaEta"), deta, 0.f); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDeltaPhi"), dphi, 0.f); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDEtaDPhi"), deta, dphi, 0.f); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDeltaEtaVsPairEta"), pairEta, deta, 0.f); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hCosTheta"), cosTheta, 0.f); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hOpeningAngle"), openingAngle, 0.f); - if (ellVal >= 0.f) - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hEllipseVal"), ellVal, 0.f); - if (pairComboIdx > 0) { - const float pcb = float(pairComboIdx); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDeltaEta"), deta, pcb); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDeltaPhi"), dphi, pcb); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDEtaDPhi"), deta, dphi, pcb); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hDeltaEtaVsPairEta"), pairEta, deta, pcb); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hCosTheta"), cosTheta, pcb); - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hOpeningAngle"), openingAngle, pcb); - if (ellVal >= 0.f) - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("QA/After/hEllipseVal"), ellVal, pcb); - } - } - } - - // =========================================================================== - // runPairing - // =========================================================================== - - template - void runPairing(TCollisions const& collisions, - TPhotons1 const& photons1, TPhotons2 const& photons2, - TSubInfos1 const&, TSubInfos2 const&, - TPreslice1 const& perCollision1, TPreslice2 const& perCollision2, - TCut1 const& cut1, TCut2 const& cut2) - { - for (const auto& collision : collisions) { - initCCDB(collision); - int ndiphoton = 0; - - const float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; - if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) - continue; - - const float ep2_arr[6] = {collision.ep2ft0m(), collision.ep2ft0a(), collision.ep2ft0c(), - collision.ep2btot(), collision.ep2bpos(), collision.ep2bneg()}; - const float ep3_arr[6] = {collision.ep3ft0m(), collision.ep3ft0a(), collision.ep3ft0c(), - collision.ep3btot(), collision.ep3bpos(), collision.ep3bneg()}; - float ep2 = ep2_arr[cfgEp2EstimatorForMix]; - float ep3 = ep3_arr[cfgEp3EstimatorForMix]; - - 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/after/hEP2_CentFT0C_forMix"), collision.centFT0C(), ep2); - - auto clampBin = [](int b, int nmax) { - return (b < 0) ? 0 : (b > nmax ? nmax : b); - }; - auto binOf = [&](const std::vector& edges, float val) { - int b = static_cast(std::lower_bound(edges.begin(), edges.end(), val) - edges.begin()) - 1; - return clampBin(b, static_cast(edges.size()) - 2); - }; - - int zbin = binOf(zvtx_bin_edges, collision.posZ()); - int centbin = binOf(cent_bin_edges, centralities[cfgCentEstimator]); - int ep2bin = binOf(ep2_bin_edges, ep2); - int ep3bin = binOf(ep3_bin_edges, ep3); - int occbin = binOf(occ_bin_edges, - cfgOccupancyEstimator == 1 - ? static_cast(collision.trackOccupancyInTimeRange()) - : collision.ft0cOccupancyInTimeRange()); - - auto key_bin = std::make_tuple(zbin, centbin, ep2bin, ep3bin, occbin); - auto key_df_collision = std::make_pair(ndf, collision.globalIndex()); - - if constexpr (pairtype == ggHBTPairType::kPCMPCM) { - - auto photons1_coll = photons1.sliceBy(perCollision1, collision.globalIndex()); - auto photons2_coll = photons2.sliceBy(perCollision2, collision.globalIndex()); - - // ── Single-photon QA ───────────────────────────────────────────────── - if (doSinglePhotonQa) { - for (const auto& g : photons1_coll) { - if (!cut1.template IsSelected(g)) - continue; - const int ci = classifyV0ComboIdx(g); // 0=Other/1/2/3 - fRegistry.fill(HIST("SinglePhoton/hPt"), g.pt(), 0.f); - fRegistry.fill(HIST("SinglePhoton/hEta"), g.eta(), 0.f); - fRegistry.fill(HIST("SinglePhoton/hPhi"), g.phi(), 0.f); - fRegistry.fill(HIST("SinglePhoton/hEtaVsPhi"), g.phi(), g.eta(), 0.f); - fRegistry.fill(HIST("SinglePhoton/hRxy"), std::hypot(g.vx(), g.vy()), 0.f); - if constexpr (requires { g.psipair(); }) { - fRegistry.fill(HIST("SinglePhoton/hPsiPair"), g.psipair(), 0.f); - } - if (ci > 0) { - const float cb = float(ci); - fRegistry.fill(HIST("SinglePhoton/hPt"), g.pt(), cb); - fRegistry.fill(HIST("SinglePhoton/hEta"), g.eta(), cb); - fRegistry.fill(HIST("SinglePhoton/hPhi"), g.phi(), cb); - fRegistry.fill(HIST("SinglePhoton/hEtaVsPhi"), g.phi(), g.eta(), cb); - fRegistry.fill(HIST("SinglePhoton/hRxy"), std::hypot(g.vx(), g.vy()), cb); - if constexpr (requires { g.psipair(); }) { - fRegistry.fill(HIST("SinglePhoton/hPsiPair"), g.psipair(), cb); - } - } - } - } - - // ── Same-event pair loop ────────────────────────────────────────────── - for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photons1_coll, photons2_coll))) { - 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()) - continue; - - const int c1 = classifyV0ComboIdx(g1); - const int c2 = classifyV0ComboIdx(g2); - const int pcb = pairComboBin(c1, c2); // 0 if either is Other - - ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); - ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); - - 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)); - ROOT::Math::XYZVector cp1(g1.vx(), g1.vy(), g1.vz()); - ROOT::Math::XYZVector cp2(g2.vx(), g2.vy(), g2.vz()); - float opa = std::acos(std::clamp(static_cast(cp1.Dot(cp2) / (std::sqrt(cp1.Mag2()) * std::sqrt(cp2.Mag2()))), CosineMin, CosineMax)); - o2::math_utils::bringTo02Pi(opa); - if (opa > o2::constants::math::PI) - opa -= o2::constants::math::PI; - float cosOA = std::cos(opa / HalfAngle); - if (dr / cosOA < ggpaircuts.cfgMinDrCosOa) - continue; - fRegistry.fill(HIST("Pair/same/hDeltaRCosOA"), dr / cosOA); - - float deta = g1.eta() - g2.eta(); - float dphi = normDphi(g1.phi() - g2.phi()); - float pairEta = PairEtaWeight * (g1.eta() + g2.eta()); - float cosTheta = std::fabs(computeCosTheta(v1, v2)); - const float sE = ggpaircuts.cfgEllipseSigEta.value; - const float sP = ggpaircuts.cfgEllipseSigPhi.value; - float ellVal = (sE > MinFloatTol && sP > MinFloatTol) - ? (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP) - : -1.f; - - // QA Before ellipse cut - fillPairQAStep<0, true>(deta, dphi, pairEta, cosTheta, opa, ellVal, pcb); - - // Ellipse cut - if (isInsideEllipse(deta, dphi)) - continue; - - // QA After ellipse cut - fillPairQAStep<0, false>(deta, dphi, pairEta, cosTheta, opa, ellVal, pcb); - - fillPairHistogram<0>(collision, v1, v2, 1.f); - ndiphoton++; - - if (std::find(used_photonIds_per_col.begin(), used_photonIds_per_col.end(), g1.globalIndex()) == used_photonIds_per_col.end()) { - EMPair g1tmp = EMPair(g1.pt(), g1.eta(), g1.phi(), 0); - g1tmp.setConversionPointXYZ(g1.vx(), g1.vy(), g1.vz()); - emh1->AddTrackToEventPool(key_df_collision, g1tmp); - used_photonIds_per_col.emplace_back(g1.globalIndex()); - } - if (std::find(used_photonIds_per_col.begin(), used_photonIds_per_col.end(), g2.globalIndex()) == used_photonIds_per_col.end()) { - EMPair g2tmp = EMPair(g2.pt(), g2.eta(), g2.phi(), 0); - g2tmp.setConversionPointXYZ(g2.vx(), g2.vy(), g2.vz()); - emh1->AddTrackToEventPool(key_df_collision, g2tmp); - used_photonIds_per_col.emplace_back(g2.globalIndex()); - } - } // end same-event pair loop - } - - used_photonIds_per_col.clear(); - used_photonIds_per_col.shrink_to_fit(); - - // ── Mixed-event loop ──────────────────────────────────────────────────── - if (!cfgDoMix || !(ndiphoton > 0)) - continue; - - auto selected_photons1_in_this_event = emh1->GetTracksPerCollision(key_df_collision); - auto collisionIds1_in_mixing_pool = emh1->GetCollisionIdsFromEventPool(key_bin); - - if constexpr (pairtype == ggHBTPairType::kPCMPCM) { - for (const auto& mix_dfId_collisionId : collisionIds1_in_mixing_pool) { - int mix_dfId = mix_dfId_collisionId.first; - int64_t mix_collisionId = mix_dfId_collisionId.second; - if (collision.globalIndex() == mix_collisionId && ndf == mix_dfId) - continue; - - auto globalBC_mix = map_mixed_eventId_to_globalBC[mix_dfId_collisionId]; - uint64_t diffBC = std::max(collision.globalBC(), globalBC_mix) - std::min(collision.globalBC(), globalBC_mix); - fRegistry.fill(HIST("Pair/mix/hDiffBC"), diffBC); - if (diffBC < ndiffBcMix) - continue; - - auto photons1_from_event_pool = emh1->GetTracksPerCollision(mix_dfId_collisionId); - - for (const auto& g1 : selected_photons1_in_this_event) { - for (const auto& g2 : photons1_from_event_pool) { - ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); - ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); - - 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)); - ROOT::Math::XYZVector cp1(g1.vx(), g1.vy(), g1.vz()); - ROOT::Math::XYZVector cp2(g2.vx(), g2.vy(), g2.vz()); - float opa = std::acos(std::clamp(static_cast(cp1.Dot(cp2) / (std::sqrt(cp1.Mag2()) * std::sqrt(cp2.Mag2()))), CosineMin, CosineMax)); - o2::math_utils::bringTo02Pi(opa); - if (opa > o2::constants::math::PI) - opa -= o2::constants::math::PI; - float cosOA = std::cos(opa / HalfAngle); - if (dr / cosOA < ggpaircuts.cfgMinDrCosOa) - continue; - fRegistry.fill(HIST("Pair/mix/hDeltaRCosOA"), dr / cosOA); - - float deta = g1.eta() - g2.eta(); - float dphi = normDphi(g1.phi() - g2.phi()); - float pairEta = PairEtaWeight * (g1.eta() + g2.eta()); - float cosTheta = std::fabs(computeCosTheta(v1, v2)); - const float sE = ggpaircuts.cfgEllipseSigEta.value; - const float sP = ggpaircuts.cfgEllipseSigPhi.value; - float ellVal = (sE > MinFloatTol && sP > MinFloatTol) - ? (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP) - : -1.f; - - fillPairQAStep<1, true>(deta, dphi, pairEta, cosTheta, opa, ellVal, 0); - - if (isInsideEllipse(deta, dphi)) - continue; - - fillPairQAStep<1, false>(deta, dphi, pairEta, cosTheta, opa, ellVal, 0); - - fillPairHistogram<1>(collision, v1, v2, 1.f); - } - } - } // end mixed event pool loop - } - - if (ndiphoton > 0) { - emh1->AddCollisionIdAtLast(key_bin, key_df_collision); - emh2->AddCollisionIdAtLast(key_bin, key_df_collision); - map_mixed_eventId_to_globalBC[key_df_collision] = collision.globalBC(); - } - } // end collision loop - } - - using MyEMH = o2::aod::pwgem::dilepton::utils::EventMixingHandler< - std::tuple, - std::pair, EMPair>; - MyEMH* emh1 = nullptr; - MyEMH* emh2 = nullptr; - std::vector used_photonIds_per_col; - std::map, uint64_t> map_mixed_eventId_to_globalBC; - - SliceCache cache; - Preslice perCollision_pcm = aod::v0photonkf::pmeventId; - - Filter collisionFilter_centrality = (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 collisionFilter_occupancy_track = eventcuts.cfgTrackOccupancyMin <= o2::aod::evsel::trackOccupancyInTimeRange && - o2::aod::evsel::trackOccupancyInTimeRange < eventcuts.cfgTrackOccupancyMax; - Filter collisionFilter_occupancy_ft0c = eventcuts.cfgFT0COccupancyMin <= o2::aod::evsel::ft0cOccupancyInTimeRange && - o2::aod::evsel::ft0cOccupancyInTimeRange < eventcuts.cfgFT0COccupancyMax; - using FilteredMyCollisions = soa::Filtered; - - int ndf = 0; - void processAnalysis(FilteredMyCollisions const& collisions, Types const&... args) - { - if constexpr (pairtype == ggHBTPairType::kPCMPCM) { - auto v0photons = std::get<0>(std::tie(args...)); - auto v0legs = std::get<1>(std::tie(args...)); - runPairing(collisions, v0photons, v0photons, v0legs, v0legs, - perCollision_pcm, perCollision_pcm, fV0PhotonCut, fV0PhotonCut); - } - ndf++; - } - PROCESS_SWITCH(PhotonHBT, processAnalysis, "pairing for analysis", false); - - void processDummy(MyCollisions const&) {} - PROCESS_SWITCH(PhotonHBT, processDummy, "Dummy function", true); -}; - -#endif // PWGEM_PHOTONMESON_CORE_PHOTONHBT_H_ diff --git a/PWGEM/PhotonMeson/Tasks/CMakeLists.txt b/PWGEM/PhotonMeson/Tasks/CMakeLists.txt index e0ce318628f..fe155e91fe1 100644 --- a/PWGEM/PhotonMeson/Tasks/CMakeLists.txt +++ b/PWGEM/PhotonMeson/Tasks/CMakeLists.txt @@ -131,11 +131,6 @@ o2physics_add_dpl_workflow(pi0eta-to-gammagamma-mc-emcemc PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(photon-hbt-pcmpcm - SOURCES PhotonHBTPCMPCM.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore - COMPONENT_NAME Analysis) - o2physics_add_dpl_workflow(tag-and-probe SOURCES TagAndProbe.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore O2Physics::MLCore @@ -200,3 +195,8 @@ o2physics_add_dpl_workflow(photon-reso-task SOURCES photonResoTask.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::EMCALBase O2::EMCALCalib O2Physics::PWGEMPhotonMesonCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(photonhbt + SOURCES photonhbt.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore + COMPONENT_NAME Analysis) diff --git a/PWGEM/PhotonMeson/Tasks/PhotonHBTPCMPCM.cxx b/PWGEM/PhotonMeson/Tasks/PhotonHBTPCMPCM.cxx deleted file mode 100644 index 0fde469c55c..00000000000 --- a/PWGEM/PhotonMeson/Tasks/PhotonHBTPCMPCM.cxx +++ /dev/null @@ -1,27 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. -// -// ======================== -// -// This code loops over v0 photons and makes pairs for photon HBT analysis. -// Please write to: daiki.sekihata@cern.ch - -#include "PWGEM/PhotonMeson/Core/PhotonHBT.h" - -#include "Framework/AnalysisDataModel.h" -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask>(cfgc, TaskName{"photon-hbt-pcmpcm"})}; -} diff --git a/PWGEM/PhotonMeson/Tasks/photonhbt.cxx b/PWGEM/PhotonMeson/Tasks/photonhbt.cxx new file mode 100644 index 00000000000..bacb008cd47 --- /dev/null +++ b/PWGEM/PhotonMeson/Tasks/photonhbt.cxx @@ -0,0 +1,784 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// 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 + +#include "PWGEM/Dilepton/Utils/EMTrack.h" +#include "PWGEM/Dilepton/Utils/EventMixingHandler.h" +#include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" +#include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" +#include "PWGEM/PhotonMeson/DataModel/EventTables.h" +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "PWGEM/PhotonMeson/Utils/EventHistograms.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/// Single-photon track-type combo. +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 +}; + +/// Photon-pair track-type combo. +enum class PairCombo : int { + Inclusive = 0, + IiXIi = 1, + IiXIt = 2, + IiXTt = 3, + ItXIt = 4, + ItXTt = 5, + TtXTt = 6, +}; + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; +using namespace o2::aod::pwgem::dilepton::utils; + +using MyCollisions = soa::Join; +using MyCollision = MyCollisions::iterator; + +using MyV0Photons = soa::Join; +using MyV0Photon = MyV0Photons::iterator; + +struct photonhbt { + + template + static inline V0Combo classifyV0Combo(TGamma const& g) + { + const auto pos = g.template posTrack_as(); + const auto neg = g.template negTrack_as(); + const bool posII = pos.hasITS() && pos.hasTPC(); + const bool posTPC = !pos.hasITS() && pos.hasTPC(); + const bool negII = neg.hasITS() && neg.hasTPC(); + const bool negTPC = !neg.hasITS() && neg.hasTPC(); + if (posII && negII) + return V0Combo::ItstpcItstpc; + if ((posII && negTPC) || (posTPC && negII)) + return V0Combo::ItstpcTpconly; + if (posTPC && negTPC) + return V0Combo::TpconlyTpconly; + return V0Combo::Inclusive; + } + + static inline PairCombo classifyPairCombo(V0Combo c1, V0Combo c2) + { + const int i1 = static_cast(c1); + const int i2 = static_cast(c2); + if (i1 <= 0 || i2 <= 0) + 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}}}; + return static_cast(kTable[lo][hi]); + } + + // --------------------------------------------------------------------------- + // Configurables: histogram axes + // --------------------------------------------------------------------------- + + // HBT physics + ConfigurableAxis confQBins{"confQBins", {60, 0, +0.3f}, "q bins for output histograms"}; + ConfigurableAxis confKtBins{"confKtBins", {VARIABLE_WIDTH, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6}, "kT bins"}; + + // 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)"}; + + // Pair QA + 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 + 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)"}; + const AxisSpec axisQout{confQBins, "q_{out} (GeV/c)"}; + const AxisSpec axisQside{confQBins, "q_{side} (GeV/c)"}; + const AxisSpec axisQlong{confQBins, "q_{long} (GeV/c)"}; + const AxisSpec axisPt{confPtBins, "p_{T} (GeV/c)"}; + const AxisSpec axisEta{confEtaBins, "#eta"}; + const AxisSpec axisPhi{confPhiBins, "#phi (rad)"}; + const AxisSpec axisDeltaEta{confDeltaEtaBins, "#Delta#eta"}; + const AxisSpec axisDeltaPhi{confDeltaPhiBins, "#Delta#phi (rad)"}; + 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)"}; + + // --------------------------------------------------------------------------- + // Configurables: QA flags + // --------------------------------------------------------------------------- + + struct : ConfigurableGroup { + std::string prefix = "qaflags_group"; + Configurable doPairQa{"doPairQa", true, "fill pair QA histograms (Before/After ellipse cut)"}; + Configurable doSinglePhotonQa{"doSinglePhotonQa", true, "fill single-photon QA histograms (pT, eta, phi)"}; + } qaflags; + + // --------------------------------------------------------------------------- + // Configurables: HBT kind + // --------------------------------------------------------------------------- + + Configurable cfgDo3D{"cfgDo3D", false, "enable 3D analysis"}; + Configurable cfgUseLCMS{"cfgUseLCMS", true, "measure relative momentum in LCMS for 1D"}; + + // --------------------------------------------------------------------------- + // Configurables: events + // --------------------------------------------------------------------------- + + Configurable cfgCentMin{"cfgCentMin", -1, "min. centrality"}; + Configurable cfgCentMax{"cfgCentMax", 999, "max. centrality"}; + Configurable maxY{"maxY", 0.9, "maximum rapidity"}; + + // --------------------------------------------------------------------------- + // 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 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"}; + + ConfigurableAxis ConfVtxBins{"ConfVtxBins", {VARIABLE_WIDTH, -10.f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; + ConfigurableAxis ConfCentBins{"ConfCentBins", {VARIABLE_WIDTH, 0.f, 5.f, 10.f, 20.f, 30.f, 40.f, 50.f, 60.f, 70.f, 80.f, 90.f, 100.f, 999.f}, "Mixing bins - centrality"}; + 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 + // --------------------------------------------------------------------------- + + 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"}; + 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"}; + } ggpaircuts; + + // --------------------------------------------------------------------------- + // Event cut + // --------------------------------------------------------------------------- + + EMPhotonEventCut fEMEventCut; + struct : ConfigurableGroup { + std::string prefix = "eventcut_group"; + Configurable cfgZvtxMin{"cfgZvtxMin", -10.f, "min. Zvtx"}; + Configurable cfgZvtxMax{"cfgZvtxMax", +10.f, "max. Zvtx"}; + Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8"}; + Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND"}; + Configurable cfgRequireNoTFB{"cfgRequireNoTFB", true, "require no TF border"}; + Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", true, "require no ITS ROF border"}; + Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup"}; + Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC"}; + Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx FT0 vs PV"}; + Configurable cfgTrackOccupancyMin{"cfgTrackOccupancyMin", -2, "min. track occupancy"}; + Configurable cfgTrackOccupancyMax{"cfgTrackOccupancyMax", 1000000000, "max. track occupancy"}; + Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -2.f, "min. FT0C occupancy"}; + Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000.f, "max. FT0C occupancy"}; + Configurable cfgRequireNoCollInTimeRangeStandard{"cfgRequireNoCollInTimeRangeStandard", false, "no coll in time range std"}; + Configurable cfgRequireNoCollInTimeRangeStrict{"cfgRequireNoCollInTimeRangeStrict", false, "no coll in time range strict"}; + Configurable cfgRequireNoCollInITSROFStandard{"cfgRequireNoCollInITSROFStandard", false, "no coll in ITS ROF std"}; + Configurable cfgRequireNoCollInITSROFStrict{"cfgRequireNoCollInITSROFStrict", false, "no coll in ITS ROF strict"}; + Configurable cfgRequireNoHighMultCollInPrevRof{"cfgRequireNoHighMultCollInPrevRof", false, "no HM coll in prev ROF"}; + Configurable cfgRequireGoodITSLayer3{"cfgRequireGoodITSLayer3", false, "ITS layer 3 OK"}; + Configurable cfgRequireGoodITSLayer0123{"cfgRequireGoodITSLayer0123", false, "ITS layers 0-3 OK"}; + Configurable cfgRequireGoodITSLayersAll{"cfgRequireGoodITSLayersAll", false, "all ITS layers OK"}; + } eventcuts; + + // --------------------------------------------------------------------------- + // 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; + + ~photonhbt() + { + delete emh1; + emh1 = nullptr; + delete emh2; + emh2 = nullptr; + mapMixedEventIdToGlobalBC.clear(); + usedPhotonIdsPerCol.clear(); + usedPhotonIdsPerCol.shrink_to_fit(); + } + + HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + static constexpr std::string_view event_pair_types[2] = {"same/", "mix/"}; + + std::mt19937 engine; + std::uniform_int_distribution dist01; + int mRunNumber; + + std::vector ztxBinEdges; + std::vector centBinEdges; + std::vector epBinEgdes; + std::vector occBinEdges; + + inline bool isInsideEllipse(float deta, float dphi) const + { + if (!ggpaircuts.cfgApplyEllipseCut.value) + return false; + const float sE = ggpaircuts.cfgEllipseSigEta.value; + const float sP = ggpaircuts.cfgEllipseSigPhi.value; + if (sE < 1e-9f || sP < 1e-9f) + return false; + return (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP) < ggpaircuts.cfgEllipseR2.value; + } + + static inline float computeCosTheta(const ROOT::Math::PtEtaPhiMVector& v1, + const ROOT::Math::PtEtaPhiMVector& v2) + { + ROOT::Math::PxPyPzEVector p1(v1), p2(v2); + ROOT::Math::PxPyPzEVector pair = p1 + p2; + ROOT::Math::Boost boost(-pair.BoostToCM()); + ROOT::Math::PxPyPzEVector p1cm = boost(p1); + ROOT::Math::XYZVector pairDir(pair.Px(), pair.Py(), pair.Pz()); + ROOT::Math::XYZVector p1cmDir(p1cm.Px(), p1cm.Py(), p1cm.Pz()); + if (pairDir.R() < 1e-9 || p1cmDir.R() < 1e-9) + return -1.f; + return static_cast(pairDir.Unit().Dot(p1cmDir.Unit())); + } + + static void parseBins(const ConfigurableAxis& cfg, std::vector& edges) + { + if (cfg.value[0] == VARIABLE_WIDTH) { + edges = std::vector(cfg.value.begin(), cfg.value.end()); + edges.erase(edges.begin()); + } else { + const int n = static_cast(cfg.value[0]); + const float xmin = static_cast(cfg.value[1]); + const float xmax = static_cast(cfg.value[2]); + edges.resize(n + 1); + for (int i = 0; i <= n; ++i) + edges[i] = xmin + (xmax - xmin) / n * i; + } + } + + /// 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; + return clampBin(b, static_cast(edges.size()) - 2); + } + + void init(InitContext& /*context*/) + { + mRunNumber = 0; + + parseBins(ConfVtxBins, ztxBinEdges); + parseBins(ConfCentBins, centBinEdges); + parseBins(ConfEPBins, epBinEgdes); + parseBins(ConfOccupancyBins, occBinEdges); + + emh1 = new MyEMH(ndepth); + emh2 = new MyEMH(ndepth); + + o2::aod::pwgem::photonmeson::utils::eventhistogram::addEventHistograms(&fRegistry); + DefineEMEventCut(); + DefinePCMCut(); + addhistograms(); + + std::random_device seedGen; + engine = std::mt19937(seedGen()); + dist01 = std::uniform_int_distribution(0, 1); + + fRegistry.add("Pair/mix/hDiffBC", + "diff. global BC in mixed event;|BC_{current} - BC_{mixed}|", + kTH1D, {{10001, -0.5, 10000.5}}, true); + } + + template + void initCCDB(TCollision const& collision) + { + if (mRunNumber == collision.runNumber()) + return; + mRunNumber = collision.runNumber(); + } + 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); + fRegistry.add((path + "hEllipseVal").c_str(), "(#Delta#eta/#sigma_{#eta})^{2}+(#Delta#phi/#sigma_{#phi})^{2};value;counts", kTH1D, {axisEllipseVal}, 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 ────────────────────────────────────────────────────────── + if (cfgDo3D) { + fRegistry.add("Pair/same/hs_3d", "diphoton correlation 3D LCMS", kTHnSparseD, {axisQout, axisQside, axisQlong, axisKt}, true); + } else { + if (cfgUseLCMS) { + fRegistry.add("Pair/same/hs_1d", "diphoton correlation 1D LCMS", kTHnSparseD, {axisQabsLcms, axisKt}, true); + } else { + fRegistry.add("Pair/same/hs_1d", "diphoton correlation 1D", kTHnSparseD, {axisQinv, axisKt}, true); + } + } + + fRegistry.add("Pair/same/hDeltaRCosOA", "distance between 2 conversion points;#Deltar/cos(#theta_{op}/2) (cm)", kTH1D, {{100, 0, 100}}, true); + + addQAHistogramsForStep("Pair/same/QA/Before/"); + addQAHistogramsForStep("Pair/same/QA/After/"); + + fRegistry.addClone("Pair/same/", "Pair/mix/"); + } + + // --------------------------------------------------------------------------- + // DefineEMEventCut / DefinePCMCut + // --------------------------------------------------------------------------- + + void DefineEMEventCut() + { + fEMEventCut = EMPhotonEventCut("fEMEventCut", "fEMEventCut"); + fEMEventCut.SetRequireSel8(eventcuts.cfgRequireSel8); + fEMEventCut.SetRequireFT0AND(eventcuts.cfgRequireFT0AND); + fEMEventCut.SetZvtxRange(eventcuts.cfgZvtxMin, eventcuts.cfgZvtxMax); + fEMEventCut.SetRequireNoTFB(eventcuts.cfgRequireNoTFB); + fEMEventCut.SetRequireNoITSROFB(eventcuts.cfgRequireNoITSROFB); + fEMEventCut.SetRequireNoSameBunchPileup(eventcuts.cfgRequireNoSameBunchPileup); + fEMEventCut.SetRequireVertexITSTPC(eventcuts.cfgRequireVertexITSTPC); + fEMEventCut.SetRequireGoodZvtxFT0vsPV(eventcuts.cfgRequireGoodZvtxFT0vsPV); + fEMEventCut.SetRequireNoCollInTimeRangeStandard(eventcuts.cfgRequireNoCollInTimeRangeStandard); + fEMEventCut.SetRequireNoCollInTimeRangeStrict(eventcuts.cfgRequireNoCollInTimeRangeStrict); + fEMEventCut.SetRequireNoCollInITSROFStandard(eventcuts.cfgRequireNoCollInITSROFStandard); + fEMEventCut.SetRequireNoCollInITSROFStrict(eventcuts.cfgRequireNoCollInITSROFStrict); + fEMEventCut.SetRequireNoHighMultCollInPrevRof(eventcuts.cfgRequireNoHighMultCollInPrevRof); + fEMEventCut.SetRequireGoodITSLayer3(eventcuts.cfgRequireGoodITSLayer3); + fEMEventCut.SetRequireGoodITSLayer0123(eventcuts.cfgRequireGoodITSLayer0123); + fEMEventCut.SetRequireGoodITSLayersAll(eventcuts.cfgRequireGoodITSLayersAll); + } + + void DefinePCMCut() + { + fV0PhotonCut = V0PhotonCut("fV0PhotonCut", "fV0PhotonCut"); + fV0PhotonCut.SetV0PtRange(pcmcuts.cfgMinPtV0, 1e10f); + fV0PhotonCut.SetV0EtaRange(-pcmcuts.cfgMaxEtaV0, +pcmcuts.cfgMaxEtaV0); + fV0PhotonCut.SetMinCosPA(pcmcuts.cfgMinCosPA); + fV0PhotonCut.SetMaxPCA(pcmcuts.cfgMaxPCA); + fV0PhotonCut.SetMaxChi2KF(pcmcuts.cfgMaxChi2KF); + fV0PhotonCut.SetRxyRange(pcmcuts.cfgMinV0Radius, pcmcuts.cfgMaxV0Radius); + fV0PhotonCut.SetAPRange(pcmcuts.cfgMaxAlphaAP, pcmcuts.cfgMaxQtAP); + fV0PhotonCut.RejectITSib(pcmcuts.cfgRejectV0OnITSIB); + fV0PhotonCut.SetMinNClustersTPC(pcmcuts.cfgMinNClusterTPC); + fV0PhotonCut.SetMinNCrossedRowsTPC(pcmcuts.cfgMinNCrossedRows); + fV0PhotonCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); + fV0PhotonCut.SetMaxFracSharedClustersTPC(pcmcuts.cfgMaxFracSharedClustersTPC); + fV0PhotonCut.SetChi2PerClusterTPC(0.0, pcmcuts.cfgMaxChi2TPC); + fV0PhotonCut.SetTPCNsigmaElRange(pcmcuts.cfgMinTPCNsigmaEl, pcmcuts.cfgMaxTPCNsigmaEl); + fV0PhotonCut.SetChi2PerClusterITS(-1e+10, pcmcuts.cfgMaxChi2ITS); + fV0PhotonCut.SetDisableITSonly(pcmcuts.cfgDisableITSOnlyTrack); + fV0PhotonCut.SetDisableTPConly(pcmcuts.cfgDisableTPCOnlyTrack); + fV0PhotonCut.SetNClustersITS(0, 7); + fV0PhotonCut.SetMeanClusterSizeITSob(0.0, 16.0); + fV0PhotonCut.SetRequireITSTPC(pcmcuts.cfgRequireV0WithITSTPC); + fV0PhotonCut.SetRequireITSonly(pcmcuts.cfgRequireV0WithITSOnly); + fV0PhotonCut.SetRequireTPConly(pcmcuts.cfgRequireV0WithTPCOnly); + } + + template + void fillPairHistogram(TCollision const&, + 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); + + if (cfgDo3D) { + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_3d"), + std::fabs(qout_lcms), std::fabs(qside_lcms), std::fabs(qlong_lcms), kt, weight); + } else { + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_1d"), + cfgUseLCMS ? qabs_lcms : qinv, kt, weight); + } + } + + template + inline void fillPairQAStep(float deta, float dphi, float pairEta, + float cosTheta, float openingAngle) + { + if (!qaflags.doPairQa) + return; + + const float sE = ggpaircuts.cfgEllipseSigEta.value; + const float sP = ggpaircuts.cfgEllipseSigPhi.value; + + 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)); + } + } + + template + void runPairing(TCollisions const& collisions, + TPhotons1 const& photons1, TPhotons2 const& photons2, + TSubInfos1 const&, TSubInfos2 const&, + TPreslice1 const& perCollision1, TPreslice2 const& perCollision2, + TCut1 const& cut1, TCut2 const& cut2) + { + 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()}; + 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/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()); + + auto keyBin = std::make_tuple(zbin, centbin, epbin, occbin); + auto keyDFCollision = std::make_pair(ndf, collision.globalIndex()); + + auto photons1Coll = photons1.sliceBy(perCollision1, collision.globalIndex()); + auto photons2Coll = photons2.sliceBy(perCollision2, collision.globalIndex()); + + // ── 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()); + } + } + + // ── 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()) + 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) + 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)) + continue; + + // ── QA: After ellipse cut ─────────────────────────────────────── + fillPairQAStep<0, false>(deta, dphi, pairEta, cosTheta, openingAngle); + + // ── Physics ───────────────────────────────────────────────── + fillPairHistogram<0>(collision, v1, 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); + 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 + } + + usedPhotonIdsPerCol.clear(); + usedPhotonIdsPerCol.shrink_to_fit(); + + // ── Mixed-event loop ──────────────────────────────────────────────────── + 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; + + uint64_t bcMix = mapMixedEventIdToGlobalBC[mixID]; + 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) { + + 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) + 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; + + // QA Before cut — mix/QA/Before/ histograms. + fillPairQAStep<1, true>(deta, dphi, pairEta, cosTheta, openingAngle); + + // Apply ellipse cut + if (isInsideEllipse(deta, dphi)) + continue; + + // QA After cut + fillPairQAStep<1, false>(deta, dphi, pairEta, cosTheta, openingAngle); + + fillPairHistogram<1>(collision, v1, v2, 1.f); + } + } + } + + if (ndiphoton > 0) { + emh1->AddCollisionIdAtLast(keyBin, keyDFCollision); + emh2->AddCollisionIdAtLast(keyBin, keyDFCollision); + mapMixedEventIdToGlobalBC[keyDFCollision] = collision.globalBC(); + } + } // end collision loop + } + + using MyEMH = o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, EMPair>; + MyEMH* emh1 = nullptr; + MyEMH* emh2 = nullptr; + std::vector 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; + using FilteredMyCollisions = soa::Filtered; + + int ndf = 0; + void processAnalysis(FilteredMyCollisions const& collisions, MyV0Photons const& v0photons, aod::V0Legs const& v0legs) + { + runPairing(collisions, v0photons, v0photons, v0legs, v0legs, + perCollisionPCM, perCollisionPCM, fV0PhotonCut, fV0PhotonCut); + ndf++; + } + PROCESS_SWITCH(photonhbt, processAnalysis, "pairing for analysis", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"photonhbt"})}; +}