ATLAS Offline Software
RCJetSubstructureAug.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 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 <vector>
18 
19 #include "fastjet/PseudoJet.hh"
20 #include "fastjet/ClusterSequence.hh"
21 #include "fastjet/JetDefinition.hh"
22 #include "fastjet/contrib/SoftDrop.hh"
23 #include "fastjet/tools/Filter.hh"
24 #include "fastjet/contrib/Nsubjettiness.hh"
25 #include "JetSubStructureUtils/Qw.h"
26 #include "JetSubStructureUtils/Nsubjettiness.h"
27 #include "JetSubStructureUtils/KtSplittingScale.h"
28 #include "JetSubStructureUtils/EnergyCorrelator.h"
29 
30 // Constructor
31 DerivationFramework::RCJetSubstructureAug::RCJetSubstructureAug(const std::string& t, const std::string& n, const IInterface* p) :
32 base_class(t, n, p) {}
33 
34 // Destructor
35 DerivationFramework::RCJetSubstructureAug::~RCJetSubstructureAug()
36 {}
37 
38 // Athena initialize and finalize
39 StatusCode DerivationFramework::RCJetSubstructureAug::initialize()
40 {
41  ATH_MSG_VERBOSE("initialize() ...");
42  if (m_jetKey.key().empty()) {
43  ATH_MSG_FATAL("No jet collection provided for augmentation.");
44  return StatusCode::FAILURE;
45  }
46  ATH_CHECK( m_jetKey.initialize() );
47 
48  // Set up the text-parsing machinery for selectiong the jet directly according to user cuts
49  if (!m_selectionString.empty()) {
50  ATH_CHECK( initializeParser( m_selectionString ));
51  }
52 
53  // Init moment struct
54  m_moments = new moments_t(m_suffix);
55 
56  // Define the trimmer
57  m_trimmer = std::make_unique<fastjet::Filter>(fastjet::JetDefinition(fastjet::kt_algorithm, m_rclus), fastjet::SelectorPtFractionMin(m_ptfrac));
58  m_softdropper = std::make_unique<fastjet::contrib::SoftDrop>(m_beta, m_zcut, m_R0);
59 
60  return StatusCode::SUCCESS;
61 }
62 
63 StatusCode DerivationFramework::RCJetSubstructureAug::finalize()
64 {
65  ATH_MSG_VERBOSE("finalize() ...");
66  delete m_moments;
67  return StatusCode::SUCCESS;
68 }
69 
70 StatusCode DerivationFramework::RCJetSubstructureAug::addBranches() const
71 {
72  const EventContext& ctx = Gaudi::Hive::currentContext();
73 
74  SG::ReadHandle<xAOD::JetContainer> jets(m_jetKey,ctx);
75  if (!jets.isValid()) {
76  ATH_MSG_ERROR("No jet collection with name " << m_jetKey.key() << " found in StoreGate!");
77  return StatusCode::FAILURE;
78  }
79  unsigned int nJets(jets->size());
80  std::vector<const xAOD::Jet*> jetToCheck; jetToCheck.clear();
81 
82  // Execute the text parser if requested
83  if (!m_selectionString.empty()) {
84  std::vector<int> entries = m_parser->evaluateAsVector();
85  unsigned int nEntries = entries.size();
86  // check the sizes are compatible
87  if (nJets != nEntries ) {
88  ATH_MSG_ERROR("Sizes incompatible! Are you sure your selection string used jets??");
89  return StatusCode::FAILURE;
90  } else {
91  // identify which jets to keep for the thinning check
92  for (unsigned int i=0; i<nJets; ++i) if (entries[i]==1) jetToCheck.push_back((*jets)[i]);
93  }
94  } else {
95  for (unsigned int i=0; i<nJets; ++i) jetToCheck.push_back((*jets)[i]);
96  }
97 
98  std::vector<const xAOD::IParticle*> mergedGhostConstits;
99  std::vector<const xAOD::IParticle*> ghosts;
100  std::vector<fastjet::PseudoJet> constituents;
101 
102  for (const xAOD::Jet* jet : jetToCheck) {
103 
104  // Get list of ghost constituents
105  mergedGhostConstits.clear();
106  for (std::string ghostName: m_ghostNames) {
107  ghosts.clear();
108  ATH_CHECK(jet->getAssociatedObjects(ghostName, ghosts));
109  mergedGhostConstits.insert(mergedGhostConstits.end(), ghosts.begin(), ghosts.end());
110  }
111 
112  // Construct list of constituent PseudoJets from constituents and compute timing information
113  constituents.clear();
114  double eTot = 0;
115  double time = 0;
116  for (auto constit : mergedGhostConstits){
117  if (constit==nullptr) {
118  ATH_MSG_INFO("Invalid constituent link, skipping ...");
119  continue;
120  }
121  // Filter constituents on negative energy
122  if (constit->e()<0) {ATH_MSG_INFO("################## Negative energy cluster"); continue;}
123  constituents.push_back( fastjet::PseudoJet(constit->p4()) );
124 
125  // Timing info of constituent only for calo constituents
126  if (constit->type() == xAOD::Type::CaloCluster) {
127  auto caloConstit = dynamic_cast<const xAOD::CaloCluster*> (constit);
128  double eConstit = caloConstit->e()* caloConstit->e();
129  time += caloConstit->time()* eConstit;
130  eTot += eConstit;
131  }
132  }
133 
134  // Fill timing info
135  if (eTot==0) {
136  m_moments->dec_timing(*jet) = 0;
137  } else {
138  m_moments->dec_timing(*jet) = time/eTot;
139  }
140 
141  // Get number of constituents
142  int nConstits = constituents.size();
143 
144  // Run clustering on constituents if not empty
145  fastjet::PseudoJet groomed_jet;
146  if (nConstits!=0) {
147  auto jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm, 1.5);
148  fastjet::ClusterSequence cs(constituents, jet_def);
149  fastjet::PseudoJet recluster_jet = cs.inclusive_jets(0.0).front();
150 
151  // Apply grooming to reclustered jet
152  if (m_grooming=="Trimming"){
153  groomed_jet = m_trimmer->result(recluster_jet);
154  } else if (m_grooming=="SoftDrop"){
155  groomed_jet = m_softdropper->result(recluster_jet);
156  } else {
157  ATH_MSG_DEBUG(" No grooming requested or wrong one, will not apply grooming");
158  groomed_jet = recluster_jet;
159  }
160 
161  // update nConstit
162  nConstits = groomed_jet.constituents().size();
163  }
164 
165  // Fill substructure vars with default values when no constituents in jet
166  if (nConstits==0) {
167  m_moments->dec_Qw(*jet) = -999;
168 
169  m_moments->dec_Tau1(*jet) = -999;
170  m_moments->dec_Tau2(*jet) = -999;
171  m_moments->dec_Tau3(*jet) = -999;
172  m_moments->dec_Tau4(*jet) = -999;
173 
174  m_moments->dec_Tau21(*jet) = -999;
175  m_moments->dec_Tau32(*jet) = -999;
176 
177  m_moments->dec_Split12(*jet) = -999;
178  m_moments->dec_Split23(*jet) = -999;
179  m_moments->dec_Split34(*jet) = -999;
180 
181  m_moments->dec_ECF1(*jet) = -999;
182  m_moments->dec_ECF2(*jet) = -999;
183  m_moments->dec_ECF3(*jet) = -999;
184  m_moments->dec_ECF4(*jet) = -999;
185 
186  m_moments->dec_C2(*jet) = -999;
187  m_moments->dec_D2(*jet) = -999;
188 
189  m_moments->dec_pT(*jet) = 0;
190  m_moments->dec_m(*jet) = 0;
191  m_moments->dec_NConstits(*jet) = 0;
192  m_moments->dec_eta(*jet) = -999;
193  m_moments->dec_phi(*jet) = -999;
194 
195  return StatusCode::SUCCESS;
196  }
197 
198  // Save reclustered jet infos
199  m_moments->dec_pT(*jet) = groomed_jet.pt();
200  m_moments->dec_m(*jet) = groomed_jet.m();
201  m_moments->dec_NConstits(*jet) = nConstits;
202  m_moments->dec_eta(*jet) = groomed_jet.eta();
203  m_moments->dec_phi(*jet) = groomed_jet.phi() - M_PI;
204 
205  // Qw
206  JetSubStructureUtils::Qw qw;
207  m_moments->dec_Qw(*jet) = qw.result(groomed_jet);
208 
209  // Nsubjetiness
210  fastjet::contrib::WTA_KT_Axes wta_kt_axes;
211  fastjet::contrib::NormalizedCutoffMeasure normalized_measure(1.0, 1.0, 1000000);
212  JetSubStructureUtils::Nsubjettiness Tau1_wta(1, wta_kt_axes, normalized_measure);
213  JetSubStructureUtils::Nsubjettiness Tau2_wta(2, wta_kt_axes, normalized_measure);
214  JetSubStructureUtils::Nsubjettiness Tau3_wta(3, wta_kt_axes, normalized_measure);
215  JetSubStructureUtils::Nsubjettiness Tau4_wta(4, wta_kt_axes, normalized_measure);
216  float tau1_wta = Tau1_wta.result(groomed_jet);
217  float tau2_wta = Tau2_wta.result(groomed_jet);
218  float tau3_wta = Tau3_wta.result(groomed_jet);
219  float tau4_wta = Tau4_wta.result(groomed_jet);
220  float tau21_wta = -999;
221  float tau32_wta = -999;
222 
223  if( tau1_wta > 1e-8 ) {
224  tau21_wta = tau2_wta / tau1_wta;
225  }
226  if( tau2_wta > 1e-8 ) {
227  tau32_wta = tau3_wta / tau2_wta;
228  }
229 
230  m_moments->dec_Tau1(*jet) = tau1_wta;
231  m_moments->dec_Tau2(*jet) = tau2_wta;
232  m_moments->dec_Tau3(*jet) = tau3_wta;
233  m_moments->dec_Tau4(*jet) = tau4_wta;
234  m_moments->dec_Tau21(*jet) = tau21_wta;
235  m_moments->dec_Tau32(*jet) = tau32_wta;
236 
237  // KtSplittingScale
238  JetSubStructureUtils::KtSplittingScale Split12(1);
239  JetSubStructureUtils::KtSplittingScale Split23(2);
240  JetSubStructureUtils::KtSplittingScale Split34(3);
241  float split12 = Split12.result(groomed_jet);
242  float split23 = Split23.result(groomed_jet);
243  float split34 = Split34.result(groomed_jet);
244 
245  m_moments->dec_Split12(*jet) = split12;
246  m_moments->dec_Split23(*jet) = split23;
247  m_moments->dec_Split34(*jet) = split34;
248 
249  // EnergyCorrelator
250  JetSubStructureUtils::EnergyCorrelator ECF1(1, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
251  JetSubStructureUtils::EnergyCorrelator ECF2(2, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
252  JetSubStructureUtils::EnergyCorrelator ECF3(3, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
253  JetSubStructureUtils::EnergyCorrelator ECF4(4, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
254  float ecf1 = ECF1.result(groomed_jet);
255  float ecf2 = ECF2.result(groomed_jet);
256  float ecf3 = ECF3.result(groomed_jet);
257  float ecf4 = ECF4.result(groomed_jet);
258  float c2 = -999;
259  float d2 = -999;
260 
261  if( ecf2 > 1e-8 ) {
262  c2 = ecf3 * ecf1 / pow( ecf2, 2.0 );
263  d2 = ecf3 * pow( ecf1, 3.0 ) / pow( ecf2, 3.0 );
264  }
265 
266  m_moments->dec_ECF1(*jet) = ecf1;
267  m_moments->dec_ECF2(*jet) = ecf2;
268  m_moments->dec_ECF3(*jet) = ecf3;
269  m_moments->dec_ECF4(*jet) = ecf4;
270  m_moments->dec_C2(*jet) = c2;
271  m_moments->dec_D2(*jet) = d2;
272 
273  }
274 
275  return StatusCode::SUCCESS;
276 }
277 
278 
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
DerivationFramework::RCJetSubstructureAug::finalize
virtual StatusCode finalize() override
Definition: RCJetSubstructureAug.cxx:63
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:70
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
M_PI
#define M_PI
Definition: ActiveFraction.h:11
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:116
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:59
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:39
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
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Nsubjettiness.h
DerivationFramework::RCJetSubstructureAug::addBranches
virtual StatusCode addBranches() const override
Definition: RCJetSubstructureAug.cxx:70
DerivationFramework::RCJetSubstructureAug::moments_t
Definition: RCJetSubstructureAug.h:73
python.StandardJetMods.qw
qw
Definition: StandardJetMods.py:286
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:31
DerivationFramework::RCJetSubstructureAug::~RCJetSubstructureAug
virtual ~RCJetSubstructureAug()
Definition: RCJetSubstructureAug.cxx:35
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
defineDB.jets
list jets
Definition: JetTagCalibration/share/defineDB.py:24
dqBeamSpot.nEntries
int nEntries
Definition: dqBeamSpot.py:73
EnergyCorrelator.h
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
RCJetSubstructureAug.h
Generated on Thu Nov 7 2024 21:24:26 for ATLAS Offline Software by   doxygen 1.8.18