diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h index f871586d5fa..9964f1d0333 100644 --- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h +++ b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h @@ -61,6 +61,7 @@ #include #include +#include #include #include #include @@ -104,6 +105,36 @@ struct Pi0EtaToGammaGamma { o2::framework::Configurable cfgAlphaMesonA{"cfgAlphaMesonA", 0.65, "photon energy asymmetry distribution parameter A for pT dependent cut (A * tanh(B*pT))"}; o2::framework::Configurable cfgAlphaMesonB{"cfgAlphaMesonB", 1.2, "photon energy asymmetry distribution parameter B for pT dependent cut (A * tanh(B*pT))"}; + o2::framework::Configurable cfgDoPhotonClassPairCut{"cfgDoPhotonClassPairCut", false, "apply photon-class pair selection A x B (leg track composition, PCM-PCM only)"}; + + struct : o2::framework::ConfigurableGroup { + std::string prefix = "photonclassA_group"; + o2::framework::Configurable cfgMinNLegsITSTPC{"cfgMinNLegsITSTPC", 0, "min number of ITS-TPC legs (0-2)"}; + o2::framework::Configurable cfgMaxNLegsITSTPC{"cfgMaxNLegsITSTPC", 2, "max number of ITS-TPC legs (0-2)"}; + o2::framework::Configurable cfgMinNLegsITSOnly{"cfgMinNLegsITSOnly", 0, "min number of ITS-only legs (0-2)"}; + o2::framework::Configurable cfgMaxNLegsITSOnly{"cfgMaxNLegsITSOnly", 2, "max number of ITS-only legs (0-2)"}; + o2::framework::Configurable cfgMinNLegsTPCOnly{"cfgMinNLegsTPCOnly", 0, "min number of TPC-only legs (0-2)"}; + o2::framework::Configurable cfgMaxNLegsTPCOnly{"cfgMaxNLegsTPCOnly", 2, "max number of TPC-only legs (0-2)"}; + o2::framework::Configurable cfgMinNLegsTRD{"cfgMinNLegsTRD", 0, "min number of legs with TRD (0-2)"}; + o2::framework::Configurable cfgMaxNLegsTRD{"cfgMaxNLegsTRD", 2, "max number of legs with TRD (0-2)"}; + o2::framework::Configurable cfgMinNLegsTOF{"cfgMinNLegsTOF", 0, "min number of legs with TOF (0-2)"}; + o2::framework::Configurable cfgMaxNLegsTOF{"cfgMaxNLegsTOF", 2, "max number of legs with TOF (0-2)"}; + } photonclassA; + + struct : o2::framework::ConfigurableGroup { + std::string prefix = "photonclassB_group"; + o2::framework::Configurable cfgMinNLegsITSTPC{"cfgMinNLegsITSTPC", 0, "min number of ITS-TPC legs (0-2)"}; + o2::framework::Configurable cfgMaxNLegsITSTPC{"cfgMaxNLegsITSTPC", 2, "max number of ITS-TPC legs (0-2)"}; + o2::framework::Configurable cfgMinNLegsITSOnly{"cfgMinNLegsITSOnly", 0, "min number of ITS-only legs (0-2)"}; + o2::framework::Configurable cfgMaxNLegsITSOnly{"cfgMaxNLegsITSOnly", 2, "max number of ITS-only legs (0-2)"}; + o2::framework::Configurable cfgMinNLegsTPCOnly{"cfgMinNLegsTPCOnly", 0, "min number of TPC-only legs (0-2)"}; + o2::framework::Configurable cfgMaxNLegsTPCOnly{"cfgMaxNLegsTPCOnly", 2, "max number of TPC-only legs (0-2)"}; + o2::framework::Configurable cfgMinNLegsTRD{"cfgMinNLegsTRD", 0, "min number of legs with TRD (0-2)"}; + o2::framework::Configurable cfgMaxNLegsTRD{"cfgMaxNLegsTRD", 2, "max number of legs with TRD (0-2)"}; + o2::framework::Configurable cfgMinNLegsTOF{"cfgMinNLegsTOF", 0, "min number of legs with TOF (0-2)"}; + o2::framework::Configurable cfgMaxNLegsTOF{"cfgMaxNLegsTOF", 2, "max number of legs with TOF (0-2)"}; + } photonclassB; + EMPhotonEventCut fEMEventCut; struct : o2::framework::ConfigurableGroup { std::string prefix = "eventcut_group"; @@ -371,6 +402,9 @@ struct Pi0EtaToGammaGamma { } }; + o2::aod::pwgem::photonmeson::utils::pairutil::V0PhotonClassSelection mPhotonClassSelA{}; + o2::aod::pwgem::photonmeson::utils::pairutil::V0PhotonClassSelection mPhotonClassSelB{}; + void init(o2::framework::InitContext&) { zvtx_bin_edges = std::vector(ConfVtxBins.value.begin(), ConfVtxBins.value.end()); @@ -397,6 +431,9 @@ struct Pi0EtaToGammaGamma { } DefineEMEventCut(); DefinePCMCut(); + + mPhotonClassSelA = o2::aod::pwgem::photonmeson::utils::pairutil::buildV0PhotonClassSelection(photonclassA); + mPhotonClassSelB = o2::aod::pwgem::photonmeson::utils::pairutil::buildV0PhotonClassSelection(photonclassB); DefineDileptonCut(); DefineEMCCut(); DefinePHOSCut(); @@ -799,6 +836,7 @@ struct Pi0EtaToGammaGamma { continue; } } + auto pos1 = g1.template posTrack_as(); auto ele1 = g1.template negTrack_as(); ROOT::Math::PtEtaPhiMVector v_gamma(g1.pt(), g1.eta(), g1.phi(), 0.); @@ -882,6 +920,20 @@ struct Pi0EtaToGammaGamma { } } + o2::aod::pwgem::photonmeson::utils::pairutil::V0PhotonLegCounts legCounts1{}, legCounts2{}; + if constexpr (std::is_same_v && std::is_same_v) { + auto pos1 = g1.template posTrack_as(); + auto neg1 = g1.template negTrack_as(); + auto pos2 = g2.template posTrack_as(); + auto neg2 = g2.template negTrack_as(); + legCounts1 = o2::aod::pwgem::photonmeson::utils::pairutil::getV0PhotonLegCounts(pos1, neg1); + legCounts2 = o2::aod::pwgem::photonmeson::utils::pairutil::getV0PhotonLegCounts(pos2, neg2); + if (cfgDoPhotonClassPairCut.value && + !o2::aod::pwgem::photonmeson::utils::pairutil::isPairPhotonClassSelected(legCounts1, legCounts2, mPhotonClassSelA, mPhotonClassSelB)) { + continue; + } + } + ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; @@ -920,11 +972,19 @@ struct Pi0EtaToGammaGamma { fRegistry.fill(HIST("Pair/same/hs"), v12.M(), v12.Pt(), wpair); if (std::find(used_photonIds_per_col.begin(), used_photonIds_per_col.end(), g1.globalIndex()) == used_photonIds_per_col.end()) { - emh1->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::photonmeson::utils::EMPhoton(g1.pt(), g1.eta(), g1.phi(), 0)); + auto emphoton1 = o2::aod::pwgem::photonmeson::utils::EMPhoton(g1.pt(), g1.eta(), g1.phi(), 0); + if constexpr (std::is_same_v && std::is_same_v) { + emphoton1.setLegCounts(legCounts1); + } + emh1->AddTrackToEventPool(key_df_collision, emphoton1); 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()) { - emh2->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::photonmeson::utils::EMPhoton(g2.pt(), g2.eta(), g2.phi(), 0)); + auto emphoton2 = o2::aod::pwgem::photonmeson::utils::EMPhoton(g2.pt(), g2.eta(), g2.phi(), 0); + if constexpr (std::is_same_v && std::is_same_v) { + emphoton2.setLegCounts(legCounts1); + } + emh2->AddTrackToEventPool(key_df_collision, emphoton2); used_photonIds_per_col.emplace_back(g2.globalIndex()); } ndiphoton++; @@ -969,6 +1029,14 @@ struct Pi0EtaToGammaGamma { for (const auto& g1 : selected_photons1_in_this_event) { for (const auto& g2 : photons1_from_event_pool) { + if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPCM) { + if (cfgDoPhotonClassPairCut.value) { + if (!g1.hasLegCounts() || !g2.hasLegCounts() || + !o2::aod::pwgem::photonmeson::utils::pairutil::isPairPhotonClassSelected(g1.legCounts(), g2.legCounts(), mPhotonClassSelA, mPhotonClassSelB)) { + continue; + } + } + } ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h index 95873165f71..f8ee7db2af0 100644 --- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h +++ b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h @@ -59,6 +59,7 @@ #include #include +#include #include #include #include @@ -95,6 +96,36 @@ struct Pi0EtaToGammaGammaMC { o2::framework::Configurable cfgAlphaMesonB{"cfgAlphaMesonB", 1.2, "photon energy asymmetry distribution parameter B for pT dependent cut (A * tanh(B*pT))"}; o2::framework::Configurable cfgGGContaCheck{"cfgGGContaCheck", false, "check gamma gamma contamination of dalitz"}; + o2::framework::Configurable cfgDoPhotonClassPairCut{"cfgDoPhotonClassPairCut", false, "apply photon-class pair selection A x B (leg track composition, PCM-PCM only)"}; + + struct : o2::framework::ConfigurableGroup { + std::string prefix = "photonclassA_group"; + o2::framework::Configurable cfgMinNLegsITSTPC{"cfgMinNLegsITSTPC", 0, "min number of ITS-TPC legs (0-2)"}; + o2::framework::Configurable cfgMaxNLegsITSTPC{"cfgMaxNLegsITSTPC", 2, "max number of ITS-TPC legs (0-2)"}; + o2::framework::Configurable cfgMinNLegsITSOnly{"cfgMinNLegsITSOnly", 0, "min number of ITS-only legs (0-2)"}; + o2::framework::Configurable cfgMaxNLegsITSOnly{"cfgMaxNLegsITSOnly", 2, "max number of ITS-only legs (0-2)"}; + o2::framework::Configurable cfgMinNLegsTPCOnly{"cfgMinNLegsTPCOnly", 0, "min number of TPC-only legs (0-2)"}; + o2::framework::Configurable cfgMaxNLegsTPCOnly{"cfgMaxNLegsTPCOnly", 2, "max number of TPC-only legs (0-2)"}; + o2::framework::Configurable cfgMinNLegsTRD{"cfgMinNLegsTRD", 0, "min number of legs with TRD (0-2)"}; + o2::framework::Configurable cfgMaxNLegsTRD{"cfgMaxNLegsTRD", 2, "max number of legs with TRD (0-2)"}; + o2::framework::Configurable cfgMinNLegsTOF{"cfgMinNLegsTOF", 0, "min number of legs with TOF (0-2)"}; + o2::framework::Configurable cfgMaxNLegsTOF{"cfgMaxNLegsTOF", 2, "max number of legs with TOF (0-2)"}; + } photonclassA; + + struct : o2::framework::ConfigurableGroup { + std::string prefix = "photonclassB_group"; + o2::framework::Configurable cfgMinNLegsITSTPC{"cfgMinNLegsITSTPC", 0, "min number of ITS-TPC legs (0-2)"}; + o2::framework::Configurable cfgMaxNLegsITSTPC{"cfgMaxNLegsITSTPC", 2, "max number of ITS-TPC legs (0-2)"}; + o2::framework::Configurable cfgMinNLegsITSOnly{"cfgMinNLegsITSOnly", 0, "min number of ITS-only legs (0-2)"}; + o2::framework::Configurable cfgMaxNLegsITSOnly{"cfgMaxNLegsITSOnly", 2, "max number of ITS-only legs (0-2)"}; + o2::framework::Configurable cfgMinNLegsTPCOnly{"cfgMinNLegsTPCOnly", 0, "min number of TPC-only legs (0-2)"}; + o2::framework::Configurable cfgMaxNLegsTPCOnly{"cfgMaxNLegsTPCOnly", 2, "max number of TPC-only legs (0-2)"}; + o2::framework::Configurable cfgMinNLegsTRD{"cfgMinNLegsTRD", 0, "min number of legs with TRD (0-2)"}; + o2::framework::Configurable cfgMaxNLegsTRD{"cfgMaxNLegsTRD", 2, "max number of legs with TRD (0-2)"}; + o2::framework::Configurable cfgMinNLegsTOF{"cfgMinNLegsTOF", 0, "min number of legs with TOF (0-2)"}; + o2::framework::Configurable cfgMaxNLegsTOF{"cfgMaxNLegsTOF", 2, "max number of legs with TOF (0-2)"}; + } photonclassB; + EMPhotonEventCut fEMEventCut; struct : o2::framework::ConfigurableGroup { std::string prefix = "eventcut_group"; @@ -116,6 +147,9 @@ struct Pi0EtaToGammaGammaMC { o2::framework::Configurable onlyKeepWeightedEvents{"onlyKeepWeightedEvents", false, "flag to keep only weighted events (for JJ MCs) and remove all MB events (with weight = 1)"}; } eventcuts; + o2::aod::pwgem::photonmeson::utils::pairutil::V0PhotonClassSelection mPhotonClassSelA{}; + o2::aod::pwgem::photonmeson::utils::pairutil::V0PhotonClassSelection mPhotonClassSelB{}; + V0PhotonCut fV0PhotonCut; struct : o2::framework::ConfigurableGroup { std::string prefix = "pcmcut_group"; @@ -249,6 +283,8 @@ struct Pi0EtaToGammaGammaMC { } DefineEMEventCut(); DefinePCMCut(); + mPhotonClassSelA = o2::aod::pwgem::photonmeson::utils::pairutil::buildV0PhotonClassSelection(photonclassA); + mPhotonClassSelB = o2::aod::pwgem::photonmeson::utils::pairutil::buildV0PhotonClassSelection(photonclassB); DefineDileptonCut(); DefineEMCCut(); DefinePHOSCut(); @@ -650,6 +686,14 @@ struct Pi0EtaToGammaGammaMC { auto pos2 = g2.template posTrack_as(); auto ele2 = g2.template negTrack_as(); + if (cfgDoPhotonClassPairCut.value && + !o2::aod::pwgem::photonmeson::utils::pairutil::isPairPhotonClassSelected( + o2::aod::pwgem::photonmeson::utils::pairutil::getV0PhotonLegCounts(pos1, ele1), + o2::aod::pwgem::photonmeson::utils::pairutil::getV0PhotonLegCounts(pos2, ele2), + mPhotonClassSelA, mPhotonClassSelB)) { + continue; + } + auto pos1mc = pos1.template emmcparticle_as(); auto ele1mc = ele1.template emmcparticle_as(); auto pos2mc = pos2.template emmcparticle_as(); diff --git a/PWGEM/PhotonMeson/Tasks/photonhbt.cxx b/PWGEM/PhotonMeson/Tasks/photonhbt.cxx index ae1d1017870..0ae897ca70f 100644 --- a/PWGEM/PhotonMeson/Tasks/photonhbt.cxx +++ b/PWGEM/PhotonMeson/Tasks/photonhbt.cxx @@ -20,6 +20,7 @@ #include "PWGEM/PhotonMeson/DataModel/EventTables.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/EventHistograms.h" +#include "PWGEM/PhotonMeson/Utils/PairUtilities.h" #include "Common/Core/RecoDecay.h" #include "Common/DataModel/Centrality.h" @@ -74,6 +75,7 @@ using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::soa; using namespace o2::aod::pwgem::dilepton::utils; +namespace pairutil = o2::aod::pwgem::photonmeson::utils::pairutil; // ─── Event Information Tables ──────────────────────────────────────────── @@ -102,11 +104,13 @@ static constexpr float kCmPerM = 100.f; [[nodiscard]] static float legPhiStar(float phi, float pt, int charge, float bzT, float radiusM) { - if (pt < kMinSigma) + if (pt < kMinSigma) { return kPhiStarInvalid; + } const float arg = -0.3f * bzT * static_cast(charge) * radiusM / (2.f * pt); - if (std::fabs(arg) >= 1.f) - return kPhiStarInvalid; // leg curls up before reaching this radius -> cannot merge there + if (std::fabs(arg) >= 1.f) { + return kPhiStarInvalid; // leg curls up before + } return RecoDecay::constrainAngle(phi + std::asin(arg), 0.f); } @@ -169,6 +173,8 @@ struct Photonhbt { [[nodiscard]] float legEta(int i) const { return fLegEta[i]; } [[nodiscard]] float legPhi(int i) const { return fLegPhi[i]; } [[nodiscard]] float legPt(int i) const { return fLegPt[i]; } + pairutil::V0PhotonLegCounts fLegCounts{}; + [[nodiscard]] pairutil::V0PhotonLegCounts const& legCounts() const { return fLegCounts; } }; struct PhotonMCInfo { @@ -196,7 +202,7 @@ struct Photonhbt { // CONFIGURABLES /*************************************************/ - Service ccdb; + Service ccdb{}; Configurable cfgCcdbUrl{"cfgCcdbUrl", "http://alice-ccdb.cern.ch", "CCDB url"}; Configurable cfgBzOverrideT{"cfgBzOverrideT", -999.f, "Bz in Tesla; used instead of CCDB if > -100"}; float mBzT{0.f}; // signed, Tesla @@ -220,6 +226,39 @@ struct Photonhbt { Configurable cfgUseLCMS{"cfgUseLCMS", false, "measure 1D relative momentum in LCMS"}; } hbtanalysis; + // ----- Photon Leg Classification + + Configurable cfgDoPhotonClassPairCut{"cfgDoPhotonClassPairCut", false, + "apply photon-class pair selection A x B (leg track composition)"}; + + struct : ConfigurableGroup { + std::string prefix = "photonclassA_group"; + Configurable cfgMinNLegsITSTPC{"cfgMinNLegsITSTPC", 0, "min number of ITS-TPC legs (0-2)"}; + Configurable cfgMaxNLegsITSTPC{"cfgMaxNLegsITSTPC", 2, "max number of ITS-TPC legs (0-2)"}; + Configurable cfgMinNLegsITSOnly{"cfgMinNLegsITSOnly", 0, "min number of ITS-only legs (0-2)"}; + Configurable cfgMaxNLegsITSOnly{"cfgMaxNLegsITSOnly", 2, "max number of ITS-only legs (0-2)"}; + Configurable cfgMinNLegsTPCOnly{"cfgMinNLegsTPCOnly", 0, "min number of TPC-only legs (0-2)"}; + Configurable cfgMaxNLegsTPCOnly{"cfgMaxNLegsTPCOnly", 2, "max number of TPC-only legs (0-2)"}; + Configurable cfgMinNLegsTRD{"cfgMinNLegsTRD", 0, "min number of legs with TRD (0-2)"}; + Configurable cfgMaxNLegsTRD{"cfgMaxNLegsTRD", 2, "max number of legs with TRD (0-2)"}; + Configurable cfgMinNLegsTOF{"cfgMinNLegsTOF", 0, "min number of legs with TOF (0-2)"}; + Configurable cfgMaxNLegsTOF{"cfgMaxNLegsTOF", 2, "max number of legs with TOF (0-2)"}; + } photonclassA; + + struct : ConfigurableGroup { + std::string prefix = "photonclassB_group"; + Configurable cfgMinNLegsITSTPC{"cfgMinNLegsITSTPC", 0, "min number of ITS-TPC legs (0-2)"}; + Configurable cfgMaxNLegsITSTPC{"cfgMaxNLegsITSTPC", 2, "max number of ITS-TPC legs (0-2)"}; + Configurable cfgMinNLegsITSOnly{"cfgMinNLegsITSOnly", 0, "min number of ITS-only legs (0-2)"}; + Configurable cfgMaxNLegsITSOnly{"cfgMaxNLegsITSOnly", 2, "max number of ITS-only legs (0-2)"}; + Configurable cfgMinNLegsTPCOnly{"cfgMinNLegsTPCOnly", 0, "min number of TPC-only legs (0-2)"}; + Configurable cfgMaxNLegsTPCOnly{"cfgMaxNLegsTPCOnly", 2, "max number of TPC-only legs (0-2)"}; + Configurable cfgMinNLegsTRD{"cfgMinNLegsTRD", 0, "min number of legs with TRD (0-2)"}; + Configurable cfgMaxNLegsTRD{"cfgMaxNLegsTRD", 2, "max number of legs with TRD (0-2)"}; + Configurable cfgMinNLegsTOF{"cfgMinNLegsTOF", 0, "min number of legs with TOF (0-2)"}; + Configurable cfgMaxNLegsTOF{"cfgMaxNLegsTOF", 2, "max number of legs with TOF (0-2)"}; + } photonclassB; + // ─── Event mixing ───────────────────────────────────────────────────────────── struct : ConfigurableGroup { std::string prefix = "mixing_group"; @@ -372,30 +411,30 @@ struct Photonhbt { // 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)"}; - const AxisSpec axisR{confRBins, "R_{conv} (cm)"}; - const AxisSpec axisDeltaR{confDeltaRBins, "|R_{1}-R_{2}| (cm)"}; - const AxisSpec axisDeltaR3D{confDeltaR3DBins, "|#vec{r}_{1}-#vec{r}_{2}| (cm)"}; - const AxisSpec axisDeltaRxy{confDeltaRxyBins, "#Delta r_{xy} (cm)"}; - const AxisSpec axisZConv{confZConvBins, "z_{conv} (cm)"}; - const AxisSpec axisDeltaZ{confDeltaZBins, "#Delta z (cm)"}; - const AxisSpec axisOccupancy{confOccupancyQA, "occupancy"}; - const AxisSpec axisCentQA{confCentQABins, "centrality (%)"}; - const AxisSpec axisLegPt{confLegPtBins, "p_{T,leg} (GeV/c)"}; - const AxisSpec axisLegDR{confLegDRBins, "#DeltaR_{legs}"}; + AxisSpec axisKt{confKtBins, "k_{T} (GeV/c)"}; + AxisSpec axisQinv{confQBins, "q_{inv} (GeV/c)"}; + AxisSpec axisQabsLcms{confQBins, "|#bf{q}|^{LCMS} (GeV/c)"}; + AxisSpec axisQout{confQBins, "q_{out} (GeV/c)"}; + AxisSpec axisQside{confQBins, "q_{side} (GeV/c)"}; + AxisSpec axisQlong{confQBins, "q_{long} (GeV/c)"}; + AxisSpec axisPt{confPtBins, "p_{T} (GeV/c)"}; + AxisSpec axisEta{confEtaBins, "#eta"}; + AxisSpec axisPhi{confPhiBins, "#phi (rad)"}; + AxisSpec axisDeltaEta{confDeltaEtaBins, "#Delta#eta"}; + AxisSpec axisDeltaPhi{confDeltaPhiBins, "#Delta#phi (rad)"}; + AxisSpec axisEllipseVal{confEllipseValBins, "(#Delta#eta/#sigma_{#eta})^{2}+(#Delta#phi/#sigma_{#phi})^{2}"}; + AxisSpec axisCosTheta{confCosThetaBins, "cos(#theta*)"}; + AxisSpec axisOpeningAngle{confOpeningAngleBins, "Opening angle (rad)"}; + AxisSpec axisR{confRBins, "R_{conv} (cm)"}; + AxisSpec axisDeltaR{confDeltaRBins, "|R_{1}-R_{2}| (cm)"}; + AxisSpec axisDeltaR3D{confDeltaR3DBins, "|#vec{r}_{1}-#vec{r}_{2}| (cm)"}; + AxisSpec axisDeltaRxy{confDeltaRxyBins, "#Delta r_{xy} (cm)"}; + AxisSpec axisZConv{confZConvBins, "z_{conv} (cm)"}; + AxisSpec axisDeltaZ{confDeltaZBins, "#Delta z (cm)"}; + AxisSpec axisOccupancy{confOccupancyQA, "occupancy"}; + AxisSpec axisCentQA{confCentQABins, "centrality (%)"}; + AxisSpec axisLegPt{confLegPtBins, "p_{T,leg} (GeV/c)"}; + AxisSpec axisLegDR{confLegDRBins, "#DeltaR_{legs}"}; /*************************************************/ // MAIN SPECS @@ -407,8 +446,6 @@ struct Photonhbt { HistogramRegistry fRegistryPairMC{"pairMC", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; HistogramRegistry fRegistryTruthMC{"truthMC", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; - static constexpr std::array event_pair_types = {"same/", "mix/"}; - EMPhotonEventCut fEMEventCut; V0PhotonCut fV0PhotonCut; @@ -417,6 +454,8 @@ struct Photonhbt { int mRunNumber{0}; bool isMC = false; int ndf = 0; + pairutil::V0PhotonClassSelection mPhotonClassSelA{}; + pairutil::V0PhotonClassSelection mPhotonClassSelB{}; std::vector ztxBinEdges; std::vector centBinEdges; @@ -465,6 +504,8 @@ struct Photonhbt { emh2 = std::make_shared(mixing.ndepth); DefineEMEventCut(); DefinePCMCut(); + mPhotonClassSelA = pairutil::buildV0PhotonClassSelection(photonclassA); + mPhotonClassSelB = pairutil::buildV0PhotonClassSelection(photonclassB); addhistograms(); std::random_device seedGen; engine = std::mt19937(seedGen()); @@ -477,8 +518,9 @@ struct Photonhbt { template void initCCDB(TCollision const& collision) { - if (mRunNumber == collision.runNumber()) + if (mRunNumber == collision.runNumber()) { return; + } mRunNumber = collision.runNumber(); if (cfgBzOverrideT.value > -100.f) { // o2-linter: disable=magic-number (number in case B-field is overridden in case not fetched from CCDB) mBzT = cfgBzOverrideT.value; @@ -543,24 +585,37 @@ struct Photonhbt { inline bool isInsideEllipse(float deta, float dphi) const { - if (!ggpaircuts.cfgDoEllipseCut.value) + if (!ggpaircuts.cfgDoEllipseCut.value) { return false; + } const float sE = ggpaircuts.cfgEllipseSigEta.value; const float sP = ggpaircuts.cfgEllipseSigPhi.value; - if (sE < kMinSigma || sP < kMinSigma) + if (sE < kMinSigma || sP < kMinSigma) { return false; + } return (deta / sE) * (deta / sE) + (dphi / sP) * (dphi / sP) < ggpaircuts.cfgEllipseR2.value; } inline bool passRZCut(float deltaR, float deltaZ) const { - if (ggpaircuts.cfgDoRCut.value && deltaR < ggpaircuts.cfgMinDeltaR.value) + if (ggpaircuts.cfgDoRCut.value && deltaR < ggpaircuts.cfgMinDeltaR.value) { return false; - if (ggpaircuts.cfgDoZCut.value && std::fabs(deltaZ) < ggpaircuts.cfgMinDeltaZ.value) + } + if (ggpaircuts.cfgDoZCut.value && std::fabs(deltaZ) < ggpaircuts.cfgMinDeltaZ.value) { return false; + } return true; } + [[nodiscard]] inline bool passPhotonClassPairCut(pairutil::V0PhotonLegCounts const& c1, + pairutil::V0PhotonLegCounts const& c2) const + { + if (!cfgDoPhotonClassPairCut.value) { + return true; + } + return pairutil::isPairPhotonClassSelected(c1, c2, mPhotonClassSelA, mPhotonClassSelB); + } + inline bool passAsymmetryCut(float pt1, float pt2) const { if (ggpaircuts.cfgMaxAsymmetry.value < 0.f) { @@ -998,8 +1053,9 @@ struct Photonhbt { 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() < kMinSigma || p1cmDir.R() < kMinSigma) + if (pairDir.R() < kMinSigma || p1cmDir.R() < kMinSigma) { return -1.f; + } return static_cast(pairDir.Unit().Dot(p1cmDir.Unit())); } @@ -1013,8 +1069,9 @@ struct Photonhbt { const auto xmin = static_cast(cfg.value[1]); const auto xmax = static_cast(cfg.value[2]); edges.resize(n + 1); - for (int i = 0; i <= n; ++i) + for (int i = 0; i <= n; ++i) { edges[i] = xmin + (xmax - xmin) / n * i; + } } } @@ -1051,8 +1108,9 @@ struct Photonhbt { template static constexpr const char* fullRangePrefix() { - if constexpr (ev_id == 0) + if constexpr (ev_id == 0) { return "Pair/same/FullRange/"; + } return "Pair/mix/FullRange/"; } @@ -1060,31 +1118,35 @@ struct Photonhbt { static constexpr const char* mcDirPrefix() { if constexpr (!IsMix) { - if constexpr (TruthT == PairTruthType::TrueTrueDistinct) + if constexpr (TruthT == PairTruthType::TrueTrueDistinct) { return "Pair/same/MC/TrueTrueDistinct/"; - else if constexpr (TruthT == PairTruthType::TrueTrueSamePhoton) + } else if constexpr (TruthT == PairTruthType::TrueTrueSamePhoton) { return "Pair/same/MC/TrueTrueSamePhoton/"; - else if constexpr (TruthT == PairTruthType::SharedMcLeg) + } else if constexpr (TruthT == PairTruthType::SharedMcLeg) { return "Pair/same/MC/SharedMcLeg/"; - else if constexpr (TruthT == PairTruthType::TrueFake) + } else if constexpr (TruthT == PairTruthType::TrueFake) { return "Pair/same/MC/TrueFake/"; - else if constexpr (TruthT == PairTruthType::FakeFake) + } else if constexpr (TruthT == PairTruthType::FakeFake) { return "Pair/same/MC/FakeFake/"; - else + } else { return "Pair/same/MC/Pi0Daughters/"; + } } else { - if constexpr (TruthT == PairTruthType::TrueTrueDistinct) + if constexpr (TruthT == PairTruthType::TrueTrueDistinct) { return "Pair/mix/MC/TrueTrueDistinct/"; - else if constexpr (TruthT == PairTruthType::TrueTrueSamePhoton) + + } else if constexpr (TruthT == PairTruthType::TrueTrueSamePhoton) { return "Pair/mix/MC/TrueTrueSamePhoton/"; - else if constexpr (TruthT == PairTruthType::SharedMcLeg) + } else if constexpr (TruthT == PairTruthType::SharedMcLeg) { return "Pair/mix/MC/SharedMcLeg/"; - else if constexpr (TruthT == PairTruthType::TrueFake) + } else if constexpr (TruthT == PairTruthType::TrueFake) { return "Pair/mix/MC/TrueFake/"; - else if constexpr (TruthT == PairTruthType::FakeFake) + + } else if constexpr (TruthT == PairTruthType::FakeFake) { return "Pair/mix/MC/FakeFake/"; - else + } else { return "Pair/mix/MC/Pi0Daughters/"; + } } } @@ -1128,6 +1190,7 @@ struct Photonhbt { p.fLegPt = {static_cast(pos.pt()), static_cast(ele.pt())}; p.fLegEta = {static_cast(pos.eta()), static_cast(ele.eta())}; p.fLegPhi = {static_cast(pos.phi()), static_cast(ele.phi())}; + p.fLegCounts = pairutil::getV0PhotonLegCounts(pos, ele); const float rConvCm = std::hypot(p.fVx, p.fVy); for (size_t ir = 0; ir < kPhiStarRadiiM.size(); ++ir) { if (kPhiStarRadiiM[ir] * kCmPerM < rConvCm) { // leg does not exist below R_conv @@ -1155,8 +1218,9 @@ struct Photonhbt { int n = 0; for (size_t ir = 0; ir < kPhiStarRadiiM.size(); ++ir) { const float p1 = a.fLegPhiStar[i][ir], p2 = b.fLegPhiStar[i][ir]; - if (p1 < -100.f || p2 < -100.f) // o2-linter: disable=magic-number (skip if non-sensical value is chosen) + if (p1 < -100.f || p2 < -100.f) { // o2-linter: disable=magic-number (skip if non-sensical value is chosen) continue; + } sumDp += RecoDecay::constrainAngle(p1 - p2, -o2::constants::math::PI); ++n; } @@ -1175,8 +1239,9 @@ struct Photonhbt { [[nodiscard]] inline bool passLegSepCut(PairSep const& s) const { - if (!pairsep.cfgDoLegSepCut.value) + if (!pairsep.cfgDoLegSepCut.value) { return true; + } return s.u > pairsep.cfgUCut.value; } @@ -1194,17 +1259,20 @@ struct Photonhbt { int n = 0; for (const float& r : kPhiStarRadiiM) { if ((g1.rTrue >= 0.f && r * kCmPerM < g1.rTrue) || - (g2.rTrue >= 0.f && r * kCmPerM < g2.rTrue)) + (g2.rTrue >= 0.f && r * kCmPerM < g2.rTrue)) { continue; + } const float p1 = legPhiStar(g1.legPhiTrue[i], g1.legPtTrue[i], q, mBzT, r); const float p2 = legPhiStar(g2.legPhiTrue[i], g2.legPtTrue[i], q, mBzT, r); - if (p1 < -100.f || p2 < -100.f) // o2-linter: disable=magic-number (skip if non-sensical value is chosen) + if (p1 < -100.f || p2 < -100.f) { // o2-linter: disable=magic-number (skip if non-sensical value is chosen) continue; + } sumDp += RecoDecay::constrainAngle(p1 - p2, -o2::constants::math::PI); ++n; } - if (n == 0) + if (n == 0) { return 999.f; + } return std::hypot(dEta / sE, (sumDp / static_cast(n)) / sP); }; return std::min(uSameSign(0), uSameSign(1)); @@ -1229,11 +1297,13 @@ struct Photonhbt { template inline void fillPairSep(PairSep const& s, PairQAObservables const& obs) { - if (!pairsep.cfgDoPairSepQA.value) + if (!pairsep.cfgDoPairSepQA.value) { return; + } const float limit = pairsep.cfgPairSepMaxQinv.value; - if (limit > 0.f && obs.qinv >= limit) + if (limit > 0.f && obs.qinv >= limit) { return; + } constexpr auto dir = pairSepPrefix(); fRegistryCF.fill(HIST(dir) + HIST("hDEtaDPhiKt_PP"), s.dEtaPP, s.dPhiPP, obs.kt); fRegistryCF.fill(HIST(dir) + HIST("hDEtaDPhiKt_NN"), s.dEtaNN, s.dPhiNN, obs.kt); @@ -1257,8 +1327,9 @@ struct Photonhbt { template inline void fillSinglePhotonQAStep(TPhoton const& g) { - if (!qaflags.doSinglePhotonQa) + if (!qaflags.doSinglePhotonQa) { return; + } constexpr auto base = singlePhotonQAPrefix(); const float r = std::sqrt(g.vx() * g.vx() + g.vy() * g.vy()); fRegistryPairQA.fill(HIST(base) + HIST("hEtaVsPhiPt"), g.phi(), g.eta(), g.pt()); @@ -1288,14 +1359,12 @@ struct Photonhbt { float qside_lcms = q3_lcms.Dot(uv_side); float qlong_lcms = q3_lcms.Dot(uv_long); if (hbtanalysis.cfgDo3D) { - fRegistryCF.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("CF_3D"), - std::fabs(qout_lcms), std::fabs(qside_lcms), std::fabs(qlong_lcms), kt, weight); - if (hbtanalysis.cfgDo2D) - fRegistryCF.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("CF_2D"), - std::fabs(qout_lcms), std::fabs(qinv), kt, weight); + fRegistryPairMC.fill(HIST("Pair/same/MC/NoLabel/CF_3D"), std::fabs(qout_lcms), std::fabs(qside_lcms), std::fabs(qlong_lcms), kt); + if (hbtanalysis.cfgDo2D) { + fRegistryPairMC.fill(HIST("Pair/same/MC/NoLabel/CF_2D"), std::fabs(qout_lcms), std::fabs(qinv), kt); + } } else { - fRegistryCF.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("CF_1D"), - hbtanalysis.cfgUseLCMS ? qabs_lcms : qinv, kt, weight); + fRegistryPairMC.fill(HIST("Pair/same/MC/NoLabel/CF_1D"), hbtanalysis.cfgUseLCMS ? qabs_lcms : qinv, kt); } float deta_pair = v1.Eta() - v2.Eta(); float dphi_pair = v1.Phi() - v2.Phi(); @@ -1333,8 +1402,9 @@ struct Photonhbt { if (hbtanalysis.cfgDo3D) { fRegistryPairMC.fill(HIST(mcDir) + HIST("CF_3D"), std::fabs(qout_lcms), std::fabs(qside_lcms), std::fabs(qlong_lcms), kt, weight); - if (hbtanalysis.cfgDo2D) + if (hbtanalysis.cfgDo2D) { fRegistryPairMC.fill(HIST(mcDir) + HIST("CF_2D"), std::fabs(qout_lcms), std::fabs(qinv), kt, weight); + } } else { fRegistryPairMC.fill(HIST(mcDir) + HIST("CF_1D"), hbtanalysis.cfgUseLCMS ? qabs_lcms : qinv, kt, weight); } @@ -1343,8 +1413,9 @@ struct Photonhbt { template inline void fillPairQAStep(PairQAObservables const& o, float /*cent*/, float /*occupancy*/) { - if (!qaflags.doPairQa) + if (!qaflags.doPairQa) { return; + } constexpr auto base = qaPrefix(); ///// Delta Eta QA @@ -1374,15 +1445,17 @@ struct Photonhbt { fRegistryPairQA.fill(HIST(base) + HIST("hDeltaR3DKt"), o.deltaR3D, o.kt); const float sE = ggpaircuts.cfgEllipseSigEta.value, sP = ggpaircuts.cfgEllipseSigPhi.value; - if (sE > kMinSigma && sP > kMinSigma) + if (sE > kMinSigma && sP > kMinSigma) { fRegistryPairQA.fill(HIST(base) + HIST("hEllipseVal"), (o.deta / sE) * (o.deta / sE) + (o.dphi / sP) * (o.dphi / sP)); + } } template inline void fillLegPairQAStep(LegPairObservables const& lo, float kt) { - if (!qaflags.doPairQa) + if (!qaflags.doPairQa) { return; + } constexpr auto base = qaPrefix<0, step_id>(); fRegistryPairQA.fill(HIST(base) + HIST("LegLS/hDEtaDPhi_PP"), lo.dEtaPP, lo.dPhiPP); @@ -1456,8 +1529,9 @@ struct Photonhbt { fRegistryPairMC.fill(HIST(base) + HIST("hDeltaZ"), obs.deltaZ); fRegistryPairMC.fill(HIST(base) + HIST("hDeltaR3D"), obs.deltaR3D); } - if (doSparse) + if (doSparse) { fRegistryPairMC.fill(HIST(base) + HIST("hSparse_DEtaDPhi_kT"), obs.deta, obs.dphi, obs.kt); + } constexpr auto summaryDir = IsMix ? "Pair/mix/MC/" : "Pair/same/MC/"; if (doMCQA) { fRegistryPairMC.fill(HIST(summaryDir) + HIST("hTruthTypeVsQinv"), obs.qinv, static_cast(TruthT)); @@ -1499,20 +1573,29 @@ struct Photonhbt { fRegistryPairMC.fill(HIST(base) + HIST("hDeltaRVsQinv"), obs.qinv, obs.deltaR); fRegistryPairMC.fill(HIST(base) + HIST("hDeltaZVsQinv"), obs.qinv, obs.deltaZ); fRegistryPairMC.fill(HIST(base) + HIST("hDeltaR3DVsQinv"), obs.qinv, obs.deltaR3D); - const bool fillDRDZ = ((TruthT == PairTruthType::TrueTrueDistinct) ? mctruthSparse.cfgFillDRDZQinvTrueTrueDistinct.value : (TruthT == PairTruthType::TrueTrueSamePhoton) ? mctruthSparse.cfgFillDRDZQinvTrueTrueSamePhoton.value - : (TruthT == PairTruthType::SharedMcLeg) ? mctruthSparse.cfgFillDRDZQinvSharedMcLeg.value - : (TruthT == PairTruthType::TrueFake) ? mctruthSparse.cfgFillDRDZQinvTrueFake.value - : (TruthT == PairTruthType::FakeFake) ? mctruthSparse.cfgFillDRDZQinvFakeFake.value - : mctruthSparse.cfgFillDRDZQinvPi0Daughters.value); - if (fillDRDZ) + + const bool fillDRDZ = + (TruthT == PairTruthType::TrueTrueDistinct) ? mctruthSparse.cfgFillDRDZQinvTrueTrueDistinct.value + : (TruthT == PairTruthType::TrueTrueSamePhoton) ? mctruthSparse.cfgFillDRDZQinvTrueTrueSamePhoton.value + : (TruthT == PairTruthType::SharedMcLeg) ? mctruthSparse.cfgFillDRDZQinvSharedMcLeg.value + : (TruthT == PairTruthType::TrueFake) ? mctruthSparse.cfgFillDRDZQinvTrueFake.value + : (TruthT == PairTruthType::FakeFake) ? mctruthSparse.cfgFillDRDZQinvFakeFake.value + : mctruthSparse.cfgFillDRDZQinvPi0Daughters.value; + + if (fillDRDZ) { fRegistryPairMC.fill(HIST(base) + HIST("hSparseDeltaRDeltaZQinv"), obs.deltaR, obs.deltaZ, obs.qinv); - const bool enabled = ((TruthT == PairTruthType::TrueTrueDistinct) ? mctruthSparse.cfgFillDEtaDPhiVsQinvTrueTrueDistinct.value : (TruthT == PairTruthType::TrueTrueSamePhoton) ? mctruthSparse.cfgFillDEtaDPhiVsQinvTrueTrueSamePhoton.value - : (TruthT == PairTruthType::SharedMcLeg) ? mctruthSparse.cfgFillDEtaDPhiVsQinvSharedMcLeg.value - : (TruthT == PairTruthType::TrueFake) ? mctruthSparse.cfgFillDEtaDPhiVsQinvTrueFake.value - : (TruthT == PairTruthType::FakeFake) ? mctruthSparse.cfgFillDEtaDPhiVsQinvFakeFake.value - : mctruthSparse.cfgFillDEtaDPhiVsQinvPi0Daughters.value); - if (enabled) + } + + const bool enabled = + (TruthT == PairTruthType::TrueTrueDistinct) ? mctruthSparse.cfgFillDEtaDPhiVsQinvTrueTrueDistinct.value + : (TruthT == PairTruthType::TrueTrueSamePhoton) ? mctruthSparse.cfgFillDEtaDPhiVsQinvTrueTrueSamePhoton.value + : (TruthT == PairTruthType::SharedMcLeg) ? mctruthSparse.cfgFillDEtaDPhiVsQinvSharedMcLeg.value + : (TruthT == PairTruthType::TrueFake) ? mctruthSparse.cfgFillDEtaDPhiVsQinvTrueFake.value + : (TruthT == PairTruthType::FakeFake) ? mctruthSparse.cfgFillDEtaDPhiVsQinvFakeFake.value + : mctruthSparse.cfgFillDEtaDPhiVsQinvPi0Daughters.value; + if (enabled) { fRegistryPairMC.fill(HIST(base) + HIST("hDEtaDPhiVsQinv"), obs.deta, obs.dphi, obs.qinv); + } } template @@ -1564,10 +1647,10 @@ struct Photonhbt { float qside_lcms = q3_lcms.Dot(uv_side); float qlong_lcms = q3_lcms.Dot(uv_long); if (hbtanalysis.cfgDo3D) { - fRegistryPairMC.fill(HIST("Pair/same/MC/NoLabel/CF_3D"), - std::fabs(qout_lcms), std::fabs(qside_lcms), std::fabs(qlong_lcms), kt); - if (hbtanalysis.cfgDo2D) - fRegistryPairMC.fill(HIST("Pair/same/MC/NoLabel/CF_2D"), std::fabs(qout_lcms), std::fabs(qinv), kt); + fRegistryPairMC.fill(HIST("Pair/same/MC/NoLabel/CF_3D"), std::fabs(qout_lcms), std::fabs(qside_lcms), std::fabs(qlong_lcms), kt); + } + if (hbtanalysis.cfgDo2D) { + fRegistryPairMC.fill(HIST("Pair/same/MC/NoLabel/CF_2D"), std::fabs(qout_lcms), std::fabs(qinv), kt); } else { fRegistryPairMC.fill(HIST("Pair/same/MC/NoLabel/CF_1D"), hbtanalysis.cfgUseLCMS ? qabs_lcms : qinv, kt); } @@ -1606,8 +1689,9 @@ struct Photonhbt { cosPA = std::clamp(cosPA, -1.f, 1.f); o.opa = std::acos(cosPA); o2::math_utils::bringTo02Pi(o.opa); - if (o.opa > o2::constants::math::PI) - o.opa -= o2::constants::math::PI; + if (o.opa > o2::constants::math::PI) { + o.opa = 2.f * o2::constants::math::PI - o.opa; + } o.cosOA = std::cos(o.opa / 2.f); o.drOverCosOA = (std::fabs(o.cosOA) < kMinCosine) ? 1e12f : (o.deltaR3D / o.cosOA); o.v1 = ROOT::Math::PtEtaPhiMVector(g1.pt(), g1.eta(), g1.phi(), 0.f); @@ -1630,18 +1714,21 @@ struct Photonhbt { PhotonMCInfo info{}; const auto pos = g.template posTrack_as(); const auto neg = g.template negTrack_as(); - if (!pos.has_emmcparticle() || !neg.has_emmcparticle()) + if (!pos.has_emmcparticle() || !neg.has_emmcparticle()) { return info; + } info.hasMC = true; info.mcPosId = pos.emmcparticleId(); info.mcNegId = neg.emmcparticleId(); const auto mcPos = pos.template emmcparticle_as(); const auto mcNeg = neg.template emmcparticle_as(); - if (!mcPos.has_mothers() || !mcNeg.has_mothers()) + if (!mcPos.has_mothers() || !mcNeg.has_mothers()) { return info; + } const int mothIdPos = mcPos.mothersIds()[0], mothIdNeg = mcNeg.mothersIds()[0]; - if (mothIdPos != mothIdNeg) + if (mothIdPos != mothIdNeg) { return info; + } info.sameMother = true; info.motherId = mothIdPos; const auto mother = mcParticles.iteratorAt(mothIdPos); @@ -1657,15 +1744,19 @@ struct Photonhbt { const bool t2 = m2.hasMC && m2.sameMother && m2.isTruePhoton; if (m1.hasMC && m2.hasMC) { if ((m1.mcPosId >= 0 && (m1.mcPosId == m2.mcPosId || m1.mcPosId == m2.mcNegId)) || - (m1.mcNegId >= 0 && (m1.mcNegId == m2.mcPosId || m1.mcNegId == m2.mcNegId))) + (m1.mcNegId >= 0 && (m1.mcNegId == m2.mcPosId || m1.mcNegId == m2.mcNegId))) { return PairTruthType::SharedMcLeg; + } } - if (!t1 && !t2) + if (!t1 && !t2) { return PairTruthType::FakeFake; - if (t1 != t2) + } + if (t1 != t2) { return PairTruthType::TrueFake; - if (m1.motherId >= 0 && m1.motherId == m2.motherId) + } + if (m1.motherId >= 0 && m1.motherId == m2.motherId) { return PairTruthType::TrueTrueSamePhoton; + } return PairTruthType::TrueTrueDistinct; } @@ -1673,15 +1764,18 @@ struct Photonhbt { static bool isPi0DaughterPair(PhotonMCInfo const& m1, PhotonMCInfo const& m2, TMCParticles const& mcParticles) { - if (!m1.isTruePhoton || !m2.isTruePhoton || m1.motherId < 0 || m2.motherId < 0) + if (!m1.isTruePhoton || !m2.isTruePhoton || m1.motherId < 0 || m2.motherId < 0) { return false; + } const auto ph1 = mcParticles.iteratorAt(m1.motherId); const auto ph2 = mcParticles.iteratorAt(m2.motherId); - if (!ph1.has_mothers() || !ph2.has_mothers()) + if (!ph1.has_mothers() || !ph2.has_mothers()) { return false; + } const int gm1 = ph1.mothersIds()[0], gm2 = ph2.mothersIds()[0]; - if (gm1 != gm2) + if (gm1 != gm2) { return false; + } return (std::abs(mcParticles.iteratorAt(gm1).pdgCode()) == kPi0); } @@ -1723,15 +1817,17 @@ struct Photonhbt { initCCDB(collision); int ndiphoton = 0; const std::array cent = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; - if (cent[mixing.cfgCentEstimator] < centralitySelection.cfgCentMin || centralitySelection.cfgCentMax < cent[mixing.cfgCentEstimator]) + if (cent[mixing.cfgCentEstimator] < centralitySelection.cfgCentMin || centralitySelection.cfgCentMax < cent[mixing.cfgCentEstimator]) { continue; + } const std::array epArr = {collision.ep2ft0m(), collision.ep2ft0a(), collision.ep2ft0c(), collision.ep2fv0a(), collision.ep2btot(), collision.ep2bpos(), collision.ep2bneg()}; const float ep2 = epArr[mixing.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)) + if (!fEMEventCut.IsSelected(collision)) { continue; + } o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(&fRegistry, collision, 1.f); fRegistry.fill(HIST("Event/before/hCollisionCounter"), 12.0); fRegistry.fill(HIST("Event/after/hCollisionCounter"), 12.0); @@ -1756,18 +1852,22 @@ struct Photonhbt { std::unordered_set idsAfterPairCuts; for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photons1Coll, photons2Coll))) { if (!cut1.template IsSelected(g1) || - !cut2.template IsSelected(g2)) + !cut2.template IsSelected(g2)) { continue; + } const auto pos1 = g1.template posTrack_as(), ele1 = g1.template negTrack_as(); const auto pos2 = g2.template posTrack_as(), ele2 = g2.template negTrack_as(); if (pos1.trackId() == pos2.trackId() || pos1.trackId() == ele2.trackId() || - ele1.trackId() == pos2.trackId() || ele1.trackId() == ele2.trackId()) + ele1.trackId() == pos2.trackId() || ele1.trackId() == ele2.trackId()) { continue; - if (!passAsymmetryCut(g1.pt(), g2.pt())) + } + if (!passAsymmetryCut(g1.pt(), g2.pt())) { continue; + } auto obs = buildPairQAObservables(g1, g2); - if (!obs.valid) + if (!obs.valid) { continue; + } const bool doQA = passQinvQAGate(obs.qinv), doFR = passQinvFullRangeGate(obs.qinv); const auto legObs = buildLegPairObservables(g1, g2, pos1, ele1, pos2, ele2); const auto pwl1 = makePhotonWithLegs(g1, pos1, ele1); @@ -1780,18 +1880,27 @@ struct Photonhbt { fillPairQAStep<0, 0>(obs, centForQA, occupancy); fillLegPairQAStep<0>(legObs, obs.kt); } - if (doFR) + if (doFR) { fillFullRangeDeltaRCosOA<0>(obs.qinv, obs.drOverCosOA); + } // ─── Pair cuts ──────────────────────────────────────────────────── - if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA) + if (!passPhotonClassPairCut(pwl1.legCounts(), pwl2.legCounts())) { continue; - if (!passLegSepCut(sep)) + } + if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA) { continue; - if (!passRZCut(obs.deltaR, obs.deltaZ)) + } + if (!passLegSepCut(sep)) { + continue; + } + if (!passRZCut(obs.deltaR, obs.deltaZ)) { continue; - if (isInsideEllipse(obs.deta, obs.dphi)) + } + + if (isInsideEllipse(obs.deta, obs.dphi)) { continue; + } // ───after pair cuts ────────────────────────────────────── fillPairSep<0, true>(sep, obs); @@ -1801,8 +1910,9 @@ struct Photonhbt { fillPairQAStep<0, 1>(obs, centForQA, occupancy); fillLegPairQAStep<1>(legObs, obs.kt); } - if (doFR) + if (doFR) { fillFullRangeQA<0>(obs, centForQA, occupancy); + } fillPairHistogram<0>(collision, obs.v1, obs.v2, 1.f); ndiphoton++; fRegistryCF.fill(HIST("Pair/same/hPhi_lowerPtV0"), (g1.pt() < g2.pt()) ? g1.phi() : g2.phi()); @@ -1825,51 +1935,67 @@ struct Photonhbt { } } usedPhotonIdsPerCol.clear(); - if (!mixing.cfgDoMix || ndiphoton == 0) + if (!mixing.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) + if (mixID.second == collision.globalIndex() && mixID.first == ndf) { continue; + } const uint64_t bcMix = mapMixedEventIdToGlobalBC[mixID]; const uint64_t diffBC = std::max(collision.globalBC(), bcMix) - std::min(collision.globalBC(), bcMix); fRegistryCF.fill(HIST("Pair/mix/hDiffBC"), diffBC); - if (diffBC < mixing.ndiffBCMix) + if (diffBC < mixing.ndiffBCMix) { continue; + } auto poolPhotons = emh1->GetTracksPerCollision(mixID); for (const auto& g1 : selectedPhotons) { for (const auto& g2 : poolPhotons) { - if (!passAsymmetryCut(g1.pt(), g2.pt())) + if (!passAsymmetryCut(g1.pt(), g2.pt())) { continue; + } auto obs = buildPairQAObservables(g1, g2); - if (!obs.valid) + if (!obs.valid) { continue; + } const bool doQA = passQinvQAGate(obs.qinv), doFR = passQinvFullRangeGate(obs.qinv); const auto sep = computePairSep(g1, g2); // ───before pair cuts ───────────────────────────────── fillPairSep<1>(sep, obs); - if (doQA) + if (doQA) { fillPairQAStep<1, 0>(obs, centForQA, occupancy); - if (doFR) + } + if (doFR) { fillFullRangeDeltaRCosOA<1>(obs.qinv, obs.drOverCosOA); + } // ─── Pair cuts ──────────────────────────────────────────────── - if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA) + if (!passPhotonClassPairCut(g1.legCounts(), g2.legCounts())) { continue; - if (!passLegSepCut(sep)) + } + if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA) { continue; - if (!passRZCut(obs.deltaR, obs.deltaZ)) + } + if (!passLegSepCut(sep)) { continue; - if (isInsideEllipse(obs.deta, obs.dphi)) + } + if (!passRZCut(obs.deltaR, obs.deltaZ)) { + continue; + } + if (isInsideEllipse(obs.deta, obs.dphi)) { continue; + } fillPairSep<1, true>(sep, obs); - if (doQA) + if (doQA) { fillPairQAStep<1, 1>(obs, centForQA, occupancy); - if (doFR) + } + if (doFR) { fillFullRangeQA<1>(obs, centForQA, occupancy); + } fillPairHistogram<1>(collision, obs.v1, obs.v2, 1.f); fRegistryCF.fill(HIST("Pair/mix/hPhi_lowerPtV0"), (g1.pt() < g2.pt()) ? g1.phi() : g2.phi()); } @@ -1898,15 +2024,17 @@ struct Photonhbt { initCCDB(collision); int ndiphoton = 0; const std::array cent = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; - if (cent[mixing.cfgCentEstimator] < centralitySelection.cfgCentMin || centralitySelection.cfgCentMax < cent[mixing.cfgCentEstimator]) + if (cent[mixing.cfgCentEstimator] < centralitySelection.cfgCentMin || centralitySelection.cfgCentMax < cent[mixing.cfgCentEstimator]) { continue; + } const std::array epArr = {collision.ep2ft0m(), collision.ep2ft0a(), collision.ep2ft0c(), collision.ep2fv0a(), collision.ep2btot(), collision.ep2bpos(), collision.ep2bneg()}; const float ep2 = epArr[mixing.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)) + if (!fEMEventCut.IsSelected(collision)) { continue; + } o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(&fRegistry, collision, 1.f); fRegistry.fill(HIST("Event/before/hCollisionCounter"), 12.0); fRegistry.fill(HIST("Event/after/hCollisionCounter"), 12.0); @@ -1929,23 +2057,28 @@ struct Photonhbt { } std::unordered_set idsAfterPairCuts; for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photonsColl, photonsColl))) { - if (!cut.template IsSelected(g1) || !cut.template IsSelected(g2)) + if (!cut.template IsSelected(g1) || !cut.template IsSelected(g2)) { continue; + } const auto pos1 = g1.template posTrack_as(), ele1 = g1.template negTrack_as(); const auto pos2 = g2.template posTrack_as(), ele2 = g2.template negTrack_as(); if (pos1.trackId() == pos2.trackId() || pos1.trackId() == ele2.trackId() || - ele1.trackId() == pos2.trackId() || ele1.trackId() == ele2.trackId()) + ele1.trackId() == pos2.trackId() || ele1.trackId() == ele2.trackId()) { continue; + } const auto mc1 = buildPhotonMCInfo(g1, mcParticles); const auto mc2 = buildPhotonMCInfo(g2, mcParticles); auto truthType = classifyPairTruth(mc1, mc2); - if (truthType == PairTruthType::TrueTrueDistinct && isPi0DaughterPair(mc1, mc2, mcParticles)) + if (truthType == PairTruthType::TrueTrueDistinct && isPi0DaughterPair(mc1, mc2, mcParticles)) { truthType = PairTruthType::Pi0Daughters; - if (!passAsymmetryCut(g1.pt(), g2.pt())) + } + if (!passAsymmetryCut(g1.pt(), g2.pt())) { continue; + } auto obs = buildPairQAObservables(g1, g2); - if (!obs.valid) + if (!obs.valid) { continue; + } const bool doQA = passQinvQAGate(obs.qinv), doFR = passQinvFullRangeGate(obs.qinv); const auto legObs = buildLegPairObservables(g1, g2, pos1, ele1, pos2, ele2); const auto pwl1 = makePhotonWithLegs(g1, pos1, ele1); @@ -1958,18 +2091,25 @@ struct Photonhbt { fillPairQAStep<0, 0>(obs, centForQA, occupancy); fillLegPairQAStep<0>(legObs, obs.kt); } - if (doFR) + if (doFR) { fillFullRangeDeltaRCosOA<0>(obs.qinv, obs.drOverCosOA); + } // ─── Pair cuts ──────────────────────────────────────────────────── - if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA) + if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA) { continue; - if (!passLegSepCut(sep)) + } + + if (!passLegSepCut(sep)) { continue; - if (!passRZCut(obs.deltaR, obs.deltaZ)) + } + + if (!passRZCut(obs.deltaR, obs.deltaZ)) { continue; - if (isInsideEllipse(obs.deta, obs.dphi)) + } + if (isInsideEllipse(obs.deta, obs.dphi)) { continue; + } // ─── after pair cuts ────────────────────────────────────── fillPairSep<0, true>(sep, obs); @@ -1979,8 +2119,9 @@ struct Photonhbt { fillPairQAStep<0, 1>(obs, centForQA, occupancy); fillLegPairQAStep<1>(legObs, obs.kt); } - if (doFR) + if (doFR) { fillFullRangeQA<0>(obs, centForQA, occupancy); + } fillPairHistogram<0>(collision, obs.v1, obs.v2, 1.f); fRegistryCF.fill(HIST("Pair/same/hPhi_lowerPtV0"), (g1.pt() < g2.pt()) ? g1.phi() : g2.phi()); ndiphoton++; @@ -1995,8 +2136,9 @@ struct Photonhbt { const bool doMCQA = passQinvMCQAGate(obs.qinv); fillMCPairQA(truthType, obs, doQA, doMCQA); fillLegPairMC(truthType, legObs, obs.kt); - if (doFR) + if (doFR) { fillMCPairQAFullRange(truthType, obs); + } const bool isTruePair = (truthType == PairTruthType::TrueTrueDistinct || truthType == PairTruthType::TrueTrueSamePhoton || truthType == PairTruthType::Pi0Daughters); @@ -2056,52 +2198,65 @@ struct Photonhbt { } } usedPhotonIdsPerCol.clear(); - if (!mixing.cfgDoMix || ndiphoton == 0) + if (!mixing.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) + if (mixID.second == collision.globalIndex() && mixID.first == ndf) { continue; + } const uint64_t bcMix = mapMixedEventIdToGlobalBC[mixID]; const uint64_t diffBC = std::max(collision.globalBC(), bcMix) - std::min(collision.globalBC(), bcMix); fRegistryCF.fill(HIST("Pair/mix/hDiffBC"), diffBC); - if (diffBC < mixing.ndiffBCMix) + if (diffBC < mixing.ndiffBCMix) { continue; + } auto poolPhotons = emh1->GetTracksPerCollision(mixID); for (const auto& g1 : selectedPhotons) { for (const auto& g2 : poolPhotons) { - if (!passAsymmetryCut(g1.pt(), g2.pt())) + if (!passAsymmetryCut(g1.pt(), g2.pt())) { continue; + } auto obs = buildPairQAObservables(g1, g2); - if (!obs.valid) + if (!obs.valid) { continue; + } const bool doQA = passQinvQAGate(obs.qinv), doFR = passQinvFullRangeGate(obs.qinv); const auto sep = computePairSep(g1, g2); // ───before pair cuts ───────────────────────────────── fillPairSep<1>(sep, obs); - if (doQA) + if (doQA) { fillPairQAStep<1, 0>(obs, centForQA, occupancy); - if (doFR) + } + if (doFR) { fillFullRangeDeltaRCosOA<1>(obs.qinv, obs.drOverCosOA); + } // ─── Pair cuts ──────────────────────────────────────────────── - if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA) + if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA) { continue; - if (!passLegSepCut(sep)) + } + if (!passLegSepCut(sep)) { continue; - if (!passRZCut(obs.deltaR, obs.deltaZ)) + } + if (!passRZCut(obs.deltaR, obs.deltaZ)) { continue; - if (isInsideEllipse(obs.deta, obs.dphi)) + } + if (isInsideEllipse(obs.deta, obs.dphi)) { continue; + } // ─after pair cuts ────────────────────────────────── fillPairSep<1, true>(sep, obs); - if (doQA) + if (doQA) { fillPairQAStep<1, 1>(obs, centForQA, occupancy); - if (doFR) + } + if (doFR) { fillFullRangeQA<1>(obs, centForQA, occupancy); + } fillPairHistogram<1>(collision, obs.v1, obs.v2, 1.f); fRegistryCF.fill(HIST("Pair/mix/hPhi_lowerPtV0"), (g1.pt() < g2.pt()) ? g1.phi() : g2.phi()); } @@ -2134,20 +2289,24 @@ struct Photonhbt { std::unordered_set mcIdsWithAnyV0Leg; mcIdsWithAnyV0Leg.reserve(v0legs.size() * 2); for (const auto& leg : v0legs) { - if (leg.has_emmcparticle()) + if (leg.has_emmcparticle()) { mcIdsWithAnyV0Leg.insert(leg.emmcparticleId()); + } } for (const auto& collision : collisions) { initCCDB(collision); - if (!fEMEventCut.IsSelected(collision)) + if (!fEMEventCut.IsSelected(collision)) { continue; + } const std::array cent = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; if (cent[mixing.cfgCentEstimator] < centralitySelection.cfgCentMin || - centralitySelection.cfgCentMax < cent[mixing.cfgCentEstimator]) + centralitySelection.cfgCentMax < cent[mixing.cfgCentEstimator]) { continue; - if (!collision.has_emmcevent()) + } + if (!collision.has_emmcevent()) { continue; + } const int64_t thisCollisionId = collision.globalIndex(); const int mcEventId = collision.template emmcevent_as().globalIndex(); @@ -2158,14 +2317,18 @@ struct Photonhbt { std::unordered_set legIdsThisCollision; legIdsThisCollision.reserve(legsColl.size()); - for (const auto& leg : legsColl) - if (leg.has_emmcparticle()) + for (const auto& leg : legsColl) { + if (leg.has_emmcparticle()) { legIdsThisCollision.insert(leg.emmcparticleId()); + } + } std::unordered_map> crossBuildMap; struct PhotonRecoInfo { bool hasV0 = false, passesCut = false; + bool hasCounts = false; + pairutil::V0PhotonLegCounts legCounts{}; }; std::unordered_map gammaRecoMap; gammaRecoMap.reserve(recoPhotonsColl.size()); @@ -2173,14 +2336,17 @@ struct Photonhbt { for (const auto& g : recoPhotonsColl) { const auto pos = g.template posTrack_as(); const auto neg = g.template negTrack_as(); - if (pos.collisionId() != thisCollisionId || neg.collisionId() != thisCollisionId) + if (pos.collisionId() != thisCollisionId || neg.collisionId() != thisCollisionId) { continue; - if (!pos.has_emmcparticle() || !neg.has_emmcparticle()) + } + if (!pos.has_emmcparticle() || !neg.has_emmcparticle()) { continue; + } const auto mcPos = pos.template emmcparticle_as(); const auto mcNeg = neg.template emmcparticle_as(); - if (!mcPos.has_mothers() || !mcNeg.has_mothers()) + if (!mcPos.has_mothers() || !mcNeg.has_mothers()) { continue; + } const int posMotherId = mcPos.mothersIds()[0], negMotherId = mcNeg.mothersIds()[0]; if (posMotherId != negMotherId) { const auto posMother = emmcParticles.iteratorAt(posMotherId); @@ -2192,12 +2358,17 @@ struct Photonhbt { continue; } const int gammaId = posMotherId; - if (emmcParticles.iteratorAt(gammaId).pdgCode() != kGamma) + if (emmcParticles.iteratorAt(gammaId).pdgCode() != kGamma) { continue; + } const bool passes = cut.template IsSelected, TLegs>(g); auto& info = gammaRecoMap[gammaId]; info.hasV0 = true; info.passesCut = info.passesCut || passes; + if (passes && !info.hasCounts) { + info.legCounts = pairutil::getV0PhotonLegCounts(pos, neg); + info.hasCounts = true; + } } // ─── Build true gamma list ──────────────────────────────────────────────── @@ -2205,19 +2376,24 @@ struct Photonhbt { trueGammas.reserve(32); for (const auto& g : emmcPartsColl) { - if (g.pdgCode() != kGamma) + if (g.pdgCode() != kGamma) { continue; - if (!g.isPhysicalPrimary() && !g.producedByGenerator()) + } + if (!g.isPhysicalPrimary() && !g.producedByGenerator()) { continue; - if (std::fabs(g.eta()) > pcmcuts.cfgMaxEtaV0.value) + } + if (std::fabs(g.eta()) > pcmcuts.cfgMaxEtaV0.value) { continue; + } const float mcV0PtMin = (mctruth.cfgMCMinV0Pt.value > 0.f) ? mctruth.cfgMCMinV0Pt.value : pcmcuts.cfgMinPtV0.value; - if (g.pt() < mcV0PtMin) + if (g.pt() < mcV0PtMin) { continue; - if (!g.has_daughters()) + } + if (!g.has_daughters()) { continue; + } int posId = -1, negId = -1; float rTrue = -1.f; @@ -2242,16 +2418,17 @@ struct Photonhbt { if (mctruth.cfgMCMinLegPt.value > 0.f && (static_cast(mcPosE.pt()) < mctruth.cfgMCMinLegPt.value || - static_cast(mcNegE.pt()) < mctruth.cfgMCMinLegPt.value)) + static_cast(mcNegE.pt()) < mctruth.cfgMCMinLegPt.value)) { continue; + } - const float deTrE = static_cast(mcPosE.eta() - mcNegE.eta()); - const float dpTrE = wrapPhi(static_cast(mcPosE.phi() - mcNegE.phi())); - const float legDRt = std::sqrt(deTrE * deTrE + dpTrE * dpTrE); + const auto deTrE = static_cast(mcPosE.eta() - mcNegE.eta()); + const auto dpTrE = wrapPhi(static_cast(mcPosE.phi() - mcNegE.phi())); + const auto legDRt = std::sqrt(deTrE * deTrE + dpTrE * dpTrE); - const float pxG = static_cast(g.px()), pyG = static_cast(g.py()), - pzG = static_cast(g.pz()); - const float magG = std::sqrt(pxG * pxG + pyG * pyG + pzG * pzG); + const auto pxG = static_cast(g.px()), pyG = static_cast(g.py()), + pzG = static_cast(g.pz()); + const auto magG = std::sqrt(pxG * pxG + pyG * pyG + pzG * pzG); float alphaTrue = 0.f; if (magG > kMinSigma) { const float ux = pxG / magG, uy = pyG / magG, uz = pzG / magG; @@ -2262,8 +2439,9 @@ struct Photonhbt { static_cast(mcNegE.py()) * uy + static_cast(mcNegE.pz()) * uz; const float sumPL = pLpos + pLneg; - if (std::fabs(sumPL) > kMinSigma) + if (std::fabs(sumPL) > kMinSigma) { alphaTrue = (pLpos - pLneg) / sumPL; + } } TruthGamma tg; @@ -2289,8 +2467,9 @@ struct Photonhbt { for (const auto& tg : trueGammas) { const auto it = gammaRecoMap.find(tg.id); if (it != gammaRecoMap.end() && it->second.hasV0 && - (mcIdsWithAnyV0Leg.contains(tg.posId) == 0 || mcIdsWithAnyV0Leg.count(tg.negId) == 0)) + (!mcIdsWithAnyV0Leg.contains(tg.posId) || !mcIdsWithAnyV0Leg.contains(tg.negId))) { ++nBad; + } } fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hStageConsistency"), static_cast(nBad)); } @@ -2302,8 +2481,10 @@ struct Photonhbt { for (const auto& [legId, legFound] : std::initializer_list>{{tg.posId, posFound}, {tg.negId, negFound}}) { - if (legId < 0) + if (legId < 0) { continue; + } + if (legFound) { fRegistryTruthMC.fill(HIST("MC/LegDiag/hLegDRtrue_vs_pt_legFound"), tg.pt, tg.legDRtrue); fRegistryTruthMC.fill(HIST("MC/LegDiag/hLegDEta_legFound_vs_pt"), tg.pt, std::fabs(tg.legDEta)); @@ -2336,24 +2517,27 @@ struct Photonhbt { const float px2 = g2.pt * std::cos(g2.phi), py2 = g2.pt * std::sin(g2.phi); const float kt = 0.5f * std::sqrt((px1 + px2) * (px1 + px2) + (py1 + py2) * (py1 + py2)); - if (mctruth.cfgMCMinKt > 0.f && kt < mctruth.cfgMCMinKt) + if (mctruth.cfgMCMinKt > 0.f && kt < mctruth.cfgMCMinKt) { continue; - if (mctruth.cfgMCMaxKt > 0.f && kt > mctruth.cfgMCMaxKt) + } + if (mctruth.cfgMCMaxKt > 0.f && kt > mctruth.cfgMCMaxKt) { continue; + } const float e1 = g1.pt * std::cosh(g1.eta), e2 = g2.pt * std::cosh(g2.eta); - const float dot = e1 * e2 - (px1 * px2 + py1 * py2 + - g1.pt * std::sinh(g1.eta) * g2.pt * std::sinh(g2.eta)); + const float dot = e1 * e2 - (px1 * px2 + py1 * py2 + g1.pt * std::sinh(g1.eta) * g2.pt * std::sinh(g2.eta)); const float qinv_true = std::sqrt(std::max(0.f, 2.f * dot)); - if (mctruth.cfgMCMaxQinv > 0.f && qinv_true > mctruth.cfgMCMaxQinv) + if (mctruth.cfgMCMaxQinv > 0.f && qinv_true > mctruth.cfgMCMaxQinv) { continue; + } auto it1 = gammaRecoMap.find(g1.id), it2 = gammaRecoMap.find(g2.id); const bool g1Built = (it1 != gammaRecoMap.end()) && it1->second.hasV0; const bool g2Built = (it2 != gammaRecoMap.end()) && it2->second.hasV0; const bool g1Sel = (it1 != gammaRecoMap.end()) && it1->second.passesCut; const bool g2Sel = (it2 != gammaRecoMap.end()) && it2->second.passesCut; + const bool pairClassSel = !cfgDoPhotonClassPairCut.value || (g1Sel && g2Sel && pairutil::isPairPhotonClassSelected(it1->second.legCounts, it2->second.legCounts, mPhotonClassSelA, mPhotonClassSelB)); const bool pairAll4LegsThisColl = legIdsThisCollision.contains(g1.posId) > 0 && legIdsThisCollision.contains(g1.negId) > 0 && @@ -2370,13 +2554,13 @@ struct Photonhbt { fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hQinvVsKt_all4LegsThisColl"), kt, qinv_true); fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hDEtaDPhi_all4LegsThisColl"), deta, dphi); } - if (g1Built && g2Built) { + if (g1Sel && g2Sel) { fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hSparse_DEtaDPhi_qinv_bothPhotonsBuilt"), deta, dphi, qinv_true); fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hSparse_DEtaDPhi_kT_bothPhotonsBuilt"), deta, dphi, kt); fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hQinvVsKt_bothPhotonsBuilt"), kt, qinv_true); fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hDEtaDPhi_bothPhotonsBuilt"), deta, dphi); } - if (g1Sel && g2Sel) { + if (g1Sel && g2Sel && pairClassSel) { fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hSparse_DEtaDPhi_qinv_bothPhotonsSelected"), deta, dphi, qinv_true); fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hSparse_DEtaDPhi_kT_bothPhotonsSelected"), deta, dphi, kt); fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hQinvVsKt_bothPhotonsSelected"), kt, qinv_true); @@ -2385,8 +2569,9 @@ struct Photonhbt { if (g1.rTrue >= 0.f && g2.rTrue >= 0.f) { fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hRconv1_vs_Rconv2_truthConverted"), g1.rTrue, g2.rTrue); - if (g1Sel && g2Sel) + if (g1Sel && g2Sel && pairClassSel) { fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hRconv1_vs_Rconv2_bothPhotonsSelected"), g1.rTrue, g2.rTrue); + } } const float ut = uTrue(g1, g2); fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hUtrueVsKt_truthConverted"), kt, ut); @@ -2399,7 +2584,7 @@ struct Photonhbt { fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hUtrueVsKt_bothPhotonsBuilt"), kt, ut); fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hUtrueVsQinv_bothPhotonsBuilt"), qinv_true, ut); } - if (g1Sel && g2Sel) { + if (g1Sel && g2Sel && pairClassSel) { fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hUtrueVsKt_bothPhotonsSelected"), kt, ut); fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hUtrueVsQinv_bothPhotonsSelected"), qinv_true, ut); } @@ -2408,20 +2593,24 @@ struct Photonhbt { : (g1.rTrue >= 0.f ? g1.rTrue : g2.rTrue); if (minRconv >= 0.f) { fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hMinRconv_vs_kT_truthConverted"), kt, minRconv); - if (g1Sel && g2Sel) + if (g1Sel && g2Sel && pairClassSel) { fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hMinRconv_vs_kT_bothPhotonsSelected"), kt, minRconv); + } } fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hStage_vs_kT"), kt, 0.f); - if (pairAll4LegsThisColl) + if (pairAll4LegsThisColl) { fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hStage_vs_kT"), kt, 1.f); - if (g1Built && g2Built) + } + if (g1Built && g2Built) { fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hStage_vs_kT"), kt, 2.f); - if (g1Sel && g2Sel) + } + if (g1Sel && g2Sel && pairClassSel) { fRegistryTruthMC.fill(HIST("MC/TruthAO2D/hStage_vs_kT"), kt, 3.f); + } const auto itCB = crossBuildMap.find(g1.id); - if (itCB != crossBuildMap.end() && itCB->second.count(g2.id) > 0) { + if (itCB != crossBuildMap.end() && itCB->second.contains(g2.id)) { fRegistryTruthMC.fill(HIST("MC/PairCrossBuild/hSparse_DEtaDPhi_kT"), deta, dphi, kt); fRegistryTruthMC.fill(HIST("MC/PairCrossBuild/hStageOut_vs_kT"), kt, static_cast((g1Built ? 1 : 0) + (g2Built ? 1 : 0))); @@ -2437,13 +2626,15 @@ struct Photonhbt { for (const auto& [tgRef, legId] : std::initializer_list, int>>{ {std::cref(g1), g1.posId}, {std::cref(g1), g1.negId}, {std::cref(g2), g2.posId}, {std::cref(g2), g2.negId}}) { - if (legId < 0 || legIdsThisCollision.contains(legId) > 0) + if (legId < 0 || legIdsThisCollision.contains(legId) > 0) { continue; + } const auto& tg_parent = tgRef.get(); const auto legPtTrue = static_cast(emmcParticles.iteratorAt(legId).pt()); fRegistryTruthMC.fill(HIST("MC/LegDiag/hMissingLegPt_vs_kT"), kt, legPtTrue); - if (tg_parent.rTrue >= 0.f) + if (tg_parent.rTrue >= 0.f) { fRegistryTruthMC.fill(HIST("MC/LegDiag/hMissingLegRconv_vs_kT"), kt, tg_parent.rTrue); + } } } } @@ -2454,8 +2645,9 @@ struct Photonhbt { for (size_t j = i + 1; j < trueGammas.size(); ++j) { const auto& g1 = trueGammas[i]; const auto& g2 = trueGammas[j]; - if (!passAsymmetryCut(g1.pt, g2.pt)) + if (!passAsymmetryCut(g1.pt, g2.pt)) { continue; + } const float deta = g1.eta - g2.eta; const float dphi = wrapPhi(g1.phi - g2.phi); const float px1 = g1.pt * std::cos(g1.phi), py1 = g1.pt * std::sin(g1.phi); @@ -2487,8 +2679,9 @@ struct Photonhbt { for (const auto& poolEvent : truthGammaPool[keyBin]) { for (const auto& g1 : trueGammas) { for (const auto& g2 : poolEvent) { - if (!passAsymmetryCut(g1.pt, g2.pt)) + if (!passAsymmetryCut(g1.pt, g2.pt)) { continue; + } const float deta = g1.eta - g2.eta; const float dphi = wrapPhi(g1.phi - g2.phi); const float px1 = g1.pt * std::cos(g1.phi), py1 = g1.pt * std::sin(g1.phi); @@ -2508,8 +2701,9 @@ struct Photonhbt { if (!trueGammas.empty()) { auto& poolBin = truthGammaPool[keyBin]; poolBin.push_back(trueGammas); - if (static_cast(poolBin.size()) > mctruth.cfgTruthMixDepth.value) + if (static_cast(poolBin.size()) > mctruth.cfgTruthMixDepth.value) { poolBin.pop_front(); + } } } // end cfgDoTruthMix } // end collision loop @@ -2546,7 +2740,7 @@ struct Photonhbt { PROCESS_SWITCH(Photonhbt, processMC, "MC CF + truth efficiency maps for CF correction", false); }; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +WorkflowSpec defineDataProcessing(ConfigContext const& context) { - return WorkflowSpec{adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(context)}; } diff --git a/PWGEM/PhotonMeson/Utils/EMPhoton.h b/PWGEM/PhotonMeson/Utils/EMPhoton.h index cc50f02855e..fb24ab96b52 100644 --- a/PWGEM/PhotonMeson/Utils/EMPhoton.h +++ b/PWGEM/PhotonMeson/Utils/EMPhoton.h @@ -15,41 +15,48 @@ #ifndef PWGEM_PHOTONMESON_UTILS_EMPHOTON_H_ #define PWGEM_PHOTONMESON_UTILS_EMPHOTON_H_ -#include +#include "PWGEM/PhotonMeson/Utils/PairUtilities.h" -#include +#include namespace o2::aod::pwgem::photonmeson::utils { class EMPhoton { public: - EMPhoton(float pt, float eta, float phi, float mass) + EMPhoton(float pt, float eta, float phi, float mass) : fPt(pt), fEta(eta), fPhi(phi), fMass(mass) { - fPt = pt; - fEta = eta; - fPhi = phi; - fMass = mass; } - ~EMPhoton() {} + ~EMPhoton() = default; + + [[nodiscard]] float pt() const { return fPt; } + [[nodiscard]] float eta() const { return fEta; } + [[nodiscard]] float phi() const { return fPhi; } - float pt() const { return fPt; } - float eta() const { return fEta; } - float phi() const { return fPhi; } + [[nodiscard]] float p() const { return fPt * std::cosh(fEta); } + [[nodiscard]] float px() const { return fPt * std::cos(fPhi); } + [[nodiscard]] float py() const { return fPt * std::sin(fPhi); } + [[nodiscard]] float pz() const { return fPt * std::sinh(fEta); } + [[nodiscard]] float mass() const { return fMass; } + [[nodiscard]] float rapidity() const { return std::log((std::sqrt(std::pow(fMass, 2) + std::pow(fPt * std::cosh(fEta), 2)) + fPt * std::sinh(fEta)) / std::sqrt(std::pow(fMass, 2) + std::pow(fPt, 2))); } - float p() const { return fPt * std::cosh(fEta); } - float px() const { return fPt * std::cos(fPhi); } - float py() const { return fPt * std::sin(fPhi); } - float pz() const { return fPt * std::sinh(fEta); } - float mass() const { return fMass; } - float rapidity() const { return std::log((std::sqrt(std::pow(fMass, 2) + std::pow(fPt * std::cosh(fEta), 2)) + fPt * std::sinh(fEta)) / std::sqrt(std::pow(fMass, 2) + std::pow(fPt, 2))); } + // optional leg track-composition info for pair-level photon-class selections (PCM only) + void setLegCounts(o2::aod::pwgem::photonmeson::utils::pairutil::V0PhotonLegCounts const& c) + { + fLegCounts = c; + fHasLegCounts = true; + } + [[nodiscard]] o2::aod::pwgem::photonmeson::utils::pairutil::V0PhotonLegCounts const& legCounts() const { return fLegCounts; } + [[nodiscard]] bool hasLegCounts() const { return fHasLegCounts; } protected: - float fPt; - float fEta; - float fPhi; - float fMass; + float fPt{0.f}; + float fEta{0.f}; + float fPhi{0.f}; + float fMass{0.f}; + o2::aod::pwgem::photonmeson::utils::pairutil::V0PhotonLegCounts fLegCounts{}; + bool fHasLegCounts{false}; }; } // namespace o2::aod::pwgem::photonmeson::utils diff --git a/PWGEM/PhotonMeson/Utils/PairUtilities.h b/PWGEM/PhotonMeson/Utils/PairUtilities.h index 8a1ff20de88..e19b8db1458 100644 --- a/PWGEM/PhotonMeson/Utils/PairUtilities.h +++ b/PWGEM/PhotonMeson/Utils/PairUtilities.h @@ -19,6 +19,7 @@ #include #include +#include #include @@ -37,14 +38,14 @@ namespace o2::aod::pwgem::photonmeson::photonpair { enum PairType { kPCMPCM = 0, - kPHOSPHOS = 1, - kEMCEMC = 2, - kPCMPHOS = 3, - kPCMEMC = 4, - kPCMDalitzEE = 5, - kPCMDalitzMuMu = 6, - kPHOSEMC = 7, - kEEEE = 8, // dielectron-dielectron + kPHOSPHOS, + kEMCEMC, + kPCMPHOS, + kPCMEMC, + kPCMDalitzEE, + kPCMDalitzMuMu, + kPHOSEMC, + kEEEE, // dielectron-dielectron kNpair, }; @@ -69,4 +70,84 @@ bool DoesV0LegMatchWithCluster(TV0Leg const& v0leg, TCluster const& cluster, con } } // namespace o2::aod::pwgem::photonmeson::photonpair +namespace o2::aod::pwgem::photonmeson::utils::pairutil +{ +// ─── photon-class selection by leg track composition ────────────────────── + +struct V0PhotonLegCounts { + int nITSTPC{0}; + int nITSOnly{0}; + int nTPCOnly{0}; + int nTRD{0}; + int nTOF{0}; +}; + +template +V0PhotonLegCounts getV0PhotonLegCounts(TLeg const& pos, TLeg const& ele) +{ + V0PhotonLegCounts c; + auto countLeg = [&c](auto const& l) { + const bool its = l.hasITS(); + const bool tpc = l.hasTPC(); + if (its && tpc) { + c.nITSTPC++; + } else if (its) { + c.nITSOnly++; + } else { + c.nTPCOnly++; + } + if (l.hasTRD()) { + c.nTRD++; + } + if (l.hasTOF()) { + c.nTOF++; + } + }; + countLeg(pos); + countLeg(ele); + return c; +} + +struct V0PhotonClassSelection { + int minITSTPC{0}, maxITSTPC{2}; + int minITSOnly{0}, maxITSOnly{2}; + int minTPCOnly{0}, maxTPCOnly{2}; + int minTRD{0}, maxTRD{2}; + int minTOF{0}, maxTOF{2}; + + [[nodiscard]] bool isSelected(V0PhotonLegCounts const& c) const + { + return c.nITSTPC >= minITSTPC && c.nITSTPC <= maxITSTPC && + c.nITSOnly >= minITSOnly && c.nITSOnly <= maxITSOnly && + c.nTPCOnly >= minTPCOnly && c.nTPCOnly <= maxTPCOnly && + c.nTRD >= minTRD && c.nTRD <= maxTRD && + c.nTOF >= minTOF && c.nTOF <= maxTOF; + } +}; + +inline bool isPairPhotonClassSelected(V0PhotonLegCounts const& c1, V0PhotonLegCounts const& c2, + V0PhotonClassSelection const& selA, V0PhotonClassSelection const& selB) +{ + return (selA.isSelected(c1) && selB.isSelected(c2)) || + (selA.isSelected(c2) && selB.isSelected(c1)); +} + +template +V0PhotonClassSelection buildV0PhotonClassSelection(TGroup const& g) +{ + V0PhotonClassSelection s; + s.minITSTPC = g.cfgMinNLegsITSTPC.value; + s.maxITSTPC = g.cfgMaxNLegsITSTPC.value; + s.minITSOnly = g.cfgMinNLegsITSOnly.value; + s.maxITSOnly = g.cfgMaxNLegsITSOnly.value; + s.minTPCOnly = g.cfgMinNLegsTPCOnly.value; + s.maxTPCOnly = g.cfgMaxNLegsTPCOnly.value; + s.minTRD = g.cfgMinNLegsTRD.value; + s.maxTRD = g.cfgMaxNLegsTRD.value; + s.minTOF = g.cfgMinNLegsTOF.value; + s.maxTOF = g.cfgMaxNLegsTOF.value; + return s; +} +} // namespace o2::aod::pwgem::photonmeson::utils::pairutil + #endif // PWGEM_PHOTONMESON_UTILS_PAIRUTILITIES_H_