[PWGJE] Adding Matched MC function + dPhi fixes (#15551)

Author: MattOckleton

PWGJE/Tasks/jetCorrelationD0.cxx

@@ -14,24 +14,18 @@
 /// \author Matthew Ockleton matthew.ockleton@cern.ch, University of Liverpool
 
 #include "PWGJE/Core/JetDerivedDataUtilities.h"
+#include "PWGJE/Core/JetHFUtilities.h"
 #include "PWGJE/DataModel/Jet.h"
 #include "PWGJE/DataModel/JetReducedData.h"
 
 #include "Common/Core/RecoDecay.h"
 
-#include <CommonConstants/MathConstants.h>
-#include <Framework/ASoA.h>
-#include <Framework/AnalysisDataModel.h>
-#include <Framework/AnalysisHelpers.h>
-#include <Framework/AnalysisTask.h>
-#include <Framework/Configurable.h>
-#include <Framework/HistogramRegistry.h>
-#include <Framework/HistogramSpec.h>
-#include <Framework/InitContext.h>
-#include <Framework/OutputObjHeader.h>
-#include <Framework/runDataProcessing.h>
-
-#include <cstdlib>
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/HistogramRegistry.h"
+#include "Framework/Logger.h"
+#include "Framework/runDataProcessing.h"
+
 #include <string>
 #include <type_traits>
 #include <vector>
@@ -58,12 +52,15 @@ DECLARE_SOA_TABLE(McCollisionTables, "AOD", "MCCOLLINFOTABLE",
                   o2::soa::Index<>,
                   d0collisionInfo::PosZ);
 
+DECLARE_SOA_TABLE(MatchCollTables, "AOD", "MATCHCOLLTABLE",
+                  o2::soa::Index<>,
+                  d0collisionInfo::PosZ);
+
 namespace collisionInfo
 {
-// DECLARE_SOA_INDEX_COLUMN(CollisionTable, collisionTable);
 DECLARE_SOA_INDEX_COLUMN_CUSTOM(CollisionTable, collisionTable, "COLLINFOTABLES");
-// DECLARE_SOA_INDEX_COLUMN(McCollisionTable, mcCollisionTable);
 DECLARE_SOA_INDEX_COLUMN_CUSTOM(McCollisionTable, mcCollisionTable, "MCCOLLINFOTABLES");
+DECLARE_SOA_INDEX_COLUMN_CUSTOM(MatchCollTable, matchCollTable, "MATCHCOLLTABLES");
 } // namespace collisionInfo
 namespace d0Info
 {
@@ -84,7 +81,7 @@ DECLARE_SOA_COLUMN(D0PhiD, d0PhiD, float);
 DECLARE_SOA_COLUMN(D0Reflection, d0Reflection, int);
 } // namespace d0Info
 
-DECLARE_SOA_TABLE(D0DataTables, "AOD", "D0DATATABLE",
+DECLARE_SOA_TABLE(D0DataTables, "AOD", "D0TABLE",
                   o2::soa::Index<>,
                   collisionInfo::CollisionTableId,
                   d0Info::D0PromptBDT,
@@ -114,11 +111,15 @@ DECLARE_SOA_INDEX_COLUMN(D0McPTable, d0McPTable);
 DECLARE_SOA_COLUMN(JetPt, jetPt, float);
 DECLARE_SOA_COLUMN(JetEta, jetEta, float);
 DECLARE_SOA_COLUMN(JetPhi, jetPhi, float);
+DECLARE_SOA_COLUMN(PJetPt, pJetPt, float);
+DECLARE_SOA_COLUMN(PJetEta, pJetEta, float);
+DECLARE_SOA_COLUMN(PJetPhi, pJetPhi, float);
 // D0-jet
 DECLARE_SOA_COLUMN(D0JetDeltaPhi, d0JetDeltaPhi, float);
+DECLARE_SOA_COLUMN(D0JetDeltaPhiP, d0JetDeltaPhiP, float);
 } // namespace jetInfo
 
-DECLARE_SOA_TABLE_STAGED(JetDataTables, "JETDATATABLE",
+DECLARE_SOA_TABLE_STAGED(JetDataTables, "JETTABLE",
                          o2::soa::Index<>,
                          collisionInfo::CollisionTableId,
                          jetInfo::D0DataTableId,
@@ -127,37 +128,52 @@ DECLARE_SOA_TABLE_STAGED(JetDataTables, "JETDATATABLE",
                          jetInfo::JetPhi,
                          jetInfo::D0JetDeltaPhi);
 
-DECLARE_SOA_TABLE_STAGED(JetMCPTables, "JETMCPTABLE",
+DECLARE_SOA_TABLE_STAGED(JetMcPTables, "JETMCPTABLE",
                          o2::soa::Index<>,
                          collisionInfo::McCollisionTableId,
                          jetInfo::D0McPTableId,
                          jetInfo::JetPt,
                          jetInfo::JetEta,
                          jetInfo::JetPhi,
-                         jetInfo::D0JetDeltaPhi);
+                         jetInfo::D0JetDeltaPhiP);
+
+DECLARE_SOA_TABLE_STAGED(JetMatchedTables, "JETMATCHEDTABLE",
+                         o2::soa::Index<>,
+                         collisionInfo::MatchCollTableId,
+                         jetInfo::JetPt,
+                         jetInfo::JetEta,
+                         jetInfo::JetPhi,
+                         jetInfo::PJetPt,
+                         jetInfo::PJetEta,
+                         jetInfo::PJetPhi,
+                         jetInfo::D0JetDeltaPhi,
+                         jetInfo::D0JetDeltaPhiP);
 
 } // namespace o2::aod
 
 struct JetCorrelationD0 {
   // Define new table
   Produces<aod::CollisionTables> tableCollision;
+  Produces<aod::MatchCollTables> tableMatchedCollision;
   Produces<aod::McCollisionTables> tableMcCollision;
   Produces<aod::D0DataTables> tableD0;
-  Produces<aod::D0McPTables> tableD0MCParticle;
+  Produces<aod::D0McPTables> tableD0McParticle;
   Produces<aod::JetDataTables> tableJet;
-  Produces<aod::JetMCPTables> tableJetMCParticle;
+  Produces<aod::JetMcPTables> tableJetMcParticle;
+  Produces<aod::JetMatchedTables> tableJetMatched;
 
   // Configurables
   Configurable<std::string> eventSelections{"eventSelections", "sel8", "choose event selection"};
   Configurable<bool> skipMBGapEvents{"skipMBGapEvents", false, "decide to run over MB gap events or not"};
   Configurable<bool> applyRCTSelections{"applyRCTSelections", true, "decide to apply RCT selections"};
-  // Configurable<std::string> triggerMasks{"triggerMasks", "", "possible JE Trigger masks: fJetChLowPt,fJetChHighPt,fTrackLowPt,fTrackHighPt,fJetD0ChLowPt,fJetD0ChHighPt,fJetLcChLowPt,fJetLcChHighPt,fEMCALReadout,fJetFullHighPt,fJetFullLowPt,fJetNeutralHighPt,fJetNeutralLowPt,fGammaVeryHighPtEMCAL,fGammaVeryHighPtDCAL,fGammaHighPtEMCAL,fGammaHighPtDCAL,fGammaLowPtEMCAL,fGammaLowPtDCAL,fGammaVeryLowPtEMCAL,fGammaVeryLowPtDCAL"};
   Configurable<float> jetPtCutMin{"jetPtCutMin", 5.0, "minimum value of jet pt"};
   Configurable<float> d0PtCutMin{"d0PtCutMin", 1.0, "minimum value of d0 pt"};
+  Configurable<float> jetMcPtCutMin{"jetMcPtCutMin", 3.0, "minimum value of jet pt particle level"};
+  Configurable<float> d0McPtCutMin{"d0McPtCutMin", 0.5, "minimum value of d0 pt particle level"};
   Configurable<float> vertexZCut{"vertexZCut", 10.0, "Accepted z-vertex range"};
   Configurable<float> pTHatExponent{"pTHatExponent", 6.0, "exponent of the event weight for the calculation of pTHat"};
-  Configurable<float> pTHatMaxMCD{"pTHatMaxMCD", 999.0, "maximum fraction of hard scattering for jet acceptance in detector MC"};
-  Configurable<float> pTHatMaxMCP{"pTHatMaxMCP", 999.0, "maximum fraction of hard scattering for jet acceptance in particle MC"};
+  Configurable<float> pTHatMaxMcD{"pTHatMaxMcD", 999.0, "maximum fraction of hard scattering for jet acceptance in detector MC"};
+  Configurable<float> pTHatMaxMcP{"pTHatMaxMcP", 999.0, "maximum fraction of hard scattering for jet acceptance in particle MC"};
   Configurable<float> pTHatAbsoluteMin{"pTHatAbsoluteMin", -99.0, "minimum value of pTHat"};
 
   // Filters
@@ -179,13 +195,14 @@ struct JetCorrelationD0 {
     registry.fill(HIST("hD0Phi"), d0.phi());
   }
   template <typename T, typename U>
-  void fillJetHistograms(T const& jet, U const& dphi)
+  void fillJetHistograms(T const& jet, U const& dPhi)
   {
     registry.fill(HIST("hJetPt"), jet.pt());
     registry.fill(HIST("hJetEta"), jet.eta());
     registry.fill(HIST("hJetPhi"), jet.phi());
     registry.fill(HIST("hJet3D"), jet.pt(), jet.eta(), jet.phi());
-    registry.fill(HIST("h_Jet_pT_D0_Jet_dPhi"), jet.pt(), dphi);
+    registry.fill(HIST("h_Jet_D0_Jet_dPhi"), dPhi);
+    registry.fill(HIST("h_Jet_pT_D0_Jet_dPhi"), jet.pt(), dPhi);
   }
 
   template <typename T>
@@ -199,26 +216,6 @@ struct JetCorrelationD0 {
     registry.fill(HIST("hZvtxSelected"), collision.posZ());
     return true;
   }
-
-  template <typename T, typename U>
-  // Jetbase is an MCD jet. We then loop through jettagv(MCP jets) to test if they match
-  // void fillMatchedHistograms(T const& jetBase, float weight = 1.0) // float leadingTrackPtBase,
-  void fillMatchedHistograms(T const& jetsBase, U const&, float weight = 1.0, float rho = 0.0)
-  {
-    for (const auto& jetBase : jetsBase) {
-      if (jetBase.has_matchedJetGeo()) { // geometric matching
-        for (auto const& jetTag : jetBase.template matchedJetGeo_as<std::decay_t<U>>()) {
-          registry.fill(HIST("hPtMatched"), jetBase.pt() - (rho * jetBase.area()), jetTag.pt(), weight);
-          registry.fill(HIST("hPtMatched1d"), jetTag.pt(), weight);
-          registry.fill(HIST("hPhiMatched"), jetBase.phi(), jetTag.phi(), weight);
-          registry.fill(HIST("hEtaMatched"), jetBase.eta(), jetTag.eta(), weight);
-          registry.fill(HIST("hPtResolution"), jetTag.pt(), (jetTag.pt() - (jetBase.pt() - (rho * jetBase.area()))) / jetTag.pt(), weight);
-          registry.fill(HIST("hPhiResolution"), jetTag.pt(), jetTag.phi() - jetBase.phi(), weight);
-          registry.fill(HIST("hEtaResolution"), jetTag.pt(), jetTag.eta() - jetBase.eta(), weight);
-        }
-      }
-    }
-  }
   void init(InitContext const&)
   {
     eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast<std::string>(eventSelections));
@@ -246,6 +243,7 @@ struct JetCorrelationD0 {
     registry.add("hJetEta", "jet #eta;#eta_{jet};entries", HistType::kTH1F, {axisEta});
     registry.add("hJetPhi", "jet #phi;#phi_{jet};entries", HistType::kTH1F, {axisPhi});
     registry.add("hJet3D", "3D jet distribution;p_{T};#eta;#phi", {HistType::kTH3F, {{500, -100, 400}, {100, -1.0, 1.0}, {100, 0.0, o2::constants::math::TwoPI}}});
+    registry.add("h_Jet_D0_Jet_dPhi", "#Delta #phi _{D^{0}, jet}", kTH1F, {{100, 0, o2::constants::math::TwoPI}});
     registry.add("h_Jet_pT_D0_Jet_dPhi", "p_{T, jet} vs #Delta #phi _{D^{0}, jet}", kTH2F, {{100, 0, 100}, {100, 0, o2::constants::math::TwoPI}});
 
     // Matching histograms
@@ -284,23 +282,23 @@ struct JetCorrelationD0 {
         if (jet.pt() < jetPtCutMin) {
           continue;
         }
-        float dphi = RecoDecay::constrainAngle(jet.phi() - d0Candidate.phi());
-        if (std::abs(dphi - o2::constants::math::PI) > (o2::constants::math::PI / 2)) { // this is quite loose  instead of pi/2 could do 0.6
+        float dPhi = RecoDecay::constrainAngle(jet.phi() - d0Candidate.phi(), -o2::constants::math::PI);
+        if (std::abs(dPhi - o2::constants::math::PI) > (o2::constants::math::PI / 2)) {
           continue;
         }
-        fillJetHistograms(jet, dphi);
+        fillJetHistograms(jet, dPhi);
         tableJet(tableCollision.lastIndex(),
                  tableD0.lastIndex(),
                  jet.pt(),
                  jet.eta(),
                  jet.phi(),
-                 dphi);
+                 dPhi);
       }
     }
   }
   PROCESS_SWITCH(JetCorrelationD0, processData, "charged particle level jet analysis", true);
 
-  void processMCDetector(soa::Filtered<aod::JetCollisions>::iterator const& collision,
+  void processMcDetector(soa::Filtered<aod::JetCollisions>::iterator const& collision,
                          aod::CandidatesD0MCD const& d0Candidates,
                          soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents> const& jets)
   {
@@ -327,60 +325,110 @@ struct JetCorrelationD0 {
         if (jet.pt() < jetPtCutMin) {
           continue;
         }
-        float dphi = RecoDecay::constrainAngle(jet.phi() - d0Candidate.phi());
-        if (std::abs(dphi - o2::constants::math::PI) > (o2::constants::math::PI / 2)) { // this is quite loose  instead of pi/2 could do 0.6
+        float dPhi = RecoDecay::constrainAngle(jet.phi() - d0Candidate.phi(), -o2::constants::math::PI);
+        if (std::abs(dPhi - o2::constants::math::PI) > (o2::constants::math::PI / 2)) {
           continue;
         }
-        fillJetHistograms(jet, dphi);
+        fillJetHistograms(jet, dPhi);
         tableJet(tableCollision.lastIndex(),
                  tableD0.lastIndex(),
                  jet.pt(),
                  jet.eta(),
                  jet.phi(),
-                 dphi);
+                 dPhi);
       }
     }
   }
-  PROCESS_SWITCH(JetCorrelationD0, processMCDetector, "charged particle level jet analysis", false);
+  PROCESS_SWITCH(JetCorrelationD0, processMcDetector, "charged detector level jet analysis", false);
 
-  void processMCParticle(aod::JetMcCollision const& collision,
-                         aod::CandidatesD0MCP const& d0MCPCandidates,
-                         soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents>> const& jets)
+  void processMcParticle(aod::JetMcCollision const& collision,
+                         aod::CandidatesD0MCP const& d0McPCandidates,
+                         soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents> const& jets)
   {
-    if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents, applyRCTSelections)) { // build without this
+    if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents, applyRCTSelections)) {
       return;
     }
     tableMcCollision(collision.posZ());
-    for (const auto& d0MCPCandidate : d0MCPCandidates) {
-      if (d0MCPCandidate.pt() < d0PtCutMin) {
+    for (const auto& d0McPCandidate : d0McPCandidates) {
+      if (d0McPCandidate.pt() < d0McPtCutMin) {
         continue;
       }
-      tableD0MCParticle(tableCollision.lastIndex(),
-                        d0MCPCandidate.originMcGen(),
-                        d0MCPCandidate.pt(),
-                        d0MCPCandidate.eta(),
-                        d0MCPCandidate.phi(),
-                        d0MCPCandidate.y());
+      tableD0McParticle(tableMcCollision.lastIndex(),
+                        d0McPCandidate.originMcGen(),
+                        d0McPCandidate.pt(),
+                        d0McPCandidate.eta(),
+                        d0McPCandidate.phi(),
+                        d0McPCandidate.y());
 
       for (const auto& jet : jets) {
-        if (jet.pt() < jetPtCutMin) {
+        if (jet.pt() < jetMcPtCutMin) {
           continue;
         }
-        float dphi = RecoDecay::constrainAngle(jet.phi() - d0MCPCandidate.phi());
-        if (std::abs(dphi - o2::constants::math::PI) > (o2::constants::math::PI / 2)) {
+        float dPhi = RecoDecay::constrainAngle(jet.phi() - d0McPCandidate.phi(), -o2::constants::math::PI);
+        if (std::abs(dPhi - o2::constants::math::PI) > (o2::constants::math::PI / 2)) {
           continue;
         }
-        fillJetHistograms(jet, dphi);
-        tableJetMCParticle(tableCollision.lastIndex(),
-                           tableD0MCParticle.lastIndex(),
+        fillJetHistograms(jet, dPhi);
+        tableJetMcParticle(tableMcCollision.lastIndex(),
+                           tableD0McParticle.lastIndex(),
                            jet.pt(),
                            jet.eta(),
                            jet.phi(),
-                           dphi);
+                           dPhi);
+      }
+    }
+  }
+  PROCESS_SWITCH(JetCorrelationD0, processMcParticle, "charged MC Particle jets", false);
+
+  void processMcMatched(soa::Filtered<aod::JetCollisions>::iterator const& collision,
+                        aod::CandidatesD0MCD const& d0Candidates,
+                        aod::JetTracksMCD const& tracks,
+                        aod::JetParticles const& particles,
+                        soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets> const& McDJets,
+                        aod::ChargedMCParticleLevelJets const&)
+  {
+    if (!applyCollisionSelections(collision)) {
+      return;
+    }
+    tableMatchedCollision(collision.posZ());
+    for (const auto& d0Candidate : d0Candidates) {
+      if (d0Candidate.pt() < d0PtCutMin) { // once settled on a mlcut, then add the lower bound of the systematics as a cut here
+        continue;
+      }
+      bool isMatched = false;
+      const auto& d0Particle = jethfutilities::matchedHFParticle(d0Candidate, tracks, particles, isMatched);
+      if (!isMatched) {
+        continue;
+      }
+      for (const auto& McDJet : McDJets) {
+        if (McDJet.pt() < jetPtCutMin) {
+          continue;
+        }
+        float dPhiD = RecoDecay::constrainAngle(McDJet.phi() - d0Candidate.phi(), -o2::constants::math::PI);
+        if (std::abs(dPhiD - o2::constants::math::PI) > (o2::constants::math::PI / 2)) {
+          continue;
+        }
+        if (McDJet.has_matchedJetGeo()) { // geometric matching
+          for (auto const& McPJet : McDJet.template matchedJetGeo_as<aod::ChargedMCParticleLevelJets>()) {
+            float dPhiP = RecoDecay::constrainAngle(McPJet.phi() - d0Particle.phi(), -o2::constants::math::PI);
+            // if (std::abs(dPhiP - o2::constants::math::PI) > (o2::constants::math::PI / 2)) {
+            //   continue;
+            // }
+            tableJetMatched(tableMatchedCollision.lastIndex(),
+                            McDJet.pt(),
+                            McDJet.eta(),
+                            McDJet.phi(),
+                            McPJet.pt(),
+                            McPJet.eta(),
+                            McPJet.phi(),
+                            dPhiD,
+                            dPhiP);
+          }
+        }
       }
     }
   }
-  PROCESS_SWITCH(JetCorrelationD0, processMCParticle, "process MC Particle jets", false);
+  PROCESS_SWITCH(JetCorrelationD0, processMcMatched, "process matching of particle level jets to detector level jets", false);
 };
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)