ATLAS Offline Software
NSubjettinessTool.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 #include "JetSubStructureMomentTools/NSubjettinessTool.h"
6 #include "JetSubStructureUtils/Nsubjettiness.h"
7 #include "AthContainers/ConstAccessor.h"
8 #include "AsgDataHandles/WriteDecorHandle.h"
9 #include "CxxUtils/ubsan_suppress.h"
10 
11 #include "fastjet/contrib/Nsubjettiness.hh"
12 #include "fastjet/contrib/AxesDefinition.hh"
13 
14 NSubjettinessTool::NSubjettinessTool(const std::string& name) :
15  JetSubStructureMomentToolsBase(name)
16 {
17  declareProperty("Alpha", m_Alpha = 1.0);
18  declareProperty("AlphaList", m_rawAlphaVals = {});
19  declareProperty("DoDichroic", m_doDichroic = false);
20 }
21 
22 StatusCode NSubjettinessTool::initialize() {
23  if(m_jetContainerName.empty()){
24  ATH_MSG_ERROR("NSubjettinessTool needs to have its input jet container name configured!");
25  return StatusCode::FAILURE;
26  }
27 
29  ATH_CHECK( JetSubStructureMomentToolsBase::initialize() );
30 
32  m_moments.emplace_back( 1.0, moments_t(1.0, m_prefix) );
33 
35  if( std::abs(m_Alpha-1.0) > 1.0e-5 ) {
36 
38  ATH_MSG_WARNING( "The Alpha property is deprecated, please use the AlphaList property to provide a list of values" );
39 
41  m_moments.emplace_back( m_Alpha, moments_t(m_Alpha, m_prefix) );
42 
43  }
44 
46  for( float alpha : m_rawAlphaVals ) {
47 
49  float alphaFix = round( alpha * 10.0 ) / 10.0;
50  if( std::abs(alpha-alphaFix) > 1.0e-5 ) ATH_MSG_DEBUG( "alpha = " << alpha << " has been rounded to " << alphaFix );
51 
53  if( alphaFix < 0.0 ) {
54  ATH_MSG_WARNING( "alpha must be positive. Skipping alpha = " << alpha );
55  continue;
56  }
57 
59  m_moments.emplace_back( alphaFix, moments_t(alphaFix, m_prefix) );
60 
61  }
62 
63  for( auto const& moment : m_moments ) {
64  ATH_MSG_DEBUG( "Including alpha = " << moment.first );
65  }
66 
67  // Initialise WriteDecorHandleKeys
68  for( const auto& [alpha, moment] : m_moments ) {
69  m_Tau1_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
70  "Tau1" + moment.suffix);
71  m_Tau2_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
72  "Tau2" + moment.suffix);
73  m_Tau3_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
74  "Tau3" + moment.suffix);
75  m_Tau4_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
76  "Tau4" + moment.suffix);
77 
78  m_Tau2_ungroomed_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
79  "Tau2_ungroomed" + moment.suffix);
80  m_Tau3_ungroomed_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
81  "Tau3_ungroomed" + moment.suffix);
82  m_Tau4_ungroomed_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
83  "Tau4_ungroomed" + moment.suffix);
84 
85  m_Tau1_wta_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
86  "Tau1_wta" + moment.suffix);
87  m_Tau2_wta_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
88  "Tau2_wta" + moment.suffix);
89  m_Tau3_wta_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
90  "Tau3_wta" + moment.suffix);
91  m_Tau4_wta_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
92  "Tau4_wta" + moment.suffix);
93 
94  m_Tau2_wta_ungroomed_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
95  "Tau2_wta_ungroomed" + moment.suffix);
96  m_Tau3_wta_ungroomed_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
97  "Tau3_wta_ungroomed" + moment.suffix);
98  m_Tau4_wta_ungroomed_Keys.emplace_back(m_jetContainerName + "." + moment.prefix +
99  "Tau4_wta_ungroomed" + moment.suffix);
100  }
101 
102  ATH_CHECK(m_Tau1_Keys.initialize());
103  ATH_CHECK(m_Tau2_Keys.initialize());
104  ATH_CHECK(m_Tau3_Keys.initialize());
105  ATH_CHECK(m_Tau4_Keys.initialize());
106 
107  ATH_CHECK(m_Tau2_ungroomed_Keys.initialize());
108  ATH_CHECK(m_Tau3_ungroomed_Keys.initialize());
109  ATH_CHECK(m_Tau4_ungroomed_Keys.initialize());
110 
111  ATH_CHECK(m_Tau1_wta_Keys.initialize());
112  ATH_CHECK(m_Tau2_wta_Keys.initialize());
113  ATH_CHECK(m_Tau3_wta_Keys.initialize());
114  ATH_CHECK(m_Tau4_wta_Keys.initialize());
115 
116  ATH_CHECK(m_Tau2_wta_ungroomed_Keys.initialize());
117  ATH_CHECK(m_Tau3_wta_ungroomed_Keys.initialize());
118  ATH_CHECK(m_Tau4_wta_ungroomed_Keys.initialize());
119 
120  return StatusCode::SUCCESS;
121 }
122 
123 StatusCode NSubjettinessTool::modify(xAOD::JetContainer& jets) const {
124 
125  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau1;
126  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau2;
127  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau3;
128  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau4;
129 
130  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau2_ungroomed;
131  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau3_ungroomed;
132  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau4_ungroomed;
133 
134  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau1_wta;
135  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau2_wta;
136  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau3_wta;
137  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau4_wta;
138 
139  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau2_wta_ungroomed;
140  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau3_wta_ungroomed;
141  std::vector<SG::WriteDecorHandle<xAOD::JetContainer, float>> wdhs_Tau4_wta_ungroomed;
142 
143  for(unsigned int i=0; i<m_moments.size(); i++) {
144  wdhs_Tau1.emplace_back(m_Tau1_Keys[i]);
145  wdhs_Tau2.emplace_back(m_Tau2_Keys[i]);
146  wdhs_Tau3.emplace_back(m_Tau3_Keys[i]);
147  wdhs_Tau4.emplace_back(m_Tau4_Keys[i]);
148 
149  wdhs_Tau2_ungroomed.emplace_back(m_Tau2_ungroomed_Keys[i]);
150  wdhs_Tau3_ungroomed.emplace_back(m_Tau3_ungroomed_Keys[i]);
151  wdhs_Tau4_ungroomed.emplace_back(m_Tau4_ungroomed_Keys[i]);
152 
153  wdhs_Tau1_wta.emplace_back(m_Tau1_wta_Keys[i]);
154  wdhs_Tau2_wta.emplace_back(m_Tau2_wta_Keys[i]);
155  wdhs_Tau3_wta.emplace_back(m_Tau3_wta_Keys[i]);
156  wdhs_Tau4_wta.emplace_back(m_Tau4_wta_Keys[i]);
157 
158  wdhs_Tau2_wta_ungroomed.emplace_back(m_Tau2_wta_ungroomed_Keys[i]);
159  wdhs_Tau3_wta_ungroomed.emplace_back(m_Tau3_wta_ungroomed_Keys[i]);
160  wdhs_Tau4_wta_ungroomed.emplace_back(m_Tau4_wta_ungroomed_Keys[i]);
161  }
162 
163  for(const xAOD::Jet* injet : jets){
164  fastjet::PseudoJet jet;
165  fastjet::PseudoJet jet_ungroomed;
166 
168  bool calculate = SetupDecoration(jet, *injet);
169 
171  bool calculate_ungroomed = false;
172 
173  if( m_doDichroic ) {
174 
176  static const SG::ConstAccessor<ElementLink<xAOD::JetContainer> > ParentAcc ("Parent");
177  ElementLink<xAOD::JetContainer> parentLink = ParentAcc (*injet);
178 
180  if( !parentLink.isValid() ) {
181  ATH_MSG_ERROR( "Parent element link is not valid. Aborting" );
182  return StatusCode::FAILURE;
183  }
184 
185  const xAOD::Jet* parentJet = *(parentLink);
186  calculate_ungroomed = SetupDecoration(jet_ungroomed,*parentJet);
187  }
188 
189  // Supress a warning about undefined behavior in the fastjet
190  // WTA_KT_Axes ctor:
191  // .../fastjet/contrib/AxesDefinition.hh:551:43: runtime error: member access within address 0x7ffd770850d0 which does not point to an object of type 'WTA_KT_Axes'
192  // 0x7ffd770850d0: note: object has invalid vptr
193  std::once_flag oflag;
194  std::call_once (oflag, CxxUtils::ubsan_suppress,
195  []() { fastjet::contrib::WTA_KT_Axes x; });
196 
197  for(unsigned int i=0; i<m_moments.size(); i++) {
198 
199  float alpha = m_moments[i].first;
200 
202  fastjet::contrib::NormalizedCutoffMeasure normalized_measure
203  (alpha, injet->getSizeParameter(), 1000000);
204 
205  float Tau1_value = -999;
206  float Tau2_value = -999;
207  float Tau3_value = -999;
208  float Tau4_value = -999;
209 
210  float Tau2_ungroomed_value = -999;
211  float Tau3_ungroomed_value = -999;
212  float Tau4_ungroomed_value = -999;
213 
214  float Tau1_wta_value = -999;
215  float Tau2_wta_value = -999;
216  float Tau3_wta_value = -999;
217  float Tau4_wta_value = -999;
218 
219  float Tau2_wta_ungroomed_value = -999;
220  float Tau3_wta_ungroomed_value = -999;
221  float Tau4_wta_ungroomed_value = -999;
222 
223  if( calculate ) {
224 
227  fastjet::contrib::KT_Axes kt_axes;
228  JetSubStructureUtils::Nsubjettiness tau1(1, kt_axes, normalized_measure);
229  JetSubStructureUtils::Nsubjettiness tau2(2, kt_axes, normalized_measure);
230  JetSubStructureUtils::Nsubjettiness tau3(3, kt_axes, normalized_measure);
231  JetSubStructureUtils::Nsubjettiness tau4(4, kt_axes, normalized_measure);
232 
233  Tau1_value = tau1.result(jet);
234  Tau2_value = tau2.result(jet);
235  Tau3_value = tau3.result(jet);
236  Tau4_value = tau4.result(jet);
237 
238  if( calculate_ungroomed ) {
239  Tau2_ungroomed_value = tau2.result(jet_ungroomed);
240  Tau3_ungroomed_value = tau3.result(jet_ungroomed);
241  Tau4_ungroomed_value = tau4.result(jet_ungroomed);
242  }
243 
246  fastjet::contrib::WTA_KT_Axes wta_kt_axes;
247  JetSubStructureUtils::Nsubjettiness tau1_wta(1, wta_kt_axes, normalized_measure);
248  JetSubStructureUtils::Nsubjettiness tau2_wta(2, wta_kt_axes, normalized_measure);
249  JetSubStructureUtils::Nsubjettiness tau3_wta(3, wta_kt_axes, normalized_measure);
250  JetSubStructureUtils::Nsubjettiness tau4_wta(4, wta_kt_axes, normalized_measure);
251 
252  Tau1_wta_value = tau1_wta.result(jet);
253  Tau2_wta_value = tau2_wta.result(jet);
254  Tau3_wta_value = tau3_wta.result(jet);
255  Tau4_wta_value = tau4_wta.result(jet);
256 
257  if( calculate_ungroomed ) {
258  Tau2_wta_ungroomed_value = tau2_wta.result(jet_ungroomed);
259  Tau3_wta_ungroomed_value = tau3_wta.result(jet_ungroomed);
260  Tau4_wta_ungroomed_value = tau4_wta.result(jet_ungroomed);
261  }
262 
263  }
264 
265  wdhs_Tau1[i](*injet) = Tau1_value;
266  wdhs_Tau2[i](*injet) = Tau2_value;
267  wdhs_Tau3[i](*injet) = Tau3_value;
268  wdhs_Tau4[i](*injet) = Tau4_value;
269 
270  wdhs_Tau2_ungroomed[i](*injet) = Tau2_ungroomed_value;
271  wdhs_Tau3_ungroomed[i](*injet) = Tau3_ungroomed_value;
272  wdhs_Tau4_ungroomed[i](*injet) = Tau4_ungroomed_value;
273 
274  wdhs_Tau1_wta[i](*injet) = Tau1_wta_value;
275  wdhs_Tau2_wta[i](*injet) = Tau2_wta_value;
276  wdhs_Tau3_wta[i](*injet) = Tau3_wta_value;
277  wdhs_Tau4_wta[i](*injet) = Tau4_wta_value;
278 
279  wdhs_Tau2_wta_ungroomed[i](*injet) = Tau2_wta_ungroomed_value;
280  wdhs_Tau3_wta_ungroomed[i](*injet) = Tau3_wta_ungroomed_value;
281  wdhs_Tau4_wta_ungroomed[i](*injet) = Tau4_wta_ungroomed_value;
282  }
283 
284  }
285 
286  return StatusCode::SUCCESS;
287 }
NSubjettinessTool::modify
StatusCode modify(xAOD::JetContainer &jets) const override
Loop over calls to modifyJet.
Definition: NSubjettinessTool.cxx:123
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
NSubjettinessTool::m_Tau1_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau1_Keys
Definition: NSubjettinessTool.h:77
NSubjettinessTool::m_Tau4_ungroomed_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau4_ungroomed_Keys
Definition: NSubjettinessTool.h:90
add-xsec-uncert-quadrature-N.alpha
alpha
Definition: add-xsec-uncert-quadrature-N.py:110
ubsan_suppress.h
Helper for suppressing ubsan warnings.
NSubjettinessTool::m_Tau2_ungroomed_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau2_ungroomed_Keys
Definition: NSubjettinessTool.h:86
JetSubStructureUtils::Nsubjettiness::result
virtual double result(const fastjet::PseudoJet &jet) const
Definition: Nsubjettiness.h:21
NSubjettinessTool::m_Tau2_wta_ungroomed_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau2_wta_ungroomed_Keys
Definition: NSubjettinessTool.h:102
MuonGM::round
float round(const float toRound, const unsigned int decimals)
Definition: Mdt.cxx:27
NSubjettinessTool::moments_t
N-subjettiness moments structure.
Definition: NSubjettinessTool.h:53
defineDB.jets
jets
Definition: JetTagCalibration/share/defineDB.py:24
NSubjettinessTool::m_Tau2_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau2_Keys
Definition: NSubjettinessTool.h:79
SG::ConstAccessor
Helper class to provide constant type-safe access to aux data.
Definition: ConstAccessor.h:55
NSubjettinessTool::m_Tau3_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau3_Keys
Definition: NSubjettinessTool.h:81
NSubjettinessTool::m_Tau3_wta_ungroomed_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau3_wta_ungroomed_Keys
Definition: NSubjettinessTool.h:104
x
#define x
jet
Definition: JetCalibTools_PlotJESFactors.cxx:23
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
lumiFormat.i
int i
Definition: lumiFormat.py: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
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition: AthCommonDataStore.h:145
NSubjettinessTool::m_Tau1_wta_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau1_wta_Keys
Definition: NSubjettinessTool.h:93
NSubjettinessTool::m_rawAlphaVals
std::vector< float > m_rawAlphaVals
Definition: NSubjettinessTool.h:71
NSubjettinessTool::initialize
StatusCode initialize() override
Dummy implementation of the initialisation function.
Definition: NSubjettinessTool.cxx:22
NSubjettinessTool::m_Tau4_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau4_Keys
Definition: NSubjettinessTool.h:83
JetSubStructureMomentToolsBase::initialize
StatusCode initialize()
Dummy implementation of the initialisation function.
Definition: JetSubStructureMomentToolsBase.cxx:15
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Nsubjettiness.h
DataVector
Derived DataVector<T>.
Definition: DataVector.h:794
WriteDecorHandle.h
Handle class for adding a decoration to an object.
NSubjettinessTool::m_Tau4_wta_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau4_wta_Keys
Definition: NSubjettinessTool.h:99
NSubjettinessTool::m_Alpha
float m_Alpha
Configurable as properties.
Definition: NSubjettinessTool.h:70
NSubjettinessTool::m_jetContainerName
Gaudi::Property< std::string > m_jetContainerName
Definition: NSubjettinessTool.h:41
EMFourMomBuilder::calculate
void calculate(xAOD::Electron &electron)
Definition: EMFourMomBuilder.cxx:69
JetSubStructureMomentToolsBase::m_prefix
std::string m_prefix
Definition: JetSubStructureMomentToolsBase.h:30
JetSubStructureMomentToolsBase::SetupDecoration
bool SetupDecoration(fastjet::PseudoJet &pseudojet, const xAOD::Jet &jet, bool requireJetStructure=false) const
Definition: JetSubStructureMomentToolsBase.cxx:30
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
NSubjettinessTool::m_Tau4_wta_ungroomed_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau4_wta_ungroomed_Keys
Definition: NSubjettinessTool.h:106
xAOD::Jet_v1
Class describing a jet.
Definition: Jet_v1.h:57
NSubjettinessTool::NSubjettinessTool
NSubjettinessTool(const std::string &name)
Definition: NSubjettinessTool.cxx:14
JetSubStructureUtils::Nsubjettiness
Definition: Nsubjettiness.h:14
NSubjettinessTool::m_moments
std::vector< std::pair< float, moments_t > > m_moments
Map of decorators using alpha as the key.
Definition: NSubjettinessTool.h:75
NSubjettinessTool::m_Tau2_wta_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau2_wta_Keys
Definition: NSubjettinessTool.h:95
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
NSubjettinessTool::m_Tau3_ungroomed_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau3_ungroomed_Keys
Definition: NSubjettinessTool.h:88
JetSubStructureMomentToolsBase
Definition: JetSubStructureMomentToolsBase.h:18
ConstAccessor.h
Helper class to provide constant type-safe access to aux data.
NSubjettinessTool::m_Tau3_wta_Keys
SG::WriteDecorHandleKeyArray< xAOD::JetContainer > m_Tau3_wta_Keys
Definition: NSubjettinessTool.h:97
NSubjettinessTool::m_doDichroic
bool m_doDichroic
Vector of input values before cleaning.
Definition: NSubjettinessTool.h:72
CxxUtils::ubsan_suppress
void ubsan_suppress(void(*func)())
Helper for suppressing ubsan warnings.
Definition: ubsan_suppress.cxx:69
NSubjettinessTool.h
Generated on Tue Sep 2 2025 21:17:58 for ATLAS Offline Software by   doxygen 1.8.18