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