ATLAS Offline Software
ClassifyAndCalculateHFTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // //
7 // ClassifyAndCalculateHFTool.cxx //
8 // Implementation file for class ClassifyAndCalculateHFTool //
9 // Author: Adrian Berrocal Guardia <adrian.berrocal.guardia@cern.ch> //
10 // //
12 
16 
17 namespace DerivationFramework {
18 
19  /*
20  ---------------------------------------------------------------------------------------------------------------------------------------
21  ------------------------------------------------------- Constructor/Destructor --------------------------------------------------------
22  ---------------------------------------------------------------------------------------------------------------------------------------
23  */
24 
25  ClassifyAndCalculateHFTool::ClassifyAndCalculateHFTool(const std::string& t, const std::string& n, const IInterface* p) : AthAlgTool(t,n,p){
26  }
27 
29  }
30 
31  /*
32  ---------------------------------------------------------------------------------------------------------------------------------------
33  --------------------------------------------------------- Initialize/Finalize ---------------------------------------------------------
34  ---------------------------------------------------------------------------------------------------------------------------------------
35  */
36 
38 
39  ATH_MSG_INFO("Initialize");
40 
41  // Print the string variables.
42 
43  ATH_MSG_INFO("Cut on the pt of the jets: " << m_jetPtCut);
44  ATH_MSG_INFO("Cut on the eta of the jets: " << m_jetEtaCut);
45  ATH_MSG_INFO("Cut on the pt of the leading hadron: " << m_leadingHadronPtCut);
46  ATH_MSG_INFO("Cut on the ratio between the pt of the leading hadron and the pt of its associated jet: " << m_leadingHadronPtRatioCut);
47 
48  return StatusCode::SUCCESS;
49  }
50 
52  return StatusCode::SUCCESS;
53  }
54 
55  /*
56  ---------------------------------------------------------------------------------------------------------------------------------------
57  -------------------------------------------------------------- Flag Jets --------------------------------------------------------------
58  ---------------------------------------------------------------------------------------------------------------------------------------
59  */
60 
62  const std::map<const xAOD::Jet*, std::vector<xAOD::TruthParticleContainer::const_iterator>>& particleMatch,
63  const std::map<const xAOD::TruthParticle*, DerivationFramework::HadronOriginClassifier::HF_id>& hadronMap,
64  const std::string& hfDecorationName) const{
65 
66  SG::Decorator< int > decorator_flav(hfDecorationName + "_flav");
67  SG::Decorator< int > decorator_id(hfDecorationName + "_id");
68  SG::Decorator< int > decorator_count(hfDecorationName + "_count");
69 
70  for(const xAOD::Jet* jet : *jets){
71 
72  // Check if the jet passes the cuts and if it does not, then skip it.
73 
74  if(jet->p4().Pt() < m_jetPtCut || std::abs(jet->p4().Eta()) > m_jetEtaCut) {
75  decorator_flav(*jet) = -999;
76  decorator_id(*jet) = -999;
77  decorator_count(*jet) = -999;
78  continue;
79  }
80 
81  // Create a set of integer variables to save the necessary variables for the HF classifier:
82  // -flav: Flavour of the jet.
83  // -id: Origin of the most initial hadron of the jet.
84  // -count: Number of associated hadrons of the jet that have the same flavour as the jet.
85 
86  int flav=0;
87  int id=0;
88  int count=0;
89 
90  // Create a set of integer variables to save information that is required to compute the necessary variables for the HF classifier.:
91  // -bcount: It contains the number of B-hadrons.
92  // -ccount: It contains the number of C-hadrons.
93  // -bcountcut: It contains the number of B-hadron that pass the cuts.
94  // -ccountcut: It contains the number of C-hadron that pass the cuts.
95  // -bid: The greatest value of the variable "TopHadronOriginFlag" for B-hadrons (most initial B-hadron).
96  // -cid: The greatest value of the variable "TopHadronOriginFlag" for C-hadrons (most initial B-hadron).
97 
98  int bcount=0;
99  int ccount=0;
100  int bcountcut=0;
101  int ccountcut=0;
102 
103  int bid=0;
104  int cid=0;
105 
106  // Get the hadrons associated with the jet that is being considered.
107 
108  auto it = particleMatch.find (jet);
109  if (it != particleMatch.end()) {
110 
112 
113  // Create two integer variables:
114  // -hforigin: It will contain the origin of the hadron if it is a HF hadron. Otherwise, it will be 6.
115  // -pdgId: It will contain the value of the variable "pdgId" of the hadron.
116 
117  int hforigin = 6;
118  int pdgId = (*hf)->pdgId();
119 
120  // Extract the origin of the hadron.
121 
122  auto h_it = hadronMap.find(*hf);
123  if(h_it!=hadronMap.end()){
124  hforigin= static_cast<int>(h_it->second);
125  }
126 
127  // Check if hforigin is 6 and if it is the case, then hadron is not HF and it is skipped.
128 
129  if(6==hforigin) continue;
130 
131  // Compute the ratio between the pt of the hadron and the pt of its associated jet.
132 
133  float ptratio = (*hf)->p4().Pt()/jet->p4().Pt();
134 
135  // Determine if the hadron is a B-hadron or a C-hadron.
136 
137  int hftype = 0;
138 
139  if(MC::isCharmHadron(pdgId)) hftype=4; // B-hadron
140  if(MC::isBottomHadron(pdgId)) hftype=5; // C-hadron.
141 
142  // Check if hftype is 4 or 5.
143 
144  if(5==hftype){
145 
146  // In this case, hftype is 5 so the hadron is a B-hadron.
147  // Save hforigin in bid if it is greater than the current bid.
148 
149  if(bid<hforigin)bid=hforigin;
150 
151  // Add one to bcount and to bcountcut if hadron passes the cuts.
152 
153  ++bcount;
154 
155  if((*hf)->p4().Pt()>m_leadingHadronPtCut && ptratio>m_leadingHadronPtRatioCut){
156  ++bcountcut;
157  }
158  }
159  else if(4==hftype){
160 
161  // In this case, hftype is 4 so the hadron is a C-hadron.
162  // Save hforigin in cid if it is greater than the current cid.
163 
164  if(cid>hforigin)cid=hforigin;
165 
166  // Add one to ccount and to ccountcut if hadron passes the cuts.
167 
168  ++ccount;
169 
170  if((*hf)->p4().Pt()>m_leadingHadronPtCut && ptratio>m_leadingHadronPtRatioCut){
171  ++ccountcut;
172  }
173 
174  }
175  else{
176 
177  // In this case, hftype is not 4 neither 5 so print an error.
178 
179  ATH_MSG_WARNING("Hadron type '" << hftype << "' is not 4 or 5");
180 
181  }
182 
183  }
184  }
185 
186  // Check if there is at least one B-hadron or a C-hadron that passes the cuts.
187 
188  if(bcountcut){
189 
190  // In this case, at least one B-hadron passes the cuts.
191  // The "flavour" of the jet (flav) is set to 5.
192  // As a id of the jet, the greatest value of the variable "TopHadronOriginFlag" for B-hadrons is used (origin of the most initial hadron).
193  // The number of B-hadrons is saved in count.
194 
195  flav=5;
196  id=bid;
197  count=bcount;
198  }
199  else if(ccountcut){
200 
201  // In this case, no B-hadron passes the cuts, but at least one C-hadron does.
202  // The "flavour" of the jet (flav) is set to 4.
203  // As a id of the jet, the greatest value of the variable "TopHadronOriginFlag" for C-hadrons is used (origin of the most initial hadron).
204  // The number of C-hadrons is saved in count.
205 
206  flav=4;
207  id=cid;
208  count=ccount;
209  }
210 
211  // Dectorate the jet with the flav, id and count.
212  decorator_flav(*jet) = flav;
213  decorator_id(*jet) = id;
214  decorator_count(*jet) = count;
215  }
216 
217  }
218 
219  /*
220  ---------------------------------------------------------------------------------------------------------------------------------------
221  ---------------------------------------------------------- HF Classification ----------------------------------------------------------
222  ---------------------------------------------------------------------------------------------------------------------------------------
223  */
224 
225  int ClassifyAndCalculateHFTool::computeHFClassification(const xAOD::JetContainer* jets, const std::string& hfDecorationName) const{
226 
227  // Create a set of integer variables to save information that is required to compute the HF classifier.:
228  // -b: Number of jets that has just one B-hadron that passes the cuts.
229  // -B: Number of jets that has more than one B-hadron and at least one of them passes the cuts.
230  // -c: Number of jets that has just one C-hadron that passes the cuts (and no B-hadron)
231  // -C: Number of jets that has more than one C-hadron and at least one of them passes the cuts (and no B-hadron).
232  // -b_prompt: Number of jets that has just one prompt B-hadron that passes the cuts.
233  // -B_prompt: Number of jets that has more than one prompt B-hadron and at least one of them passes the cuts.
234  // -c_prompt: Number of jets that has just one prompt C-hadron that passes the cuts (and no B-hadron).
235  // -C_prompt: Number of jets that has more than one prompt C-hadron and at least one of them passes the cuts (and no B-hadron).
236  // -mpifsr_code: HF classifier for events with no prompt additional hadrons, just Multi-Parton Interaction (MPI) and/or Final State Radiation (FSR) hadrons.
237  // Note: prompt hadrons in this context are hadrons that do not come from the direct decay of top, W or Higgs but they are not from MPI or FSR.
238 
239  int b=0, B=0, c=0, C=0;
240  int b_prompt=0, B_prompt=0, c_prompt=0, C_prompt=0;
241 
242  int mpifsr_code=0;
243 
244  for(const xAOD::Jet* jet : *jets){
245 
246  // Check if the jet passes the cuts and if it does not, then skip it.
247 
248  if(jet->p4().Pt() < m_jetPtCut) continue;
249  if(std::abs(jet->p4().Eta()) > m_jetEtaCut) continue;
250 
251  // Get the flavour, the id and the number of hadrons of the considered jet.
252 
253  int flav = 0;
254  int id = 0;
255  int count = 0;
256 
257  SG::ConstAccessor<int> hfflavAcc(hfDecorationName + "_flav");
258  if(hfflavAcc.isAvailable(*jet)){
259  flav=hfflavAcc(*jet);
260  }else{
261  ATH_MSG_WARNING("variable '" + hfDecorationName + "_flav' not found.");
262  continue;
263  }
264 
265  SG::ConstAccessor<int> hfidAcc(hfDecorationName + "_id");
266  if(hfidAcc.isAvailable(*jet)){
267  id=hfidAcc(*jet);
268  }else{
269  ATH_MSG_WARNING("variable '" + hfDecorationName + "_id' not found.");
270  continue;
271  }
272 
273  SG::ConstAccessor<int> hfcountAcc(hfDecorationName + "_count");
274  if(hfcountAcc.isAvailable(*jet)){
275  count=hfcountAcc(*jet);
276  }else{
277  ATH_MSG_WARNING("variable '" + hfDecorationName + "_count' not found.");
278  continue;
279  }
280 
281  // Check the flavour and the id of the jet.
282 
283  if(flav==5 && id < 3){
284 
285  // In this case, the flavour is 5 and id is less than 3.
286  // This means that the jet has at least one B-hadron that is not from direct decay of top, W or Higgs.
287  // Hence, the jet is an additional one.
288 
289  // Check the number of B-hadrons.
290 
291  if(count > 1){
292 
293  // In this case, there is more than one B-hadron, so add 1 to B.
294 
295  B++;
296 
297  }
298  else{
299 
300  // In this case, there is just one B-hadron, so add 1 to b.
301 
302  b++;
303 
304  }
305  }
306  if(flav==4 && (id==0 || id==-1 || id==-2)){
307 
308  // In this case, the flavour is 4 and id is 0, -1 or -2.
309  // This means that the jet has no B-hadron and at least one C-hadron that is not from direct decay of top, W or Higgs.
310  // Hence, the jet is an additional one.
311 
312  // Check the number of C-hadrons.
313 
314  if(count > 1){
315 
316  // In this case, there is more than one C-hadron, so add 1 to C.
317 
318  C++;
319  }
320  else{
321 
322  // In this case, there is just one C-hadron, so add 1 to c.
323 
324  c++;
325  }
326  }
327 
328  // Check the flavour and the id of the jet but considering only prompt particles (id=0).
329 
330  if(flav==5 && id==0){
331 
332  // In this case, the flavour is 5 and id is 0.
333  // This means that the jet has at least one B-hadron that is not from direct decay of top, W or Higgs neither from MPI or FSR.
334 
335  // Check the number of B-hadrons.
336 
337  if(count > 1){
338 
339  // In this case, there is more than one B-hadron, so add 1 to B_prompt.
340 
341  B_prompt++;
342 
343  }
344  else{
345 
346  // In this case, there is jut one B-hadron, so add 1 to b_prompt.
347 
348  b_prompt++;
349 
350  }
351  }
352  if(flav==4 && id==0){
353 
354  // In this case, the flavour is 4 and id is 0.
355  // This means that the jet has no B-hadron and at least one C-hadron that is not from direct decay of top, W or Higgs neither from MPI or FSR.
356 
357  // Check the number of C-hadrons.
358 
359  if(count > 1){
360 
361  // In this case, there is more than one C-hadron, so add 1 to C_prompt.
362 
363  C_prompt++;
364  }
365  else{
366 
367  // In this case, there is just one C-hadron, so add 1 to C_prompt.
368 
369  c_prompt++;
370 
371  }
372  }
373 
374  // In the case when there is no prompt hadrons, then the HF classifier is computed with non-promt hadrons.
375  // This hadrons come from Multi-Parton Interactions (MPI) and the Final State Radiation (FSR).
376  // Compute mpifsr_code which will contain the HF classifier when there is no prompt hadrons.
377  // Each jet adds a different quantity to mpifsr_code depending on its flavour and its id.
378 
379  if(id==1 && flav==5){ // b MPI
380  mpifsr_code -= 1000;
381  }
382  else if(id==2 && flav==5){ // b FSR
383  mpifsr_code -= 100;
384  }
385  else if(id==-1 && flav==4){ // c MPI
386  mpifsr_code -= 10;
387  }
388  else if(id==-2 && flav==4) { // c FSR
389  mpifsr_code -= 1;
390  }
391  }
392 
393  // Compute the ext_code using the number of additional jets with HF hadrons.
394  // Compute the prompt_code using the number of additional jets with propmt HF hadrons.
395 
396  int ext_code = 1000*b+100*B+10*c+1*C;
397  int prompt_code = 1000*b_prompt+100*B_prompt+10*c_prompt+1*C_prompt;
398 
399  // Check if there are prompt hadrons and non-prompt hadrons using ext_code and prompt_code.
400 
401  if(prompt_code==0 && ext_code!=0){
402 
403  // In this case, there is no prompt hadrons but there are MPI and FSR hadrons.
404  // Hence, return mpifsr_code as a HF classifier, which is equal to ext_code but in negative sign.
405 
406  return mpifsr_code;
407 
408  }
409 
410  // In this case, there are prompt hadrons, so return ext_code as HF classifier.
411 
412  return ext_code;
413 
414  }
415 
417 
418  // Check the value of the HF classifier (hfclassif) which is computed with the function computeHFClassification.
419 
420  if(std::abs(hfclassif)>=100){
421 
422  // If the absolute value of hfclassif is greater than 100, then there is at least one jet with a B-hadron.
423  // In this case, return 1.
424 
425  return 1;
426 
427 
428  }
429  else if(hfclassif==0){
430 
431  // If hfclassif is 0, then there is no jet with a B-hadron or a C-hadron.
432  // In this case, return 0.
433 
434  return 0;
435  }
436 
437  // If the absolute value of hfclassif is lower than 100 and non-zero, then there is no jet with a B-hadron but at least one with a C-hadron.
438  // In this case, return -1.
439 
440  return -1;
441  }
442 }
DerivationFramework::ClassifyAndCalculateHFTool::flagJets
void flagJets(const xAOD::JetContainer *jets, const std::map< const xAOD::Jet *, std::vector< xAOD::TruthParticleContainer::const_iterator >> &particleMatch, const std::map< const xAOD::TruthParticle *, DerivationFramework::HadronOriginClassifier::HF_id > &hadronMap, const std::string &hfDecorationName) const
Definition: ClassifyAndCalculateHFTool.cxx:61
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
skel.it
it
Definition: skel.GENtoEVGEN.py:407
defineDB.jets
jets
Definition: JetTagCalibration/share/defineDB.py:24
DerivationFramework::ClassifyAndCalculateHFTool::m_leadingHadronPtCut
Gaudi::Property< float > m_leadingHadronPtCut
Definition: ClassifyAndCalculateHFTool.h:98
SG::ConstAccessor< int >
DerivationFramework::ClassifyAndCalculateHFTool::~ClassifyAndCalculateHFTool
virtual ~ClassifyAndCalculateHFTool()
Definition: ClassifyAndCalculateHFTool.cxx:28
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
DerivationFramework::ClassifyAndCalculateHFTool::initialize
virtual StatusCode initialize() override
Definition: ClassifyAndCalculateHFTool.cxx:37
isBottomHadron
bool isBottomHadron(const T &p)
Definition: AtlasPID.h:901
XMLtoHeader.count
count
Definition: XMLtoHeader.py:84
DerivationFramework::ClassifyAndCalculateHFTool::computeHFClassification
int computeHFClassification(const xAOD::JetContainer *jets, const std::string &hfDecorationName) const
Definition: ClassifyAndCalculateHFTool.cxx:225
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
jet
Definition: JetCalibTools_PlotJESFactors.cxx:23
DerivationFramework::ClassifyAndCalculateHFTool::m_jetPtCut
Gaudi::Property< float > m_jetPtCut
Definition: ClassifyAndCalculateHFTool.h:96
DerivationFramework::ClassifyAndCalculateHFTool::finalize
virtual StatusCode finalize() override
Definition: ClassifyAndCalculateHFTool.cxx:51
SG::Decorator< int >
DerivationFramework::ClassifyAndCalculateHFTool::m_leadingHadronPtRatioCut
Gaudi::Property< float > m_leadingHadronPtRatioCut
Definition: ClassifyAndCalculateHFTool.h:99
beamspotman.n
n
Definition: beamspotman.py:729
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
DerivationFramework
THE reconstruction tool.
Definition: ParticleSortingAlg.h:24
DataVector
Derived DataVector<T>.
Definition: DataVector.h:794
id
SG::auxid_t id
Definition: Control/AthContainers/Root/debug.cxx:239
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
runIDPVM.pdgId
pdgId
Definition: runIDPVM.py:91
DerivationFramework::ClassifyAndCalculateHFTool::m_jetEtaCut
Gaudi::Property< float > m_jetEtaCut
Definition: ClassifyAndCalculateHFTool.h:97
xAOD::Jet_v1
Class describing a jet.
Definition: Jet_v1.h:57
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
DerivationFramework::ClassifyAndCalculateHFTool::getSimpleClassification
int getSimpleClassification(int hfclassif) const
Definition: ClassifyAndCalculateHFTool.cxx:416
isCharmHadron
bool isCharmHadron(const T &p)
Definition: AtlasPID.h:900
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
DerivationFramework::ClassifyAndCalculateHFTool::ClassifyAndCalculateHFTool
ClassifyAndCalculateHFTool(const std::string &t, const std::string &n, const IInterface *p)
Definition: ClassifyAndCalculateHFTool.cxx:25
ClassifyAndCalculateHFTool.h
ConstAccessor.h
Helper class to provide constant type-safe access to aux data.
AthAlgTool
Definition: AthAlgTool.h:26
python.compressB64.c
def c
Definition: compressB64.py:93
HepMCHelpers.h