ATLAS Offline Software
TTbarPlusHeavyFlavorFilter.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "GeneratorFilters/TTbarPlusHeavyFlavorFilter.h"
6 #include "AtlasHepMC/Relatives.h"
7 #include "TruthUtils/MagicNumbers.h"
8 #include "GaudiKernel/MsgStream.h"
9 
10 //--------------------------------------------------------------------------
11 TTbarPlusHeavyFlavorFilter::TTbarPlusHeavyFlavorFilter(const std::string& fname,
12  ISvcLocator* pSvcLocator)
13  : GenFilter(fname,pSvcLocator)
14  {
15 
16  declareProperty("UseFinalStateHadrons",m_useFinalStateHadrons=true);
17  declareProperty("SelectB",m_selectB=false);
18  declareProperty("SelectC",m_selectC=false);
19  declareProperty("SelectL",m_selectL=false);
20  declareProperty("BpTMinCut",m_bPtMinCut=5000.);
21  declareProperty("BetaMaxCut",m_bEtaMaxCut=3.);
22  declareProperty("CpTMinCut",m_cPtMinCut=5000.);
23  declareProperty("CetaMaxCut",m_cEtaMaxCut=3.);
24  declareProperty("BMultiplicityCut",m_bMultiCut=1);
25  declareProperty("CMultiplicityCut",m_cMultiCut=1);
26  declareProperty("ExcludeBFromTTbar",m_excludeBFromTop=true);
27  declareProperty("ExcludeCFromTTbar",m_excludeCFromTop=true);
28 
29 }
30 
31 //--------------------------------------------------------------------------
32  TTbarPlusHeavyFlavorFilter::~TTbarPlusHeavyFlavorFilter(){
34 }
35 
36 //---------------------------------------------------------------------------
37 StatusCode TTbarPlusHeavyFlavorFilter::filterInitialize() {
38  return StatusCode::SUCCESS;
39 }
40 
41 //---------------------------------------------------------------------------
42 StatusCode TTbarPlusHeavyFlavorFilter::filterFinalize() {
43  ATH_MSG_INFO( m_nPass << " Events out of " << m_nPass+m_nFail << " passed the filter" );
44  return StatusCode::SUCCESS;
45 }
46 
47 
48 //---------------------------------------------------------------------------
49 StatusCode TTbarPlusHeavyFlavorFilter::filterEvent() {
50 //---------------------------------------------------------------------------
51 
52  bool pass = false;
53 
54  int nB=0;
55  int nC=0;
56 
57  int nBtop=0;
58  int nCtop=0;
59 
60  McEventCollection::const_iterator itr;
61  for (itr = events()->begin(); itr!=events()->end(); ++itr) {
62 
63  const HepMC::GenEvent* genEvt = (*itr);
64 
65  // Loop over all truth particles in the event
66  // ===========================================
67  for(const auto& part: *genEvt) {
68 
69  if(HepMC::is_simulation_particle(part)) break;
70 
71  bool isbquark=false;
72  bool iscquark=false;
73 
74  bool isbhadron=false;
75  bool ischadron=false;
76 
77  int pdgid = abs(part->pdg_id());
78 
80  if(pdgid == 5 ){
81  isbquark=true;
82  }
83  else if(pdgid == 4 ){
84  iscquark=true;
85  }
86  else if ( isBHadron(part) ){
87  isbhadron=true;
88  }
89  else if ( isCHadron(part) ){
90  ischadron=true;
91  }
92  else{
93  continue;
94  }
95 
96 
97  if( (isbquark || isbhadron) && !passBSelection(part) ) continue;
98  if( (iscquark || ischadron) && !passCSelection(part) ) continue;
99 
100  if(isbhadron || ischadron){
101  if(!isInitialHadron(part) && !m_useFinalStateHadrons) continue;
102  if(!isFinalHadron(part) && m_useFinalStateHadrons) continue;
103  }
104 
105  if(m_excludeBFromTop && isbquark){
106  if(isDirectlyFromTop(part))++nBtop;
107  if(isDirectlyFromWTop(part))++nBtop;
108  }
109  if(m_excludeCFromTop && iscquark){
110  if(isDirectlyFromTop(part))++nCtop;
111  if(isDirectlyFromWTop(part))++nCtop;
112  }
113 
114  bool ischadronfromb = isCHadronFromB(part);
115 
116  if(isbhadron) ++nB;
117  if(ischadron && !ischadronfromb) ++nC;
118 
119 
120  }
121 
122  // Coverity dislikes this break as it means the loop is poitless, commentted out in line with other filters
123  //break;
124 
125  }
126 
127 
128  int nAddB=nB;
129  if(m_excludeBFromTop){
130  nAddB-=nBtop;
131  }
132 
133  int nAddC=nC;
134  if(m_excludeCFromTop){
135  nAddC-=nCtop;
136  }
137 
138  int flavortype=0;
139 
140  if(nAddC>=m_cMultiCut){
141  flavortype=4;
142  }
143 
144  if(nAddB>=m_bMultiCut){
145  flavortype=5;
146  }
147 
148  if(m_selectB && 5 == flavortype) pass=true;
149  if(m_selectC && 4 == flavortype) pass=true;
150  if(m_selectL && 0 == flavortype) pass=true;
151 
152  setFilterPassed(pass);
153  return StatusCode::SUCCESS;
154 
155 }
156 
157 
158 bool TTbarPlusHeavyFlavorFilter::passBSelection(const HepMC::ConstGenParticlePtr& part) const{
159 
160  const HepMC::FourVector& p4 = part->momentum();
161  double pt = p4.perp();
162  double eta = std::abs(p4.eta());
163 
164  if(pt<m_bPtMinCut) return false;
165  if(eta>m_bEtaMaxCut) return false;
166 
167  return true;
168 
169 }
170 
171 bool TTbarPlusHeavyFlavorFilter::passCSelection(const HepMC::ConstGenParticlePtr& part) const{
172 
173  const HepMC::FourVector& p4 = part->momentum();
174  double pt = p4.perp();
175  double eta = std::abs(p4.eta());
176 
177  if(pt<m_cPtMinCut) return false;
178  if(eta>m_cEtaMaxCut) return false;
179 
180  return true;
181 }
182 
183 
184 int TTbarPlusHeavyFlavorFilter::hadronType(int pdgid) const{
185 
186  int rest1(std::abs(pdgid%1000));
187  int rest2(std::abs(pdgid%10000));
188 
189  if ( rest2 >= 5000 && rest2 < 6000 ) return 5;
190  if( rest1 >= 500 && rest1 < 600 ) return 5;
191 
192  if ( rest2 >= 4000 && rest2 < 5000 ) return 4;
193  if( rest1 >= 400 && rest1 < 500 ) return 4;
194 
195  return 0;
196 
197 }
198 
199 
200 bool TTbarPlusHeavyFlavorFilter::isBHadron(const HepMC::ConstGenParticlePtr& part) const{
201 
202  if(HepMC::is_simulation_particle(part)) return false;
203  int type = hadronType(part->pdg_id());
204  if(type == 5) return true;
205 
206  return false;
207 
208 }
209 
210 
211 bool TTbarPlusHeavyFlavorFilter::isCHadron(const HepMC::ConstGenParticlePtr& part) const{
212 
213  if(HepMC::is_simulation_particle(part)) return false;
214  int type = hadronType(part->pdg_id());
215  if(type == 4) return true;
216 
217  return false;
218 
219 }
220 
221 
222 
223 bool TTbarPlusHeavyFlavorFilter::isInitialHadron(const HepMC::ConstGenParticlePtr& part) const{
224 
225  auto prod = part->production_vertex();
226  if(!prod) return true;
227  int type = hadronType(part->pdg_id());
228 #ifdef HEPMC3
229  for(const auto& firstParent: prod->particles_in()){
230  int mothertype = hadronType( firstParent->pdg_id() );
231  if( mothertype == type ){
232  return false;
233  }
234  }
235 #else
236  HepMC::GenVertex::particle_iterator firstParent = prod->particles_begin(HepMC::parents);
237  HepMC::GenVertex::particle_iterator endParent = prod->particles_end(HepMC::parents);
238  for(;firstParent!=endParent; ++firstParent){
239  int mothertype = hadronType( (*firstParent)->pdg_id() );
240  if( mothertype == type ){
241  return false;
242  }
243  }
244 #endif
245 
246  return true;
247 }
248 
249 
250 bool TTbarPlusHeavyFlavorFilter::isFinalHadron(const HepMC::ConstGenParticlePtr& part) const{
251 
252  auto end = part->end_vertex();
253  if(!end) return true;
254  int type = hadronType(part->pdg_id());
255 #ifdef HEPMC3
256  for(const auto& firstChild: end->particles_in()){
257  int childtype = hadronType( firstChild->pdg_id() );
258  if( childtype == type ){
259  return false;
260  }
261  }
262 #else
263  HepMC::GenVertex::particle_iterator firstChild = end->particles_begin(HepMC::children);
264  HepMC::GenVertex::particle_iterator endChild = end->particles_end(HepMC::children);
265  for(;firstChild!=endChild; ++firstChild){
266  int childtype = hadronType( (*firstChild)->pdg_id() );
267  if( childtype == type ){
268  return false;
269  }
270  }
271 #endif
272  return true;
273 
274 }
275 
276 
277 
278 bool TTbarPlusHeavyFlavorFilter::isQuarkFromHadron(const HepMC::ConstGenParticlePtr& part) const{
279 
280  auto prod = part->production_vertex();
281  if(!prod) return false;
282 #ifdef HEPMC3
283  for(const auto& firstParent: HepMC::ancestor_particles(prod)){
284  int mothertype = hadronType( firstParent->pdg_id() );
285  if( 4 == mothertype || 5 == mothertype ){
286  return true;
287  }
288  }
289 #else
290  HepMC::GenVertex::particle_iterator firstParent = prod->particles_begin(HepMC::ancestors);
291  HepMC::GenVertex::particle_iterator endParent = prod->particles_end(HepMC::ancestors);
292  for(;firstParent!=endParent; ++firstParent){
293  int mothertype = hadronType( (*firstParent)->pdg_id() );
294  if( 4 == mothertype || 5 == mothertype ){
295  return true;
296  }
297  }
298 #endif
299  return false;
300 
301 }
302 
303 bool TTbarPlusHeavyFlavorFilter::isCHadronFromB(const HepMC::ConstGenParticlePtr& part) const{
304 
305  if(!isCHadron(part)) return false;
306 
307  auto prod = part->production_vertex();
308  if(!prod) return false;
309 #ifdef HEPMC3
310  for(const auto& firstParent:HepMC::ancestor_particles(prod)){
311  if( isBHadron(firstParent) ){
312  return true;
313  }
314  }
315 #else
316  HepMC::GenVertex::particle_iterator firstParent = prod->particles_begin(HepMC::ancestors);
317  HepMC::GenVertex::particle_iterator endParent = prod->particles_end(HepMC::ancestors);
318  for(;firstParent!=endParent; ++firstParent){
319  if( isBHadron(*firstParent) ){
320  return true;
321  }
322  }
323 #endif
324  return false;
325 }
326 
327 
328 
329 
330 
331 HepMC::ConstGenParticlePtr TTbarPlusHeavyFlavorFilter::findInitial(const HepMC::ConstGenParticlePtr& part) const{
332 
333  auto prod = part->production_vertex();
334 
335  if(!prod) return part;
336 #ifdef HEPMC3
337  for(const auto& firstParent: prod->particles_in()){
338  if( part->pdg_id() == firstParent->pdg_id() ){
339  return findInitial(firstParent);
340  }
341  }
342 #else
343  HepMC::GenVertex::particle_iterator firstParent = prod->particles_begin(HepMC::parents);
344  HepMC::GenVertex::particle_iterator endParent = prod->particles_end(HepMC::parents);
345  for(;firstParent!=endParent; ++firstParent){
346  if( part->pdg_id() == (*firstParent)->pdg_id() ){
347  return findInitial(*firstParent);
348  }
349  }
350 #endif
351 
352  return part;
353 
354 }
355 
356 bool TTbarPlusHeavyFlavorFilter::isFromTop(const HepMC::ConstGenParticlePtr& part) const{
357 
358  auto initpart = findInitial(part);
359  return isDirectlyFromTop(initpart);
360 
361 }
362 
363 bool TTbarPlusHeavyFlavorFilter::isDirectlyFromTop(const HepMC::ConstGenParticlePtr& part) const{
364 
365  auto prod = part->production_vertex();
366 
367  if(!prod) return false;
368 #ifdef HEPMC3
369  for (auto firstParent: prod->particles_in()){
370  if( std::abs( firstParent->pdg_id() ) == 6 ) return true;
371  }
372 #else
373  HepMC::GenVertex::particle_iterator firstParent = prod->particles_begin(HepMC::parents);
374  HepMC::GenVertex::particle_iterator endParent = prod->particles_end(HepMC::parents);
375  for(;firstParent!=endParent; ++firstParent){
376  if( std::abs( (*firstParent)->pdg_id() ) == 6 ) return true;
377  }
378 #endif
379 
380  return false;
381 }
382 
383 
384 
385 bool TTbarPlusHeavyFlavorFilter::isDirectlyFromWTop(const HepMC::ConstGenParticlePtr& part) const{
386 
387  auto prod = part->production_vertex();
388 
389  if(!prod) return false;
390 #ifdef HEPMC3
391  for(const auto& firstParent: prod->particles_in()){
392  if( std::abs( firstParent->pdg_id() ) == 24 ){
393  if( isFromTop(firstParent) ) return true;
394  }
395  }
396 #else
397  HepMC::GenVertex::particle_iterator firstParent = prod->particles_begin(HepMC::parents);
398  HepMC::GenVertex::particle_iterator endParent = prod->particles_end(HepMC::parents);
399  for(;firstParent!=endParent; ++firstParent){
400  if( std::abs( (*firstParent)->pdg_id() ) == 24 ){
401  if( isFromTop(*firstParent) ) return true;
402  }
403  }
404 #endif
405 
406  return false;
407 
408 
409 }
410 
411 
LArG4FSStartPointFilter.part
part
Definition: LArG4FSStartPointFilter.py:21
TTbarPlusHeavyFlavorFilter::hadronType
int hadronType(int pdgid) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:184
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
TTbarPlusHeavyFlavorFilter::isBHadron
bool isBHadron(const HepMC::ConstGenParticlePtr &part) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:200
TTbarPlusHeavyFlavorFilter::TTbarPlusHeavyFlavorFilter
TTbarPlusHeavyFlavorFilter(const std::string &fname, ISvcLocator *pSvcLocator)
Definition: TTbarPlusHeavyFlavorFilter.cxx:11
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
TTbarPlusHeavyFlavorFilter.h
filter to select ttbar+jets/HF events
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:79
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
python.DecayParser.parents
parents
print ("==> buf:",buf)
Definition: DecayParser.py:31
TTbarPlusHeavyFlavorFilter::m_excludeBFromTop
bool m_excludeBFromTop
Definition: TTbarPlusHeavyFlavorFilter.h:46
PlotCalibFromCool.begin
begin
Definition: PlotCalibFromCool.py:94
test_pyathena.pt
pt
Definition: test_pyathena.py:11
TTbarPlusHeavyFlavorFilter::isDirectlyFromTop
bool isDirectlyFromTop(const HepMC::ConstGenParticlePtr &part) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:363
TTbarPlusHeavyFlavorFilter::isCHadron
bool isCHadron(const HepMC::ConstGenParticlePtr &part) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:211
python.DataFormatRates.events
events
Definition: DataFormatRates.py:105
mergePhysValFiles.end
end
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:93
TTbarPlusHeavyFlavorFilter::filterFinalize
virtual StatusCode filterFinalize()
Definition: TTbarPlusHeavyFlavorFilter.cxx:42
TTbarPlusHeavyFlavorFilter::passCSelection
bool passCSelection(const HepMC::ConstGenParticlePtr &part) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:171
GenFilter
Base class for event generator filtering modules.
Definition: GenFilter.h:30
HepMC::is_simulation_particle
bool is_simulation_particle(const T &p)
Method to establish if a particle (or barcode) was created during the simulation (TODO update to be s...
Definition: MagicNumbers.h:299
TTbarPlusHeavyFlavorFilter::isQuarkFromHadron
bool isQuarkFromHadron(const HepMC::ConstGenParticlePtr &part) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:278
TTbarPlusHeavyFlavorFilter::passBSelection
bool passBSelection(const HepMC::ConstGenParticlePtr &part) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:158
TTbarPlusHeavyFlavorFilter::~TTbarPlusHeavyFlavorFilter
virtual ~TTbarPlusHeavyFlavorFilter()
Definition: TTbarPlusHeavyFlavorFilter.cxx:32
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
TTbarPlusHeavyFlavorFilter::m_excludeCFromTop
bool m_excludeCFromTop
Definition: TTbarPlusHeavyFlavorFilter.h:47
TTbarPlusHeavyFlavorFilter::m_useFinalStateHadrons
bool m_useFinalStateHadrons
properties
Definition: TTbarPlusHeavyFlavorFilter.h:32
TTbarPlusHeavyFlavorFilter::m_bPtMinCut
double m_bPtMinCut
Definition: TTbarPlusHeavyFlavorFilter.h:38
TTbarPlusHeavyFlavorFilter::filterEvent
virtual StatusCode filterEvent()
Definition: TTbarPlusHeavyFlavorFilter.cxx:49
TTbarPlusHeavyFlavorFilter::m_selectL
bool m_selectL
Definition: TTbarPlusHeavyFlavorFilter.h:35
TTbarPlusHeavyFlavorFilter::m_cMultiCut
int m_cMultiCut
Definition: TTbarPlusHeavyFlavorFilter.h:44
TTbarPlusHeavyFlavorFilter::m_selectB
bool m_selectB
Definition: TTbarPlusHeavyFlavorFilter.h:33
TTbarPlusHeavyFlavorFilter::m_bMultiCut
int m_bMultiCut
Definition: TTbarPlusHeavyFlavorFilter.h:43
TTbarPlusHeavyFlavorFilter::isInitialHadron
bool isInitialHadron(const HepMC::ConstGenParticlePtr &part) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:223
HepMC::ConstGenParticlePtr
const GenParticle * ConstGenParticlePtr
Definition: GenParticle.h:38
TTbarPlusHeavyFlavorFilter::m_cEtaMaxCut
double m_cEtaMaxCut
Definition: TTbarPlusHeavyFlavorFilter.h:41
TTbarPlusHeavyFlavorFilter::findInitial
HepMC::ConstGenParticlePtr findInitial(const HepMC::ConstGenParticlePtr &part) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:331
MagicNumbers.h
GenFilter::m_nPass
int m_nPass
Definition: GenFilter.h:65
TTbarPlusHeavyFlavorFilter::isFinalHadron
bool isFinalHadron(const HepMC::ConstGenParticlePtr &part) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:250
TTbarPlusHeavyFlavorFilter::m_cPtMinCut
double m_cPtMinCut
Definition: TTbarPlusHeavyFlavorFilter.h:40
TTbarPlusHeavyFlavorFilter::m_selectC
bool m_selectC
Definition: TTbarPlusHeavyFlavorFilter.h:34
python.AthDsoLogger.fname
string fname
Definition: AthDsoLogger.py:67
TTbarPlusHeavyFlavorFilter::isCHadronFromB
bool isCHadronFromB(const HepMC::ConstGenParticlePtr &part) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:303
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
Relatives.h
python.DecayParser.children
children
Definition: DecayParser.py:32
GenFilter::m_nFail
int m_nFail
Definition: GenFilter.h:66
TTbarPlusHeavyFlavorFilter::isFromTop
bool isFromTop(const HepMC::ConstGenParticlePtr &part) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:356
TTbarPlusHeavyFlavorFilter::filterInitialize
virtual StatusCode filterInitialize()
Definition: TTbarPlusHeavyFlavorFilter.cxx:37
TTbarPlusHeavyFlavorFilter::m_bEtaMaxCut
double m_bEtaMaxCut
Definition: TTbarPlusHeavyFlavorFilter.h:39
TTbarPlusHeavyFlavorFilter::isDirectlyFromWTop
bool isDirectlyFromWTop(const HepMC::ConstGenParticlePtr &part) const
Definition: TTbarPlusHeavyFlavorFilter.cxx:385
Generated on Fri Jul 5 2024 21:32:13 for ATLAS Offline Software by   doxygen 1.8.18