Loading [MathJax]/extensions/tex2jax.js
ATLAS Offline Software
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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
89 {
90  const EventContext& ctx = Gaudi::Hive::currentContext();
91 
92  SG::ReadHandle<xAOD::JetContainer> jets(m_jetKey,ctx);
93  if (!jets.isValid()) {
94  ATH_MSG_ERROR("No jet collection with name " << m_jetKey.key() << " found in StoreGate!");
95  return StatusCode::FAILURE;
96  }
97  unsigned int nJets(jets->size());
98  std::vector<const xAOD::Jet*> jetToCheck; jetToCheck.clear();
99 
100  // Execute the text parser if requested
101  if (!m_selectionString.empty()) {
102  std::vector<int> entries = m_parser->evaluateAsVector();
103  unsigned int nEntries = entries.size();
104  // check the sizes are compatible
105  if (nJets != nEntries ) {
106  ATH_MSG_ERROR("Sizes incompatible! Are you sure your selection string used jets??");
107  return StatusCode::FAILURE;
108  } else {
109  // identify which jets to keep for the thinning check
110  for (unsigned int i=0; i<nJets; ++i) if (entries[i]==1) jetToCheck.push_back((*jets)[i]);
111  }
112  } else {
113  for (unsigned int i=0; i<nJets; ++i) jetToCheck.push_back((*jets)[i]);
114  }
115 
116  using WDH = SG::WriteDecorHandle<xAOD::JetContainer, float>;
117  WDH dec_Qw (m_dec_Qw, ctx);
118  WDH dec_Tau1 (m_dec_Tau1, ctx);
119  WDH dec_Tau2 (m_dec_Tau2, ctx);
120  WDH dec_Tau3 (m_dec_Tau3, ctx);
121  WDH dec_Tau4 (m_dec_Tau4, ctx);
122  WDH dec_Tau21 (m_dec_Tau21, ctx);
123  WDH dec_Tau32 (m_dec_Tau32, ctx);
124  WDH dec_Split12 (m_dec_Split12, ctx);
125  WDH dec_Split23 (m_dec_Split23, ctx);
126  WDH dec_Split34 (m_dec_Split34, ctx);
127  WDH dec_ECF1 (m_dec_ECF1, ctx);
128  WDH dec_ECF2 (m_dec_ECF2, ctx);
129  WDH dec_ECF3 (m_dec_ECF3, ctx);
130  WDH dec_ECF4 (m_dec_ECF4, ctx);
131  WDH dec_C2 (m_dec_C2, ctx);
132  WDH dec_D2 (m_dec_D2, ctx);
133  WDH dec_pT (m_dec_pT, ctx);
134  WDH dec_m (m_dec_m, ctx);
135  WDH dec_NConstits (m_dec_NConstits, ctx);
136  WDH dec_eta (m_dec_eta, ctx);
137  WDH dec_phi (m_dec_phi, ctx);
138  WDH dec_timing (m_dec_timing, ctx);
139 
140  std::vector<const xAOD::IParticle*> mergedGhostConstits;
141  std::vector<const xAOD::IParticle*> ghosts;
142  std::vector<fastjet::PseudoJet> constituents;
143 
144  for (const xAOD::Jet* jet : jetToCheck) {
145 
146  // Get list of ghost constituents
147  mergedGhostConstits.clear();
148  for (std::string ghostName: m_ghostNames) {
149  ghosts.clear();
150  ATH_CHECK(jet->getAssociatedObjects(ghostName, ghosts));
151  mergedGhostConstits.insert(mergedGhostConstits.end(), ghosts.begin(), ghosts.end());
152  }
153 
154  // Construct list of constituent PseudoJets from constituents and compute timing information
155  constituents.clear();
156  double eTot = 0;
157  double time = 0;
158  for (auto constit : mergedGhostConstits){
159  if (constit==nullptr) {
160  ATH_MSG_INFO("Invalid constituent link, skipping ...");
161  continue;
162  }
163  // Filter constituents on negative energy
164  if (constit->e()<0) {ATH_MSG_INFO("################## Negative energy cluster"); continue;}
165  constituents.push_back( fastjet::PseudoJet(constit->p4()) );
166 
167  // Timing info of constituent only for calo constituents
168  if (constit->type() == xAOD::Type::CaloCluster) {
169  auto caloConstit = dynamic_cast<const xAOD::CaloCluster*> (constit);
170  double eConstit = caloConstit->e()* caloConstit->e();
171  time += caloConstit->time()* eConstit;
172  eTot += eConstit;
173  }
174  }
175 
176  // Fill timing info
177  if (eTot==0) {
178  dec_timing(*jet) = 0;
179  } else {
180  dec_timing(*jet) = time/eTot;
181  }
182 
183  // Get number of constituents
184  int nConstits = constituents.size();
185 
186  // Run clustering on constituents if not empty
187  fastjet::PseudoJet groomed_jet;
188  if (nConstits!=0) {
189  auto jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm, 1.5);
190  fastjet::ClusterSequence cs(constituents, jet_def);
191  fastjet::PseudoJet recluster_jet = cs.inclusive_jets(0.0).front();
192 
193  // Apply grooming to reclustered jet
194  if (m_grooming=="Trimming"){
195  groomed_jet = m_trimmer->result(recluster_jet);
196  } else if (m_grooming=="SoftDrop"){
197  groomed_jet = m_softdropper->result(recluster_jet);
198  } else {
199  ATH_MSG_DEBUG(" No grooming requested or wrong one, will not apply grooming");
200  groomed_jet = recluster_jet;
201  }
202 
203  // update nConstit
204  nConstits = groomed_jet.constituents().size();
205  }
206 
207  // Fill substructure vars with default values when no constituents in jet
208  if (nConstits==0) {
209  dec_Qw(*jet) = -999;
210 
211  dec_Tau1(*jet) = -999;
212  dec_Tau2(*jet) = -999;
213  dec_Tau3(*jet) = -999;
214  dec_Tau4(*jet) = -999;
215 
216  dec_Tau21(*jet) = -999;
217  dec_Tau32(*jet) = -999;
218 
219  dec_Split12(*jet) = -999;
220  dec_Split23(*jet) = -999;
221  dec_Split34(*jet) = -999;
222 
223  dec_ECF1(*jet) = -999;
224  dec_ECF2(*jet) = -999;
225  dec_ECF3(*jet) = -999;
226  dec_ECF4(*jet) = -999;
227 
228  dec_C2(*jet) = -999;
229  dec_D2(*jet) = -999;
230 
231  dec_pT(*jet) = 0;
232  dec_m(*jet) = 0;
233  dec_NConstits(*jet) = 0;
234  dec_eta(*jet) = -999;
235  dec_phi(*jet) = -999;
236 
237  return StatusCode::SUCCESS;
238  }
239 
240  // Save reclustered jet info
241  dec_pT(*jet) = groomed_jet.pt();
242  dec_m(*jet) = groomed_jet.m();
243  dec_NConstits(*jet) = nConstits;
244  dec_eta(*jet) = groomed_jet.eta();
245  dec_phi(*jet) = groomed_jet.phi() - M_PI;
246 
247  // Qw
248  JetSubStructureUtils::Qw qw;
249  dec_Qw(*jet) = qw.result(groomed_jet);
250 
251  // Nsubjetiness
252  fastjet::contrib::WTA_KT_Axes wta_kt_axes;
253  fastjet::contrib::NormalizedCutoffMeasure normalized_measure(1.0, 1.0, 1000000);
254  JetSubStructureUtils::Nsubjettiness Tau1_wta(1, wta_kt_axes, normalized_measure);
255  JetSubStructureUtils::Nsubjettiness Tau2_wta(2, wta_kt_axes, normalized_measure);
256  JetSubStructureUtils::Nsubjettiness Tau3_wta(3, wta_kt_axes, normalized_measure);
257  JetSubStructureUtils::Nsubjettiness Tau4_wta(4, wta_kt_axes, normalized_measure);
258  float tau1_wta = Tau1_wta.result(groomed_jet);
259  float tau2_wta = Tau2_wta.result(groomed_jet);
260  float tau3_wta = Tau3_wta.result(groomed_jet);
261  float tau4_wta = Tau4_wta.result(groomed_jet);
262  float tau21_wta = -999;
263  float tau32_wta = -999;
264 
265  if( tau1_wta > 1e-8 ) {
266  tau21_wta = tau2_wta / tau1_wta;
267  }
268  if( tau2_wta > 1e-8 ) {
269  tau32_wta = tau3_wta / tau2_wta;
270  }
271 
272  dec_Tau1(*jet) = tau1_wta;
273  dec_Tau2(*jet) = tau2_wta;
274  dec_Tau3(*jet) = tau3_wta;
275  dec_Tau4(*jet) = tau4_wta;
276  dec_Tau21(*jet) = tau21_wta;
277  dec_Tau32(*jet) = tau32_wta;
278 
279  // KtSplittingScale
280  JetSubStructureUtils::KtSplittingScale Split12(1);
281  JetSubStructureUtils::KtSplittingScale Split23(2);
282  JetSubStructureUtils::KtSplittingScale Split34(3);
283  float split12 = Split12.result(groomed_jet);
284  float split23 = Split23.result(groomed_jet);
285  float split34 = Split34.result(groomed_jet);
286 
287  dec_Split12(*jet) = split12;
288  dec_Split23(*jet) = split23;
289  dec_Split34(*jet) = split34;
290 
291  // EnergyCorrelator
292  JetSubStructureUtils::EnergyCorrelator ECF1(1, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
293  JetSubStructureUtils::EnergyCorrelator ECF2(2, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
294  JetSubStructureUtils::EnergyCorrelator ECF3(3, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
295  JetSubStructureUtils::EnergyCorrelator ECF4(4, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
296  float ecf1 = ECF1.result(groomed_jet);
297  float ecf2 = ECF2.result(groomed_jet);
298  float ecf3 = ECF3.result(groomed_jet);
299  float ecf4 = ECF4.result(groomed_jet);
300  float c2 = -999;
301  float d2 = -999;
302 
303  if( ecf2 > 1e-8 ) {
304  c2 = ecf3 * ecf1 / pow( ecf2, 2.0 );
305  d2 = ecf3 * pow( ecf1, 3.0 ) / pow( ecf2, 3.0 );
306  }
307 
308  dec_ECF1(*jet) = ecf1;
309  dec_ECF2(*jet) = ecf2;
310  dec_ECF3(*jet) = ecf3;
311  dec_ECF4(*jet) = ecf4;
312  dec_C2(*jet) = c2;
313  dec_D2(*jet) = d2;
314 
315  }
316 
317  return StatusCode::SUCCESS;
318 }
319 
320 
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:210
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:731
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 override
Definition: RCJetSubstructureAug.cxx:88
python.StandardJetMods.qw
qw
Definition: StandardJetMods.py:302
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
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:73
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 Wed Apr 23 2025 21:17:39 for ATLAS Offline Software by   doxygen 1.8.18