ATLAS Offline Software
RCJetSubstructureAug.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // RCJetSubstructureAug.cxx, (c) ATLAS Detector software
8 // Author: G. Albouy (galbouy@lpsc.in2p3.fr)
9 // This tool computes substructure variables for ReClustered jets
10 // from LCTopo clusters ghost associated to RC jets
11 // by constructing cluster jets
12 
13 #include "DerivationFrameworkLLP/RCJetSubstructureAug.h"
14 #include "xAODJet/JetContainer.h"
15 #include "xAODBase/IParticle.h"
16 #include "xAODCaloEvent/CaloCluster.h"
17 #include "StoreGate/WriteDecorHandle.h"
18 #include <vector>
19 
20 #include "fastjet/PseudoJet.hh"
21 #include "fastjet/ClusterSequence.hh"
22 #include "fastjet/JetDefinition.hh"
23 #include "fastjet/contrib/SoftDrop.hh"
24 #include "fastjet/tools/Filter.hh"
25 #include "fastjet/contrib/Nsubjettiness.hh"
26 #include "JetSubStructureUtils/Qw.h"
27 #include "JetSubStructureUtils/Nsubjettiness.h"
28 #include "JetSubStructureUtils/KtSplittingScale.h"
29 #include "JetSubStructureUtils/EnergyCorrelator.h"
30 
31 // Constructor
32 DerivationFramework::RCJetSubstructureAug::RCJetSubstructureAug(const std::string& t, const std::string& n, const IInterface* p) :
33 base_class(t, n, p) {}
34 
35 // Destructor
36 DerivationFramework::RCJetSubstructureAug::~RCJetSubstructureAug()
37 {}
38 
39 // Athena initialize
40 StatusCode DerivationFramework::RCJetSubstructureAug::initialize()
41 {
42  ATH_MSG_VERBOSE("initialize() ...");
43  if (m_jetKey.key().empty()) {
44  ATH_MSG_FATAL("No jet collection provided for augmentation.");
45  return StatusCode::FAILURE;
46  }
47  ATH_CHECK( m_jetKey.initialize() );
48 
49  auto init_dec = [&] ( WDHK& k ) {
50  k = k.key() + m_suffix;
51  return k.initialize();
52  };
53  ATH_CHECK( init_dec (m_dec_Qw) );
54  ATH_CHECK( init_dec (m_dec_Tau1) );
55  ATH_CHECK( init_dec (m_dec_Tau2) );
56  ATH_CHECK( init_dec (m_dec_Tau3) );
57  ATH_CHECK( init_dec (m_dec_Tau4) );
58  ATH_CHECK( init_dec (m_dec_Tau21) );
59  ATH_CHECK( init_dec (m_dec_Tau32) );
60  ATH_CHECK( init_dec (m_dec_Split12) );
61  ATH_CHECK( init_dec (m_dec_Split23) );
62  ATH_CHECK( init_dec (m_dec_Split34) );
63  ATH_CHECK( init_dec (m_dec_ECF1) );
64  ATH_CHECK( init_dec (m_dec_ECF2) );
65  ATH_CHECK( init_dec (m_dec_ECF3) );
66  ATH_CHECK( init_dec (m_dec_ECF4) );
67  ATH_CHECK( init_dec (m_dec_C2) );
68  ATH_CHECK( init_dec (m_dec_D2) );
69  ATH_CHECK( init_dec (m_dec_pT) );
70  ATH_CHECK( init_dec (m_dec_m) );
71  ATH_CHECK( init_dec (m_dec_NConstits) );
72  ATH_CHECK( init_dec (m_dec_eta) );
73  ATH_CHECK( init_dec (m_dec_phi) );
74  ATH_CHECK( init_dec (m_dec_timing) );
75 
76  // Set up the text-parsing machinery for selectiong the jet directly according to user cuts
77  if (!m_selectionString.empty()) {
78  ATH_CHECK( initializeParser( m_selectionString ));
79  }
80 
81  // Define the trimmer
82  m_trimmer.emplace(fastjet::JetDefinition(fastjet::kt_algorithm, m_rclus), fastjet::SelectorPtFractionMin(m_ptfrac));
83  m_softdropper.emplace(m_beta, m_zcut, m_R0);
84 
85  return StatusCode::SUCCESS;
86 }
87 
88 StatusCode DerivationFramework::RCJetSubstructureAug::addBranches(const EventContext& ctx) const
89 {
90 
91  SG::ReadHandle<xAOD::JetContainer> jets(m_jetKey,ctx);
92  if (!jets.isValid()) {
93  ATH_MSG_ERROR("No jet collection with name " << m_jetKey.key() << " found in StoreGate!");
94  return StatusCode::FAILURE;
95  }
96  unsigned int nJets(jets->size());
97  std::vector<const xAOD::Jet*> jetToCheck; jetToCheck.clear();
98 
99  // Execute the text parser if requested
100  if (!m_selectionString.empty()) {
101  std::vector<int> entries = m_parser->evaluateAsVector();
102  unsigned int nEntries = entries.size();
103  // check the sizes are compatible
104  if (nJets != nEntries ) {
105  ATH_MSG_ERROR("Sizes incompatible! Are you sure your selection string used jets??");
106  return StatusCode::FAILURE;
107  } else {
108  // identify which jets to keep for the thinning check
109  for (unsigned int i=0; i<nJets; ++i) if (entries[i]==1) jetToCheck.push_back((*jets)[i]);
110  }
111  } else {
112  for (unsigned int i=0; i<nJets; ++i) jetToCheck.push_back((*jets)[i]);
113  }
114 
115  using WDH = SG::WriteDecorHandle<xAOD::JetContainer, float>;
116  WDH dec_Qw (m_dec_Qw, ctx);
117  WDH dec_Tau1 (m_dec_Tau1, ctx);
118  WDH dec_Tau2 (m_dec_Tau2, ctx);
119  WDH dec_Tau3 (m_dec_Tau3, ctx);
120  WDH dec_Tau4 (m_dec_Tau4, ctx);
121  WDH dec_Tau21 (m_dec_Tau21, ctx);
122  WDH dec_Tau32 (m_dec_Tau32, ctx);
123  WDH dec_Split12 (m_dec_Split12, ctx);
124  WDH dec_Split23 (m_dec_Split23, ctx);
125  WDH dec_Split34 (m_dec_Split34, ctx);
126  WDH dec_ECF1 (m_dec_ECF1, ctx);
127  WDH dec_ECF2 (m_dec_ECF2, ctx);
128  WDH dec_ECF3 (m_dec_ECF3, ctx);
129  WDH dec_ECF4 (m_dec_ECF4, ctx);
130  WDH dec_C2 (m_dec_C2, ctx);
131  WDH dec_D2 (m_dec_D2, ctx);
132  WDH dec_pT (m_dec_pT, ctx);
133  WDH dec_m (m_dec_m, ctx);
134  WDH dec_NConstits (m_dec_NConstits, ctx);
135  WDH dec_eta (m_dec_eta, ctx);
136  WDH dec_phi (m_dec_phi, ctx);
137  WDH dec_timing (m_dec_timing, ctx);
138 
139  std::vector<const xAOD::IParticle*> mergedGhostConstits;
140  std::vector<const xAOD::IParticle*> ghosts;
141  std::vector<fastjet::PseudoJet> constituents;
142 
143  for (const xAOD::Jet* jet : jetToCheck) {
144 
145  // Get list of ghost constituents
146  mergedGhostConstits.clear();
147  for (std::string ghostName: m_ghostNames) {
148  ghosts.clear();
149  ATH_CHECK(jet->getAssociatedObjects(ghostName, ghosts));
150  mergedGhostConstits.insert(mergedGhostConstits.end(), ghosts.begin(), ghosts.end());
151  }
152 
153  // Construct list of constituent PseudoJets from constituents and compute timing information
154  constituents.clear();
155  double eTot = 0;
156  double time = 0;
157  for (auto constit : mergedGhostConstits){
158  if (constit==nullptr) {
159  ATH_MSG_INFO("Invalid constituent link, skipping ...");
160  continue;
161  }
162  // Filter constituents on negative energy
163  if (constit->e()<0) {ATH_MSG_INFO("################## Negative energy cluster"); continue;}
164  constituents.push_back( fastjet::PseudoJet(constit->p4()) );
165 
166  // Timing info of constituent only for calo constituents
167  if (constit->type() == xAOD::Type::CaloCluster) {
168  auto caloConstit = static_cast<const xAOD::CaloCluster*> (constit);
169  double eConstit = caloConstit->e()* caloConstit->e();
170  time += caloConstit->time()* eConstit;
171  eTot += eConstit;
172  }
173  }
174 
175  // Fill timing info
176  if (eTot==0) {
177  dec_timing(*jet) = 0;
178  } else {
179  dec_timing(*jet) = time/eTot;
180  }
181 
182  // Get number of constituents
183  int nConstits = constituents.size();
184 
185  // Run clustering on constituents if not empty
186  fastjet::PseudoJet groomed_jet;
187  if (nConstits!=0) {
188  auto jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm, 1.5);
189  fastjet::ClusterSequence cs(constituents, jet_def);
190  fastjet::PseudoJet recluster_jet = cs.inclusive_jets(0.0).front();
191 
192  // Apply grooming to reclustered jet
193  if (m_grooming=="Trimming"){
194  groomed_jet = m_trimmer->result(recluster_jet);
195  } else if (m_grooming=="SoftDrop"){
196  groomed_jet = m_softdropper->result(recluster_jet);
197  } else {
198  ATH_MSG_DEBUG(" No grooming requested or wrong one, will not apply grooming");
199  groomed_jet = recluster_jet;
200  }
201 
202  // update nConstit
203  nConstits = groomed_jet.constituents().size();
204  }
205 
206  // Fill substructure vars with default values when no constituents in jet
207  if (nConstits==0) {
208  dec_Qw(*jet) = -999;
209 
210  dec_Tau1(*jet) = -999;
211  dec_Tau2(*jet) = -999;
212  dec_Tau3(*jet) = -999;
213  dec_Tau4(*jet) = -999;
214 
215  dec_Tau21(*jet) = -999;
216  dec_Tau32(*jet) = -999;
217 
218  dec_Split12(*jet) = -999;
219  dec_Split23(*jet) = -999;
220  dec_Split34(*jet) = -999;
221 
222  dec_ECF1(*jet) = -999;
223  dec_ECF2(*jet) = -999;
224  dec_ECF3(*jet) = -999;
225  dec_ECF4(*jet) = -999;
226 
227  dec_C2(*jet) = -999;
228  dec_D2(*jet) = -999;
229 
230  dec_pT(*jet) = 0;
231  dec_m(*jet) = 0;
232  dec_NConstits(*jet) = 0;
233  dec_eta(*jet) = -999;
234  dec_phi(*jet) = -999;
235 
236  return StatusCode::SUCCESS;
237  }
238 
239  // Save reclustered jet info
240  dec_pT(*jet) = groomed_jet.pt();
241  dec_m(*jet) = groomed_jet.m();
242  dec_NConstits(*jet) = nConstits;
243  dec_eta(*jet) = groomed_jet.eta();
244  dec_phi(*jet) = groomed_jet.phi() - M_PI;
245 
246  // Qw
247  JetSubStructureUtils::Qw qw;
248  dec_Qw(*jet) = qw.result(groomed_jet);
249 
250  // Nsubjetiness
251  fastjet::contrib::WTA_KT_Axes wta_kt_axes;
252  fastjet::contrib::NormalizedCutoffMeasure normalized_measure(1.0, 1.0, 1000000);
253  JetSubStructureUtils::Nsubjettiness Tau1_wta(1, wta_kt_axes, normalized_measure);
254  JetSubStructureUtils::Nsubjettiness Tau2_wta(2, wta_kt_axes, normalized_measure);
255  JetSubStructureUtils::Nsubjettiness Tau3_wta(3, wta_kt_axes, normalized_measure);
256  JetSubStructureUtils::Nsubjettiness Tau4_wta(4, wta_kt_axes, normalized_measure);
257  float tau1_wta = Tau1_wta.result(groomed_jet);
258  float tau2_wta = Tau2_wta.result(groomed_jet);
259  float tau3_wta = Tau3_wta.result(groomed_jet);
260  float tau4_wta = Tau4_wta.result(groomed_jet);
261  float tau21_wta = -999;
262  float tau32_wta = -999;
263 
264  if( tau1_wta > 1e-8 ) {
265  tau21_wta = tau2_wta / tau1_wta;
266  }
267  if( tau2_wta > 1e-8 ) {
268  tau32_wta = tau3_wta / tau2_wta;
269  }
270 
271  dec_Tau1(*jet) = tau1_wta;
272  dec_Tau2(*jet) = tau2_wta;
273  dec_Tau3(*jet) = tau3_wta;
274  dec_Tau4(*jet) = tau4_wta;
275  dec_Tau21(*jet) = tau21_wta;
276  dec_Tau32(*jet) = tau32_wta;
277 
278  // KtSplittingScale
279  JetSubStructureUtils::KtSplittingScale Split12(1);
280  JetSubStructureUtils::KtSplittingScale Split23(2);
281  JetSubStructureUtils::KtSplittingScale Split34(3);
282  float split12 = Split12.result(groomed_jet);
283  float split23 = Split23.result(groomed_jet);
284  float split34 = Split34.result(groomed_jet);
285 
286  dec_Split12(*jet) = split12;
287  dec_Split23(*jet) = split23;
288  dec_Split34(*jet) = split34;
289 
290  // EnergyCorrelator
291  JetSubStructureUtils::EnergyCorrelator ECF1(1, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
292  JetSubStructureUtils::EnergyCorrelator ECF2(2, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
293  JetSubStructureUtils::EnergyCorrelator ECF3(3, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
294  JetSubStructureUtils::EnergyCorrelator ECF4(4, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
295  float ecf1 = ECF1.result(groomed_jet);
296  float ecf2 = ECF2.result(groomed_jet);
297  float ecf3 = ECF3.result(groomed_jet);
298  float ecf4 = ECF4.result(groomed_jet);
299  float c2 = -999;
300  float d2 = -999;
301 
302  if( ecf2 > 1e-8 ) {
303  c2 = ecf3 * ecf1 / pow( ecf2, 2.0 );
304  d2 = ecf3 * pow( ecf1, 3.0 ) / pow( ecf2, 3.0 );
305  }
306 
307  dec_ECF1(*jet) = ecf1;
308  dec_ECF2(*jet) = ecf2;
309  dec_ECF3(*jet) = ecf3;
310  dec_ECF4(*jet) = ecf4;
311  dec_C2(*jet) = c2;
312  dec_D2(*jet) = d2;
313 
314  }
315 
316  return StatusCode::SUCCESS;
317 }
318 
319 
SG::WriteDecorHandleKey< xAOD::JetContainer >
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
KtSplittingScale.h
JetSubStructureUtils::KtSplittingScale
Definition: KtSplittingScale.h:11
IParticle.h
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
xAOD::JetAlgorithmType::kt_algorithm
@ kt_algorithm
Definition: JetContainerInfo.h:31
JetSubStructureUtils::Nsubjettiness::result
virtual double result(const fastjet::PseudoJet &jet) const
Definition: Nsubjettiness.h:21
JetSubStructureUtils::EnergyCorrelator::result
virtual double result(const fastjet::PseudoJet &jet) const
Definition: EnergyCorrelator.h:20
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:67
M_PI
#define M_PI
Definition: ActiveFraction.h:11
defineDB.jets
jets
Definition: JetTagCalibration/share/defineDB.py:24
JetSubStructureUtils::EnergyCorrelator
Definition: EnergyCorrelator.h:14
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
IDTPM::eTot
float eTot(const U &p)
Accessor utility function for getting the value of Energy.
Definition: TrackParametersHelper.h:118
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
xAOD::CaloCluster
CaloCluster_v1 CaloCluster
Define the latest version of the calorimeter cluster class.
Definition: Event/xAOD/xAODCaloEvent/xAODCaloEvent/CaloCluster.h:19
xAOD::CaloCluster_v1
Description of a calorimeter cluster.
Definition: CaloCluster_v1.h:62
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
jet
Definition: JetCalibTools_PlotJESFactors.cxx:23
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
DerivationFramework::RCJetSubstructureAug::initialize
virtual StatusCode initialize() override
Definition: RCJetSubstructureAug.cxx:40
jet::ClusterSequence
fastjet::ClusterSequence ClusterSequence
Definition: ClusterSequence.h:21
lumiFormat.i
int i
Definition: lumiFormat.py:85
CaloCluster.h
beamspotman.n
n
Definition: beamspotman.py:727
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
SG::WriteDecorHandle
Handle class for adding a decoration to an object.
Definition: StoreGate/StoreGate/WriteDecorHandle.h:100
WriteDecorHandle.h
Handle class for adding a decoration to an object.
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Nsubjettiness.h
DerivationFramework::RCJetSubstructureAug::addBranches
virtual StatusCode addBranches(const EventContext &ctx) const override
Definition: RCJetSubstructureAug.cxx:88
python.StandardJetMods.qw
qw
Definition: StandardJetMods.py:313
xAOD::JetAlgorithmType::antikt_algorithm
@ antikt_algorithm
Definition: JetContainerInfo.h:33
xAOD::Jet_v1
Class describing a jet.
Definition: Jet_v1.h:57
python.DataFormatRates.c2
c2
Definition: DataFormatRates.py:123
JetSubStructureUtils::Nsubjettiness
Definition: Nsubjettiness.h:14
DerivationFramework::RCJetSubstructureAug::RCJetSubstructureAug
RCJetSubstructureAug(const std::string &t, const std::string &n, const IInterface *p)
Definition: RCJetSubstructureAug.cxx:32
DerivationFramework::RCJetSubstructureAug::~RCJetSubstructureAug
virtual ~RCJetSubstructureAug()
Definition: RCJetSubstructureAug.cxx:36
JetContainer.h
CaloSwCorrections.time
def time(flags, cells_name, *args, **kw)
Definition: CaloSwCorrections.py:242
Qw.h
dq_defect_virtual_defect_validation.d2
d2
Definition: dq_defect_virtual_defect_validation.py:81
entries
double entries
Definition: listroot.cxx:49
JetSubStructureUtils::Qw
Definition: Qw.h:34
dqBeamSpot.nEntries
int nEntries
Definition: dqBeamSpot.py:72
EnergyCorrelator.h
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
JetSubStructureUtils::KtSplittingScale::result
virtual double result(const fastjet::PseudoJet &jet) const
Definition: KtSplittingScale.cxx:12
xAOD::CaloCluster_v1::e
virtual double e() const
The total energy of the particle.
Definition: CaloCluster_v1.cxx:265
fitman.k
k
Definition: fitman.py:528
RCJetSubstructureAug.h
Generated on Sat Nov 1 2025 21:17:07 for ATLAS Offline Software by   doxygen 1.8.18