[PWGJE] Fix pTHat handling, remove HepMCXSections dependency, add R-dependent matching

- trackEfficiency: Remove HepMCXSections and JMcCollisionPIs table
  subscriptions from all MC processes; use ptHard() from JMcCollisions
  with weight-based fallback. Add applyRCTSelections configurable.
- jetSpectraCharged: Pass pTHat as parameter to fill functions instead
  of computing internally; fix pTHatMaxMCD->pTHatMaxMCP bug in
  fillGeoMatchedAreaSubHistograms; change MCD weighted process
  subscriptions to JetCollisionsMCD for mcCollision() access.
- jetMatchingMC/JetMatchingUtilities: Replace single maxMatchingDistance
  with per-R configurables (jetRadiiForMatchingDistance +
  maxMatchingDistancePerJetR) for R-dependent geometric matching.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: joonsukbae

PWGJE/Core/JetMatchingUtilities.h

@@ -293,7 +293,7 @@ std::tuple<std::vector<int>, std::vector<int>> MatchJetsGeometrically(
 }
 
 template <typename T, typename U>
-void MatchGeo(T const& jetsBasePerCollision, U const& jetsTagPerCollision, std::vector<std::vector<int>>& baseToTagMatchingGeo, std::vector<std::vector<int>>& tagToBaseMatchingGeo, float maxMatchingDistance)
+void MatchGeo(T const& jetsBasePerCollision, U const& jetsTagPerCollision, std::vector<std::vector<int>>& baseToTagMatchingGeo, std::vector<std::vector<int>>& tagToBaseMatchingGeo, float maxMatchingDistance, std::vector<double> const& jetRadiiForMatchingDistance = {}, std::vector<double> const& maxMatchingDistancePerJetR = {})
 {
   std::vector<double> jetsR;
   for (const auto& jetBase : jetsBasePerCollision) {
@@ -332,7 +332,14 @@ void MatchGeo(T const& jetsBasePerCollision, U const& jetsTagPerCollision, std::
       jetsTagEta.emplace_back(jetTag.eta());
       jetsTagGlobalIndex.emplace_back(jetTag.globalIndex());
     }
-    std::tie(baseToTagMatchingGeoIndex, tagToBaseMatchingGeoIndex) = MatchJetsGeometrically(jetsBasePhi, jetsBaseEta, jetsTagPhi, jetsTagEta, maxMatchingDistance); // change max distnace to a function call
+    float effectiveMatchingDistance = maxMatchingDistance;
+    for (std::size_t i = 0; i < jetRadiiForMatchingDistance.size() && i < maxMatchingDistancePerJetR.size(); i++) {
+      if (std::round(jetRadiiForMatchingDistance[i] * 100.0) == std::round(jetR)) {
+        effectiveMatchingDistance = maxMatchingDistancePerJetR[i];
+        break;
+      }
+    }
+    std::tie(baseToTagMatchingGeoIndex, tagToBaseMatchingGeoIndex) = MatchJetsGeometrically(jetsBasePhi, jetsBaseEta, jetsTagPhi, jetsTagEta, effectiveMatchingDistance);
     int jetBaseIndex = 0;
     int jetTagIndex = 0;
     for (const auto& jetBase : jetsBasePerCollision) {
@@ -564,11 +571,11 @@ void MatchPt(T const& jetsBasePerCollision, U const& jetsTagPerCollision, std::v
 
 // function that calls all the Match functions
 template <bool jetsBaseIsMc, bool jetsTagIsMc, typename T, typename U, typename V, typename M, typename N, typename O, typename P, typename R>
-void doAllMatching(T const& jetsBasePerCollision, U const& jetsTagPerCollision, std::vector<std::vector<int>>& baseToTagMatchingGeo, std::vector<std::vector<int>>& baseToTagMatchingPt, std::vector<std::vector<int>>& baseToTagMatchingHF, std::vector<std::vector<int>>& tagToBaseMatchingGeo, std::vector<std::vector<int>>& tagToBaseMatchingPt, std::vector<std::vector<int>>& tagToBaseMatchingHF, V const& candidatesBase, M const& tracksBase, N const& clustersBase, O const& candidatesTag, P const& tracksTag, R const& clustersTag, bool doMatchingGeo, bool doMatchingHf, bool doMatchingPt, float maxMatchingDistance, float minPtFraction)
+void doAllMatching(T const& jetsBasePerCollision, U const& jetsTagPerCollision, std::vector<std::vector<int>>& baseToTagMatchingGeo, std::vector<std::vector<int>>& baseToTagMatchingPt, std::vector<std::vector<int>>& baseToTagMatchingHF, std::vector<std::vector<int>>& tagToBaseMatchingGeo, std::vector<std::vector<int>>& tagToBaseMatchingPt, std::vector<std::vector<int>>& tagToBaseMatchingHF, V const& candidatesBase, M const& tracksBase, N const& clustersBase, O const& candidatesTag, P const& tracksTag, R const& clustersTag, bool doMatchingGeo, bool doMatchingHf, bool doMatchingPt, float maxMatchingDistance, float minPtFraction, std::vector<double> const& jetRadiiForMatchingDistance = {}, std::vector<double> const& maxMatchingDistancePerJetR = {})
 {
   // geometric matching
   if (doMatchingGeo) {
-    MatchGeo(jetsBasePerCollision, jetsTagPerCollision, baseToTagMatchingGeo, tagToBaseMatchingGeo, maxMatchingDistance);
+    MatchGeo(jetsBasePerCollision, jetsTagPerCollision, baseToTagMatchingGeo, tagToBaseMatchingGeo, maxMatchingDistance, jetRadiiForMatchingDistance, maxMatchingDistancePerJetR);
   }
   // pt matching
   if (doMatchingPt) {

PWGJE/TableProducer/Matching/jetMatchingMC.h

@@ -37,7 +37,8 @@ struct JetMatchingMc {
   o2::framework::Configurable<bool> doMatchingGeo{"doMatchingGeo", true, "Enable geometric matching"};
   o2::framework::Configurable<bool> doMatchingPt{"doMatchingPt", true, "Enable pt matching"};
   o2::framework::Configurable<bool> doMatchingHf{"doMatchingHf", false, "Enable HF matching"};
-  o2::framework::Configurable<float> maxMatchingDistance{"maxMatchingDistance", 0.24f, "Max matching distance"};
+  o2::framework::Configurable<std::vector<double>> jetRadiiForMatchingDistance{"jetRadiiForMatchingDistance", {0.4}, "Jet R values for per-R matching distance, e.g. {0.2, 0.4, 0.6}"};
+  o2::framework::Configurable<std::vector<double>> maxMatchingDistancePerJetR{"maxMatchingDistancePerJetR", {0.24}, "Max matching distance for each R in jetRadiiForMatchingDistance, e.g. {0.12, 0.24, 0.36}"};
   o2::framework::Configurable<float> minPtFraction{"minPtFraction", 0.5f, "Minimum pt fraction for pt matching"};
 
   o2::framework::Produces<JetsBasetoTagMatchingTable> jetsBasetoTagMatchingTable;
@@ -54,6 +55,9 @@ struct JetMatchingMc {
 
   void init(o2::framework::InitContext const&)
   {
+    if (jetRadiiForMatchingDistance->size() != maxMatchingDistancePerJetR->size()) {
+      LOGP(fatal, "jetRadiiForMatchingDistance and maxMatchingDistancePerJetR must have the same number of entries");
+    }
   }
 
   void processJets(o2::aod::JetMcCollisions const& mcCollisions, o2::aod::JetCollisionsMCD const& collisions,
@@ -84,7 +88,7 @@ struct JetMatchingMc {
         const auto jetsBasePerColl = jetsBase.sliceBy(baseJetsPerCollision, jetsBaseIsMc ? mcCollision.globalIndex() : collision.globalIndex());
         const auto jetsTagPerColl = jetsTag.sliceBy(tagJetsPerCollision, jetsTagIsMc ? mcCollision.globalIndex() : collision.globalIndex());
 
-        jetmatchingutilities::doAllMatching<jetsBaseIsMc, jetsTagIsMc>(jetsBasePerColl, jetsTagPerColl, jetsBasetoTagMatchingGeo, jetsBasetoTagMatchingPt, jetsBasetoTagMatchingHF, jetsTagtoBaseMatchingGeo, jetsTagtoBaseMatchingPt, jetsTagtoBaseMatchingHF, candidatesBase, tracks, clusters, candidatesTag, particles, particles, doMatchingGeo, doMatchingHf, doMatchingPt, maxMatchingDistance, minPtFraction);
+        jetmatchingutilities::doAllMatching<jetsBaseIsMc, jetsTagIsMc>(jetsBasePerColl, jetsTagPerColl, jetsBasetoTagMatchingGeo, jetsBasetoTagMatchingPt, jetsBasetoTagMatchingHF, jetsTagtoBaseMatchingGeo, jetsTagtoBaseMatchingPt, jetsTagtoBaseMatchingHF, candidatesBase, tracks, clusters, candidatesTag, particles, particles, doMatchingGeo, doMatchingHf, doMatchingPt, 0.24f, minPtFraction, jetRadiiForMatchingDistance, maxMatchingDistancePerJetR);
       }
     }
     for (auto i = 0; i < jetsBase.size(); ++i) {

PWGJE/Tasks/jetSpectraCharged.cxx

@@ -463,9 +463,8 @@ struct JetSpectraCharged {
   }
 
   template <typename TJets>
-  void fillJetHistograms(TJets const& jet, float centrality, float weight = 1.0)
+  void fillJetHistograms(TJets const& jet, float centrality, float weight = 1.0, float pTHat = 999.0)
   {
-    float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent));
     if (jet.pt() > pTHatMaxMCD * pTHat || pTHat < pTHatAbsoluteMin) {
       return;
     }
@@ -487,9 +486,8 @@ struct JetSpectraCharged {
   }
 
   template <typename TJets>
-  void fillJetAreaSubHistograms(TJets const& jet, float centrality, float rho, float weight = 1.0)
+  void fillJetAreaSubHistograms(TJets const& jet, float centrality, float rho, float weight = 1.0, float pTHat = 999.0)
   {
-    float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent));
     if (jet.pt() > pTHatMaxMCD * pTHat || pTHat < pTHatAbsoluteMin) {
       return;
     }
@@ -516,9 +514,8 @@ struct JetSpectraCharged {
   }
 
   template <typename TJets>
-  void fillMCPHistograms(TJets const& jet, float weight = 1.0)
+  void fillMCPHistograms(TJets const& jet, float weight = 1.0, float pTHat = 999.0)
   {
-    float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent));
     if (jet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) {
       return;
     }
@@ -538,9 +535,8 @@ struct JetSpectraCharged {
   }
 
   template <typename TJets>
-  void fillMCPAreaSubHistograms(TJets const& jet, float rho = 0.0, float weight = 1.0)
+  void fillMCPAreaSubHistograms(TJets const& jet, float rho = 0.0, float weight = 1.0, float pTHat = 999.0)
   {
-    float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent));
     if (jet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) {
       return;
     }
@@ -559,9 +555,8 @@ struct JetSpectraCharged {
   }
 
   template <typename TJets>
-  void fillEventWiseConstituentSubtractedHistograms(TJets const& jet, float centrality, float weight = 1.0)
+  void fillEventWiseConstituentSubtractedHistograms(TJets const& jet, float centrality, float weight = 1.0, float pTHat = 999.0)
   {
-    float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent));
     if (jet.pt() > pTHatMaxMCD * pTHat || pTHat < pTHatAbsoluteMin) {
       return;
     }
@@ -572,9 +567,8 @@ struct JetSpectraCharged {
   }
 
   template <typename TBase, typename TTag>
-  void fillMatchedHistograms(TBase const& jetMCD, float weight = 1.0)
+  void fillMatchedHistograms(TBase const& jetMCD, float weight = 1.0, float pTHat = 999.0)
   {
-    float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent));
     if (jetMCD.pt() > pTHatMaxMCD * pTHat || pTHat < pTHatAbsoluteMin) {
       return;
     }
@@ -665,15 +659,14 @@ struct JetSpectraCharged {
   }
 
   template <typename TBase, typename TTag>
-  void fillGeoMatchedAreaSubHistograms(TBase const& jetMCD, float rho, float mcrho = 0.0, float weight = 1.0)
+  void fillGeoMatchedAreaSubHistograms(TBase const& jetMCD, float rho, float mcrho = 0.0, float weight = 1.0, float pTHat = 999.0)
   {
-    float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent));
     if (jetMCD.pt() > pTHatMaxMCD * pTHat) {
       return;
     }
     if (jetMCD.has_matchedJetGeo()) {
       for (const auto& jetMCP : jetMCD.template matchedJetGeo_as<std::decay_t<TTag>>()) {
-        if (jetMCP.pt() > pTHatMaxMCD * pTHat) {
+        if (jetMCP.pt() > pTHatMaxMCP * pTHat) {
           continue;
         }
         if (jetMCD.r() == round(selectedJetsRadius * 100.0f)) {
@@ -848,17 +841,18 @@ struct JetSpectraCharged {
   }
   PROCESS_SWITCH(JetSpectraCharged, processSpectraAreaSubMCD, "jet spectra with rho-area subtraction for MCD", false);
 
-  void processSpectraMCDWeighted(soa::Filtered<aod::JetCollisions>::iterator const& collision,
+  void processSpectraMCDWeighted(soa::Filtered<aod::JetCollisionsMCD>::iterator const& collision,
                                  soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents> const& jets,
-                                 aod::JetTracks const&)
+                                 aod::JetTracks const&,
+                                 aod::JetMcCollisions const&)
   {
     bool fillHistograms = false;
     bool isWeighted = true;
     float eventWeight = collision.weight();
     if (!applyCollisionCuts(collision, fillHistograms, isWeighted, eventWeight)) {
       return;
     }
-
+    float pTHat = collision.has_mcCollision() && collision.mcCollision().ptHard() < 999.0f ? collision.mcCollision().ptHard() : 10. / (std::pow(eventWeight, 1.0 / pTHatExponent));
     float centrality = -1.0;
     checkCentFT0M ? centrality = collision.centFT0M() : centrality = collision.centFT0C();
 
@@ -869,22 +863,23 @@ struct JetSpectraCharged {
       if (!isAcceptedJet<aod::JetTracks>(jet)) {
         continue;
       }
-      fillJetHistograms(jet, centrality, eventWeight);
+      fillJetHistograms(jet, centrality, eventWeight, pTHat);
     }
   }
   PROCESS_SWITCH(JetSpectraCharged, processSpectraMCDWeighted, "jet finder QA mcd with weighted events", false);
 
-  void processSpectraAreaSubMCDWeighted(soa::Filtered<soa::Join<aod::JetCollisions, aod::BkgChargedRhos>>::iterator const& collision,
+  void processSpectraAreaSubMCDWeighted(soa::Filtered<soa::Join<aod::JetCollisionsMCD, aod::BkgChargedRhos>>::iterator const& collision,
                                         soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents> const& jets,
-                                        aod::JetTracks const&)
+                                        aod::JetTracks const&,
+                                        aod::JetMcCollisions const&)
   {
     bool fillHistograms = false;
     bool isWeighted = true;
     float eventWeight = collision.weight();
     if (!applyCollisionCuts(collision, fillHistograms, isWeighted, eventWeight)) {
       return;
     }
-
+    float pTHat = collision.has_mcCollision() && collision.mcCollision().ptHard() < 999.0f ? collision.mcCollision().ptHard() : 10. / (std::pow(eventWeight, 1.0 / pTHatExponent));
     float centrality = -1.0;
     checkCentFT0M ? centrality = collision.centFT0M() : centrality = collision.centFT0C();
 
@@ -895,7 +890,7 @@ struct JetSpectraCharged {
       if (!isAcceptedJet<aod::JetTracks>(jet)) {
         continue;
       }
-      fillJetAreaSubHistograms(jet, centrality, collision.rho(), eventWeight);
+      fillJetAreaSubHistograms(jet, centrality, collision.rho(), eventWeight, pTHat);
     }
   }
   PROCESS_SWITCH(JetSpectraCharged, processSpectraAreaSubMCDWeighted, "jet spectra with rho-area subtraction for MCD", false);
@@ -1127,22 +1122,21 @@ struct JetSpectraCharged {
     if (!applyMCCollisionCuts(mccollision, collisions, fillHistograms, isWeighted, eventWeight)) {
       return;
     }
-
+    float pTHat = mccollision.ptHard() < 999.0f ? mccollision.ptHard() : 10. / (std::pow(eventWeight, 1.0 / pTHatExponent));
     for (auto const& jet : jets) {
       if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
         continue;
       }
       if (!isAcceptedJet<aod::JetParticles>(jet, mcLevelIsParticleLevel)) {
         continue;
       }
-      double pTHat = 10. / (std::pow(eventWeight, 1.0 / pTHatExponent));
       int Nmax = 21;
       for (int N = 1; N < Nmax; N++) {
         if (jet.pt() < N * 0.25 * pTHat && jet.r() == round(selectedJetsRadius * 100.0f)) {
           registry.fill(HIST("h2_jet_ptcut_part"), jet.pt(), N * 0.25, eventWeight);
         }
       }
-      fillMCPHistograms(jet, eventWeight);
+      fillMCPHistograms(jet, eventWeight, pTHat);
     }
   }
   PROCESS_SWITCH(JetSpectraCharged, processSpectraMCPWeighted, "jet spectra for MC particle level weighted", false);
@@ -1160,6 +1154,7 @@ struct JetSpectraCharged {
     if (!applyMCCollisionCuts(mccollision, collisions, fillHistograms, isWeighted, eventWeight)) {
       return;
     }
+    float pTHat = mccollision.ptHard() < 999.0f ? mccollision.ptHard() : 10. / (std::pow(eventWeight, 1.0 / pTHatExponent));
     registry.fill(HIST("h_mccollisions_rho"), mccollision.rho(), eventWeight);
 
     for (auto const& jet : jets) {
@@ -1169,7 +1164,7 @@ struct JetSpectraCharged {
       if (!isAcceptedJet<aod::JetParticles>(jet, mcLevelIsParticleLevel)) {
         continue;
       }
-      fillMCPAreaSubHistograms(jet, mccollision.rho(), eventWeight);
+      fillMCPAreaSubHistograms(jet, mccollision.rho(), eventWeight, pTHat);
     }
   }
   PROCESS_SWITCH(JetSpectraCharged, processSpectraAreaSubMCPWeighted, "jet spectra with area-based subtraction for MC particle level", false);
@@ -1238,23 +1233,24 @@ struct JetSpectraCharged {
   }
   PROCESS_SWITCH(JetSpectraCharged, processJetsMatched, "matched mcp and mcd jets", false);
 
-  void processJetsMatchedWeighted(soa::Filtered<aod::JetCollisions>::iterator const& collision,
+  void processJetsMatchedWeighted(soa::Filtered<aod::JetCollisionsMCD>::iterator const& collision,
                                   ChargedMCDMatchedJets const& mcdjets,
                                   ChargedMCPMatchedJets const&,
-                                  aod::JetTracks const&, aod::JetParticles const&)
+                                  aod::JetTracks const&, aod::JetParticles const&,
+                                  aod::JetMcCollisions const&)
   {
     bool fillHistograms = false;
     bool isWeighted = true;
     float eventWeight = collision.weight();
     if (!applyCollisionCuts(collision, fillHistograms, isWeighted, eventWeight)) {
       return;
     }
-
+    float pTHat = collision.has_mcCollision() && collision.mcCollision().ptHard() < 999.0f ? collision.mcCollision().ptHard() : 10. / (std::pow(eventWeight, 1.0 / pTHatExponent));
     for (const auto& mcdjet : mcdjets) {
       if (!isAcceptedJet<aod::JetTracks>(mcdjet)) {
         continue;
       }
-      fillMatchedHistograms<ChargedMCDMatchedJets::iterator, ChargedMCPMatchedJets>(mcdjet, eventWeight);
+      fillMatchedHistograms<ChargedMCDMatchedJets::iterator, ChargedMCPMatchedJets>(mcdjet, eventWeight, pTHat);
     }
   }
   PROCESS_SWITCH(JetSpectraCharged, processJetsMatchedWeighted, "matched mcp and mcd jets with weighted events", false);
@@ -1292,14 +1288,15 @@ struct JetSpectraCharged {
     if (!applyCollisionCuts(collision, fillHistograms, isWeighted, eventWeight)) {
       return;
     }
+    float pTHat = collision.has_mcCollision() && collision.mcCollision().ptHard() < 999.0f ? collision.mcCollision().ptHard() : 10. / (std::pow(eventWeight, 1.0 / pTHatExponent));
 
     double mcrho = collision.has_mcCollision() ? collision.mcCollision_as<JetBkgRhoMcCollisions>().rho() : -1;
 
     for (const auto& mcdjet : mcdjets) {
       if (!isAcceptedJet<aod::JetTracks>(mcdjet)) {
         continue;
       }
-      fillGeoMatchedAreaSubHistograms<ChargedMCDMatchedJets::iterator, ChargedMCPMatchedJets>(mcdjet, collision.rho(), mcrho, eventWeight);
+      fillGeoMatchedAreaSubHistograms<ChargedMCDMatchedJets::iterator, ChargedMCPMatchedJets>(mcdjet, collision.rho(), mcrho, eventWeight, pTHat);
     }
   }
   PROCESS_SWITCH(JetSpectraCharged, processJetsMatchedAreaSubWeighted, "matched mcp and mcd jets after area-based pt subtraction with weighted events", false);