Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
nzardosh merged 13 commits into from
Mar 31, 2026
Merged

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

Changes from 11 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
4edc6cb
fixes to jetCorrelationD0.cxx task
MattOckleton Mar 12, 2026
69133a6
Updating columns for jetCorrelation D0 task
MattOckleton Mar 12, 2026
67c038f
Adding processMcMatched and mirroring dPhi dists
MattOckleton Mar 25, 2026
f8260f4
adding dPhi to processMcMatched
MattOckleton Mar 26, 2026
d5f18b2
formatting
MattOckleton Mar 26, 2026
6320a60
fixing headers
MattOckleton Mar 26, 2026
707436a
Fixes to merge conflicts
MattOckleton Mar 26, 2026
086e3ca
formatting
MattOckleton Mar 26, 2026
10add02
Fixes to jetCorrelationD0 task
MattOckleton Mar 27, 2026
48eb8e6
Formatting
MattOckleton Mar 27, 2026
3f6b4b7
Formatting
MattOckleton Mar 27, 2026
d07dbbe
dPhi fixes
MattOckleton Mar 31, 2026
5839fd9
Final dphi fixes
MattOckleton Mar 31, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
213 changes: 137 additions & 76 deletions PWGJE/Tasks/jetCorrelationD0.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand All @@ -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
{
Expand All @@ -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,
Expand Down Expand Up @@ -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,
Expand All @@ -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
Expand All @@ -179,13 +195,13 @@ 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_pT_D0_Jet_dPhi"), jet.pt(), dPhi);
}

template <typename T>
Expand All @@ -199,26 +215,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));
Expand Down Expand Up @@ -284,23 +280,26 @@ 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());
if (dPhi > o2::constants::math::PI) {
dPhi = 2 * o2::constants::math::PI - dPhi;
}
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)
{
Expand All @@ -327,60 +326,122 @@ 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());
if (dPhi > o2::constants::math::PI) {
dPhi = 2 * o2::constants::math::PI - dPhi;
}
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());
if (dPhi > o2::constants::math::PI) {
dPhi = 2 * o2::constants::math::PI - dPhi;
}
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());
if (dPhiD > o2::constants::math::PI) {
dPhiD = 2 * o2::constants::math::PI - dPhiD;
}
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());
if (dPhiP > o2::constants::math::PI) {
Copy link
Copy Markdown
Collaborator

@nzardosh nzardosh Mar 30, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i think constrain angle should already take care of the cyclical nature?

Copy link
Copy Markdown
Contributor Author

@MattOckleton MattOckleton Mar 31, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed

dPhiP = 2 * o2::constants::math::PI - dPhiP;
}
// 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)
Expand Down
Loading