[PWGLF] Add EventQualityStudy in hStrangeCorrelation.cxx

PWGLF/Tasks/Strangeness/hStrangeCorrelation.cxx

@@ -330,6 +330,9 @@
     // trimmed by skipUnderOverflowInTHn, so what you configure is what you get.
     ConfigurableAxis axisGenStudyNch{"axisGenStudyNch", {VARIABLE_WIDTH, 0.0f, 2.0f, 5.0f, 10.0f, 15.0f, 20.0f, 25.0f, 30.0f, 40.0f, 60.0f, 100.0f}, "generated charged multiplicity in |#eta| < 0.8"};
     Configurable<bool> doClosureTestStages{"doClosureTestStages", true, "create and fill the whole ClosureTest/PairLossK0 folder: the truth and any-reconstructed-object stages of the truth h-K0 pair, mirroring the first processPairLossK0MC stages"};
+    Configurable<bool> doEventQualityStudy{"doEventQualityStudy", false, "add eight EventQuality folders to the closure PairLossK0 study; requires doClosureTestStages"};
+    Configurable<int> eventQualityDeltaPhiBins{"eventQualityDeltaPhiBins", 36, "delta phi bins in the sparse event-quality histograms"};
+    Configurable<int> eventQualityDeltaEtaBins{"eventQualityDeltaEtaBins", 20, "delta eta bins in the sparse event-quality histograms"};
     Configurable<bool> applyRecoEventSelection{"applyRecoEventSelection", true, "apply the standard reconstructed-event selection in the K0 pair-loss diagnostic"};
     Configurable<bool> fillFinalPairOnce{"fillFinalPairOnce", true, "part 1: fill each truth h-K0 pair at most once at the final stage. Set false to loop over all reconstructed trigger x K0 combinations, as the data path does"};
     Configurable<bool> finalPairUseBestCollisionOnly{"finalPairUseBestCollisionOnly", true, "part 1: build the final trigger and K0 objects only in the best collision. Set false to use every associated reconstructed collision that passes the event selection, pairing within one collision, as the data path does"};
@@ -532,6 +535,23 @@
     "V0 final selection",
     "Final reconstructed pair"};
 
+  // Independent event classifications, not cumulative cuts. Values are fixed
+  // diagnostic bins (labels below), and never change the parent pair selection.
+  static constexpr int PairLossNEventQualityGroups = 8;
+  static constexpr std::array<std::string_view, PairLossNEventQualityGroups> PairLossEventQualityNames = {
+    "NContributors", "PVMaxPull", "PVPurity", "CollisionAssociation",
+    "TimeResolution", "Occupancy", "GenMultiplicity", "EventShape"};
+  static constexpr std::array<std::array<std::string_view, 6>, PairLossNEventQualityGroups> PairLossEventQualityLabels = {{
+    {"no contributors", "1-2", "3-5", "6-10", "11-20", ">=21"},
+    {"invalid covariance/pull", "max |pull| <1", "1<=max |pull|<3", "3<=max |pull|<5", "max |pull|>=5", "unused"},
+    {"no PV tracks", "missing MC labels", "purity <0.95", "0.95<=purity<1", "purity =1", "unused"},
+    {"no matched tracks", "one rec collision", "multiple: best fraction >=0.9", "multiple: 0.5<=fraction<0.9", "multiple: fraction <0.5", "unused"},
+    {"invalid time resolution", "0<sigma_t<25 ns", "25<=sigma_t<100 ns", "100<=sigma_t<500 ns", "sigma_t>=500 ns", "unused"},
+    {"unavailable", "0", "1-99", "100-499", "500-999", ">=1000"},
+    {"0", "1-4", "5-9", "10-19", "20-39", ">=40"},
+    {"fewer than 3 shape tracks", "S_T<0.3", "0.3<=S_T<0.7", "S_T>=0.7", "unused", "unused"},
+  }};
+
   struct PairLossTrackInfo {
     int64_t globalIndex = -1;
     float pt = 0.0f;
@@ -742,7 +762,7 @@
   void collectDescendantIds(TMcParticle const& mcParticle, std::vector<int64_t>& ids, int depth = 0)
   {
     ids.push_back(mcParticle.globalIndex());
     if (depth >= 3 || !mcParticle.has_daughters()) {
       return;
     }
     for (auto const& daughter : mcParticle.template daughters_as<aod::McParticles>()) {
@@ -1499,13 +1519,13 @@
                 ((-massWindowConfigurations.maxBgNSigma < assocCandidate.invMassNSigma(Index) && assocCandidate.invMassNSigma(Index) < -massWindowConfigurations.minBgNSigma) ||
                  (-massWindowConfigurations.maxPeakNSigma < assocCandidate.invMassNSigma(Index) && assocCandidate.invMassNSigma(Index) < +massWindowConfigurations.maxPeakNSigma) ||
                  (+massWindowConfigurations.minBgNSigma < assocCandidate.invMassNSigma(Index) && assocCandidate.invMassNSigma(Index) < +massWindowConfigurations.maxBgNSigma))) {
               if (std::abs(deltaphi) < 0.5) {
                 histos.fill(HIST("sameEvent/InvariantMass/") + HIST(V0names[Index]) + HIST("/hNearSide"), ptassoc, pttrigger, getV0InvariantMass<Index>(assoc));
               }
               if (std::abs(PI - deltaphi) < 0.5) {
                 histos.fill(HIST("sameEvent/InvariantMass/") + HIST(V0names[Index]) + HIST("/hAwaySide"), ptassoc, pttrigger, getV0InvariantMass<Index>(assoc));
               }
               if (deltaphi > 1.0 && deltaphi < 1.5) {
                 histos.fill(HIST("sameEvent/InvariantMass/") + HIST(V0names[Index]) + HIST("/hUnderlyingEvent"), ptassoc, pttrigger, getV0InvariantMass<Index>(assoc));
               }
             }
@@ -1756,9 +1776,9 @@
     int binnumb = binMult * nBinsVtxZ + binVtxZ;
     int hastirgorassoc = masterConfigurations.collisionHasTriggOrAssoc;
     if ((hastirgorassoc == 1 && currentCollision.trigParticles.empty()) ||
         (hastirgorassoc == 2 && currentCollision.assocParticles.empty()) ||
         (hastirgorassoc == 3 && currentCollision.trigParticles.empty() && currentCollision.assocParticles.empty()) ||
         (hastirgorassoc == 4 && (currentCollision.trigParticles.empty() || currentCollision.assocParticles.empty()))) {
       return;
     }
     for (const auto& collision : validCollisions[binnumb]) {
@@ -1795,7 +1815,7 @@
             if (Index == assoc.type && assoc.region == 1) {
               fillCorrelationHistogram(histos.get<THn>(HIST("mixedEvent/Signal/") + HIST(V0names[Index])), binFillThn, 1, efficiencyTrigg * efficiencyAssoc, totalEffUncert, 1., 0.);
             }
             if (Index == assoc.type && assoc.region == 2) {
               fillCorrelationHistogram(histos.get<THn>(HIST("mixedEvent/RightBg/") + HIST(V0names[Index])), binFillThn, 1, efficiencyTrigg * efficiencyAssoc, totalEffUncert, 1., 0.);
             }
           });
@@ -2094,8 +2114,8 @@
     int binnumb = binMult * nBinsVtxZ + binVtxZ;
     int hastirgorassoc = masterConfigurations.collisionHasTriggOrAssoc;
     if ((hastirgorassoc == 1 && currentCollision.trigParticles.empty()) ||
         (hastirgorassoc == 2 && currentCollision.assocParticles.empty()) ||
         (hastirgorassoc == 3 && currentCollision.trigParticles.empty() && currentCollision.assocParticles.empty()) ||
         (hastirgorassoc == 4 && (currentCollision.trigParticles.empty() || currentCollision.assocParticles.empty()))) {
       return;
     }
@@ -3119,6 +3139,32 @@
         histos.add("ClosureTest/PairLossK0/Final/sameEvent/K0Short", "truth h-K0 pairs whose trigger and K0 both have a fully selected reconstructed counterpart in the same reconstructed collision", kTHnF, {axisDeltaPhiNDim, axisDeltaEtaNDim, axisPtAssocNDim, axisPtTriggerNDim, axisVtxZNDim, axisMultNDim});
         histos.add("ClosureTest/PairLossK0/Final/hTrigger", "truth triggers with a fully selected reconstructed counterpart;#it{p}_{T}^{truth} (GeV/#it{c});#eta^{truth};#varphi^{truth}", kTH3F, {axesConfigurations.axisPtQA, axesConfigurations.axisEta, axesConfigurations.axisPhi});
         histos.add("ClosureTest/PairLossK0/Final/hK0Short", "truth K0s with a fully selected reconstructed counterpart;#it{p}_{T}^{truth} (GeV/#it{c});#eta^{truth};#varphi^{truth}", kTH3F, {axesConfigurations.axisPtQA, axesConfigurations.axisEta, axesConfigurations.axisPhi});
+        if (pairLossK0Configurations.doEventQualityStudy) {
+          if (pairLossK0Configurations.eventQualityDeltaPhiBins <= 0 || pairLossK0Configurations.eventQualityDeltaEtaBins <= 0) {
+            LOGF(fatal, "Event-quality angular bin counts must be positive");
+          }
+          // Six sparse axes; omit vertex and multiplicity dimensions. Each folder
+          // has exactly three histograms, including its event/trigger denominators.
+          const std::vector<AxisSpec> qualityAxes = {
+            {6, -0.5, 5.5, "event class"},
+            {3, -0.5, 2.5, "entry kind"},
+            {pairLossK0Configurations.eventQualityDeltaPhiBins, axisRanges[0][0], axisRanges[0][1], "#Delta#varphi^{truth}"},
+            {pairLossK0Configurations.eventQualityDeltaEtaBins, axisRanges[1][0], axisRanges[1][1], "#Delta#eta^{truth}"},
+            axisPtAssocNDim,
+            axisPtTriggerNDim,
+          };
+          for (int group = 0; group < PairLossNEventQualityGroups; ++group) {
+            for (auto const& stage : {"Truth", "AnyTrackBoth", "Final"}) {
+              auto histogram = histos.add<THnSparse>(fmt::format("ClosureTest/PairLossK0/EventQuality/{}/{}", PairLossEventQualityNames[group], stage), "unit-weight counts; select entry kind before projection", kTHnSparseF, qualityAxes);
+              for (size_t bin = 0; bin < PairLossEventQualityLabels[group].size(); ++bin) {
+                histogram->GetAxis(0)->SetBinLabel(bin + 1, PairLossEventQualityLabels[group][bin].data());
+              }
+              histogram->GetAxis(1)->SetBinLabel(1, "event (common baseline)");
+              histogram->GetAxis(1)->SetBinLabel(2, "trigger");
+              histogram->GetAxis(1)->SetBinLabel(3, "pair");
+            }
+          }
+        }
       }
       for (int i = 0; i < AssocParticleTypes; i++) {
         if (TESTBIT(doCorrelation, i)) {
@@ -5621,7 +5667,7 @@
     }
   }
 
-  void processClosureTest(aod::McCollision const& /*mcCollision*/,
+  void processClosureTest(aod::McCollision const& mcCollision,
                           soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSels, aod::CentFT0Ms, aod::CentFT0Cs, aod::PVMults>> const& recCollisions,
                           aod::McParticles const& mcParticles,
                           aod::V0MCCores const& v0MCCores,
@@ -5838,6 +5884,116 @@
 
         const float pairLossBestCollisionVtxZ = collision.posZ();
         const float pairLossBestCollisionMultiplicity = masterConfigurations.doPPAnalysis ? collision.centFT0M() : collision.centFT0C();
+        std::array<int, PairLossNEventQualityGroups> eventQualityClasses{};
+        if (pairLossK0Configurations.doEventQualityStudy) {
+          // Category edges correspond to the labels declared above.
+          auto classFromEdges = [](auto value, auto const& edges) {
+            return static_cast<int>(std::upper_bound(edges.begin(), edges.end(), value) - edges.begin());
+          };
+          constexpr double PurePVThreshold = 0.95;
+          constexpr double DominantCollisionFraction = 0.9;
+          constexpr double MajorityCollisionFraction = 0.5;
+          constexpr float GenEtaMax = 0.8f;
+          constexpr float ShapePtMin = 0.15f;
+          constexpr int MinShapeTracks = 3;
+          const auto nContributors = collision.numContrib();
+          eventQualityClasses[0] = classFromEdges(nContributors, std::array{1u, 3u, 6u, 11u, 21u});
+
+          // Maximum component pull, not a 3D chi-square: covariance correlations
+          // are deliberately not used. Invalid covariance gets its own category.
+          const std::array<double, 3> residuals = {collision.posX() - mcCollision.posX(), collision.posY() - mcCollision.posY(), collision.posZ() - mcCollision.posZ()};
+          const std::array<double, 3> variances = {collision.covXX(), collision.covYY(), collision.covZZ()};
+          double maxPull = 0.;
+          bool validPull = true;
+          for (size_t component = 0; component < variances.size(); ++component) {
+            if (!std::isfinite(variances[component]) || variances[component] <= 0. || !std::isfinite(residuals[component])) {
+              validPull = false;
+              break;
+            }
+            maxPull = std::max(maxPull, std::abs(residuals[component]) / std::sqrt(variances[component]));
+          }
+          eventQualityClasses[1] = !validPull || !std::isfinite(maxPull) ? 0 : 1 + classFromEdges(maxPull, std::array{1., 3., 5.});
+
+          std::unordered_set<int64_t> currentMcIds;
+          int nAnyMatched = 0;
+          int nGenCharged = 0;
+          int nShapeTracks = 0;
+          double sumPt = 0., q2x = 0., q2y = 0.;
+          for (auto const& particle : mcParticles) {
+            currentMcIds.insert(particle.globalIndex());
+            nAnyMatched += pairLossAnyTrackMcParticleIds.count(particle.globalIndex()) != 0;
+            if (!particle.isPhysicalPrimary() || !std::isfinite(particle.eta()) || std::abs(particle.eta()) >= GenEtaMax) {
+              continue;
+            }
+            auto const* pdgParticle = pdgDB->GetParticle(particle.pdgCode());
+            if (pdgParticle == nullptr || pdgParticle->Charge() == 0.) {
+              continue;
+            }
+            ++nGenCharged;
+            // Linearized transverse sphericity, using truth primary charged
+            // particles with |eta|<0.8 and pT>=0.15 GeV/c (not spherocity).
+            if (std::isfinite(particle.pt()) && particle.pt() >= ShapePtMin && std::isfinite(particle.phi())) {
+              ++nShapeTracks;
+              sumPt += particle.pt();
+              q2x += particle.pt() * std::cos(2. * particle.phi());
+              q2y += particle.pt() * std::sin(2. * particle.phi());
+            }
+          }
+          int nPVTracks = 0, nLabeledPVTracks = 0, nSameMcPVTracks = 0;
+          std::unordered_set<int64_t> bestMatchedMcIds;
+          const auto bestTracks = tracks.sliceBy(pairLossTracksPerCollision, collision.globalIndex());
+          for (auto const& track : bestTracks) {
+            const bool sameMc = track.has_mcParticle() && currentMcIds.count(track.mcParticleId()) != 0;
+            if (sameMc) {
+              bestMatchedMcIds.insert(track.mcParticleId());
+            }
+            if (track.isPVContributor()) {
+              ++nPVTracks;
+              nLabeledPVTracks += track.has_mcParticle();
+              nSameMcPVTracks += sameMc;
+            }
+          }
+          // Missing labels take precedence; purity otherwise uses all PV tracks.
+          const double purity = nPVTracks > 0 ? static_cast<double>(nSameMcPVTracks) / nPVTracks : 0.;
+          eventQualityClasses[2] = nPVTracks == 0 ? 0 : nLabeledPVTracks < nPVTracks ? 1
+                                                      : purity < PurePVThreshold     ? 2
+                                                      : nSameMcPVTracks < nPVTracks  ? 3
+                                                                                     : 4;
+          const double bestFraction = nAnyMatched > 0 ? static_cast<double>(bestMatchedMcIds.size()) / nAnyMatched : 0.;
+          eventQualityClasses[3] = nAnyMatched == 0 ? 0 : recCollisions.size() == 1                 ? 1
+                                                        : bestFraction >= DominantCollisionFraction ? 2
+                                                        : bestFraction >= MajorityCollisionFraction ? 3
+                                                                                                    : 4;
+          const double timeResolution = collision.collisionTimeRes(); // ns
+          eventQualityClasses[4] = !std::isfinite(timeResolution) || timeResolution <= 0. ? 0 : 1 + classFromEdges(timeResolution, std::array{25., 100., 500.});
+          const auto occupancy = collision.trackOccupancyInTimeRange();
+          eventQualityClasses[5] = occupancy < 0 ? 0 : 1 + classFromEdges(occupancy, std::array{1, 100, 500, 1000});
+          eventQualityClasses[6] = classFromEdges(nGenCharged, std::array{1, 5, 10, 20, 40});
+          if (nShapeTracks >= MinShapeTracks && sumPt > 0.) {
+            const double sphericity = std::clamp(1. - std::hypot(q2x, q2y) / sumPt, 0., 1.);
+            eventQualityClasses[7] = 1 + classFromEdges(sphericity, std::array{0.3, 0.7});
+          }
+        }
+        auto fillEventQuality = [&](int entryKind, float deltaPhi, float deltaEta, float assocPt, float triggerPt, bool hasAny, bool hasFinal) {
+          if (!pairLossK0Configurations.doEventQualityStudy) {
+            return;
+          }
+          static_for<0, PairLossNEventQualityGroups - 1>([&](auto i) {
+            constexpr int Index = i.value;
+            histos.fill(HIST("ClosureTest/PairLossK0/EventQuality/") + HIST(PairLossEventQualityNames[Index]) + HIST("/Truth"), eventQualityClasses[Index], entryKind, deltaPhi, deltaEta, assocPt, triggerPt);
+            if (hasAny) {
+              histos.fill(HIST("ClosureTest/PairLossK0/EventQuality/") + HIST(PairLossEventQualityNames[Index]) + HIST("/AnyTrackBoth"), eventQualityClasses[Index], entryKind, deltaPhi, deltaEta, assocPt, triggerPt);
+            }
+            if (hasFinal) {
+              histos.fill(HIST("ClosureTest/PairLossK0/EventQuality/") + HIST(PairLossEventQualityNames[Index]) + HIST("/Final"), eventQualityClasses[Index], entryKind, deltaPhi, deltaEta, assocPt, triggerPt);
+            }
+          });
+        };
+        // Unused coordinates lie inside regular bins. Select entry kind before
+        // projecting; reset pair-coordinate ranges for event/trigger counts.
+        const float countDeltaPhi = (axisRanges[0][0] + axisRanges[0][1]) * 0.5f;
+        const float countDeltaEta = (axisRanges[1][0] + axisRanges[1][1]) * 0.5f;
+        fillEventQuality(0, countDeltaPhi, countDeltaEta, axisRanges[2][0], axisRanges[3][0], true, true);
         std::vector<PairLossTruthTrackInfo> pairLossTruthTriggers;
         std::vector<PairLossTruthK0Info> pairLossTruthK0s;
 
@@ -5902,6 +6058,9 @@
         // at truth level in Truth/ and AnyTrack/, at any level in AnyTrackK0/ and
         // AnyTrackBoth/, at fully-selected level in Final/.
         for (auto const& truthTrigger : pairLossTruthTriggers) {
+          if (pairLossK0Configurations.doEventQualityStudy) {
+            fillEventQuality(1, countDeltaPhi, countDeltaEta, axisRanges[2][0], truthTrigger.pt, pairLossAnyTrackMcParticleIds.count(truthTrigger.globalIndex) != 0, pairLossHasFinalTrigger(truthTrigger.globalIndex));
+          }
           histos.fill(HIST("ClosureTest/PairLossK0/Truth/hTrigger"), truthTrigger.pt, truthTrigger.eta, truthTrigger.phi);
           histos.fill(HIST("ClosureTest/PairLossK0/AnyTrackK0/hTrigger"), truthTrigger.pt, truthTrigger.eta, truthTrigger.phi);
           if (pairLossAnyTrackMcParticleIds.find(truthTrigger.globalIndex) != pairLossAnyTrackMcParticleIds.end()) {
@@ -5952,7 +6111,9 @@
             if (triggerHasAnyTrack && k0HasAnyV0) {
               histos.fill(HIST("ClosureTest/PairLossK0/AnyTrackBoth/sameEvent/K0Short"), truthDeltaPhi, truthDeltaEta, truthK0.pt, truthTrigger.pt, pairLossBestCollisionVtxZ, pairLossBestCollisionMultiplicity);
             }
-            if (pairLossHasFinalPair(truthTrigger.globalIndex, truthK0.globalIndex)) {
+            const bool hasFinalPair = pairLossHasFinalPair(truthTrigger.globalIndex, truthK0.globalIndex);
+            fillEventQuality(2, truthDeltaPhi, truthDeltaEta, truthK0.pt, truthTrigger.pt, triggerHasAnyTrack && k0HasAnyV0, hasFinalPair);
+            if (hasFinalPair) {
               histos.fill(HIST("ClosureTest/PairLossK0/Final/sameEvent/K0Short"), truthDeltaPhi, truthDeltaEta, truthK0.pt, truthTrigger.pt, pairLossBestCollisionVtxZ, pairLossBestCollisionMultiplicity);
             }
           }