CMS 3D CMS Logo

ttHFGenFilter.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: ttHFGenFilter/ttHFGenFilter
4 // Class: ttHFGenFilter
5 //
18 //
19 // Original Author: Andrej Saibel
20 // Created: Tue, 05 Jul 2016 09:36:09 GMT
21 //
22 // VERY IMPORTANT: when running this code, you should make sure, that the GenHFHadronMatcher runs in onlyJetClusteredHadrons = cms.bool(False) mode
23 // If you want to filter all events with additional b-hadrons use: OnlyHardProcessBHadrons = cms.bool(False)
24 
25 
26 // system include files
27 #include <memory>
28 
29 // user include files
32 
36 
40 
41 
42 //
43 // class declaration
44 //
45 
47  public:
48  explicit ttHFGenFilter(const edm::ParameterSet&);
49  ~ttHFGenFilter() override;
50 
51  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
52 
53  private:
54  void beginStream(edm::StreamID) override;
55  bool filter(edm::Event&, const edm::EventSetup&) override;
56  void endStream() override;
57 
58  virtual bool HasAdditionalBHadron(const std::vector<int>&, const std::vector<int>&,const std::vector<reco::GenParticle>&,std::vector<const reco::Candidate*>&);
59  virtual bool analyzeMothersRecursive(const reco::Candidate*,std::vector<const reco::Candidate*>& AllTopMothers);
60  virtual std::vector<const reco::Candidate*> GetTops(const std::vector<reco::GenParticle> &,std::vector<const reco::Candidate*>& AllTopMothers);
61  virtual void FindAllTopMothers(const reco::Candidate* particle,std::vector<const reco::Candidate*>&);
62 
71 
72 };
73 
74 //
75 // constants, enums and typedefs
76 //
77 
78 //
79 // static data member definitions
80 //
81 
82 //
83 // constructors and destructor
84 //
86 genParticlesToken_(consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genParticles"))),
87 genBHadFlavourToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadFlavour"))),
88 genBHadFromTopWeakDecayToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadFromTopWeakDecay"))),
89 genBHadPlusMothersToken_(consumes<std::vector<reco::GenParticle> >(iConfig.getParameter<edm::InputTag>("genBHadPlusMothers"))),
90 genBHadPlusMothersIndicesToken_(consumes<std::vector<std::vector<int> > >(iConfig.getParameter<edm::InputTag>("genBHadPlusMothersIndices"))),
91 genBHadIndexToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadIndex")))
92 {
93  //now do what ever initialization is needed
94  OnlyHardProcessBHadrons_ = iConfig.getParameter<bool> ( "OnlyHardProcessBHadrons" );
95  taggingMode_ = iConfig.getParameter<bool>("taggingMode");
96 
97  produces<bool>();
98 }
99 
100 
102 {
103 
104  // do anything here that needs to be done at destruction time
105  // (e.g. close files, deallocate resources etc.)
106 
107 }
108 
109 
110 //
111 // member functions
112 //
113 
114 // ------------ method called on each new Event ------------
115 bool
117 {
118  using namespace std;
119 
120  //get GenParticleCollection
122  iEvent.getByToken(genParticlesToken_, genParticles);
123 
124  //get Information on B-Hadrons
125  // the information is calculated in GenHFHadronMatcher
126 
128  iEvent.getByToken(genBHadFlavourToken_, genBHadFlavour);
129 
131  iEvent.getByToken(genBHadFromTopWeakDecayToken_, genBHadFromTopWeakDecay);
132 
134  iEvent.getByToken(genBHadPlusMothersToken_, genBHadPlusMothers);
135 
137  iEvent.getByToken(genBHadPlusMothersIndicesToken_, genBHadPlusMothersIndices);
138 
140  iEvent.getByToken(genBHadIndexToken_, genBHadIndex);
141 
142 
143  //check whether the event has additional b-hadrons not coming from top/topbar decay
144 
145  std::vector<const reco::Candidate*> AllTopMothers;
146  std::vector<const reco::Candidate*> Tops = GetTops(*genParticles,AllTopMothers);
147 
148  bool pass = HasAdditionalBHadron(*genBHadIndex,*genBHadFlavour,*genBHadPlusMothers,AllTopMothers);
149 
150  iEvent.put(std::make_unique<bool>(pass));
151 
152  return taggingMode_ || pass;
153 
154 }
155 
156 bool ttHFGenFilter::HasAdditionalBHadron(const std::vector<int>& genBHadIndex, const std::vector<int>& genBHadFlavour,const std::vector<reco::GenParticle>& genBHadPlusMothers,std::vector<const reco::Candidate*>& AllTopMothers){
157 
158  for(unsigned int i=0; i<genBHadIndex.size();i++){
159 
160  const reco::GenParticle* bhadron = genBHadIndex[i]>=0&&genBHadIndex[i]<int(genBHadPlusMothers.size()) ? &(genBHadPlusMothers[genBHadIndex[i]]) : nullptr;
161  int motherflav = genBHadFlavour[i];
162  bool from_tth=(abs(motherflav)==6||abs(motherflav)==25); //b-hadron comes from top or higgs decay
163  bool fromhp=false;
164 
166  if(bhadron!=nullptr&&!from_tth){
167  return true;
168  }
169  }
170 
172  if(bhadron!=nullptr&&!from_tth){
173  //std::cout << "PT: " << bhadron->pt() << " , eta: " << bhadron->eta() << std::endl;
174 
175  fromhp=analyzeMothersRecursive(bhadron,AllTopMothers);
176  if(fromhp){
177  return true;
178  }
179  }
180  if(i==genBHadIndex.size()-1){
181  return false;
182  }
183  }
184 }
185 
186  return false;
187 }
188 
189 //recursive function, that loops over all chain particles until it finds the protons
190 bool ttHFGenFilter::analyzeMothersRecursive(const reco::Candidate* particle,std::vector<const reco::Candidate*>& AllTopMothers){
191  if(particle->status()>20&&particle->status()<30){ //particle comes from hardest process in event
192  return true;
193  }
194  for(unsigned int k=0; k<AllTopMothers.size();k++){
195  if(particle==AllTopMothers[k]){ //particle comes from ISR
196  return true;
197  }
198  }
199  bool isFromHardProcess=false;
200  for(unsigned int i=0;i<particle->numberOfMothers();i++){
201  const reco::Candidate* mother = particle->mother(i);
202  isFromHardProcess=analyzeMothersRecursive(mother,AllTopMothers);
203  if(isFromHardProcess){
204  return true;
205  }
206  }
207  return isFromHardProcess;
208 }
209 
210 std::vector< const reco::Candidate*> ttHFGenFilter::GetTops(const std::vector<reco::GenParticle>& genParticles, std::vector<const reco::Candidate*>& AllTopMothers){
211  // std::vector<const reco::Candidate*> tops;
212  const reco::GenParticle* firstTop = nullptr;
213  const reco::GenParticle* firstTopBar = nullptr;
214  bool foundTop = false;
215  bool foundTopBar = false;
216  std::vector<const reco::GenParticle*> Tops;
217 
218  //loop over all genParticles and find a Top and AntiTop quark
219  //then find all mothers of the Tops
220  //this is then used to if the b-hadron is an inital state radiation particle
221 
222  for(reco::GenParticleCollection::const_iterator i_particle = genParticles.begin(); i_particle != genParticles.end(); ++i_particle){
223  const reco::GenParticle* thisParticle = &*i_particle;
224  if(!foundTop && thisParticle->pdgId()==6){
225  firstTop = thisParticle;
226  for(unsigned int i=0; i<firstTop->numberOfMothers();i++){
227  FindAllTopMothers(thisParticle->mother(i),AllTopMothers);
228  }
229  foundTop = true;
230  }
231  else if(!foundTopBar && thisParticle->pdgId()==-6){
232  firstTopBar = thisParticle;
233  for(unsigned int i=0; i<firstTopBar->numberOfMothers();i++){
234  FindAllTopMothers(firstTopBar->mother(i),AllTopMothers);
235  }
236  foundTopBar = true;
237  }
238  if(foundTopBar&&foundTop){
239  //tops.push_back(firstTop);
240  //tops.push_back(firstTopBar);
241  //return tops;
242  break;
243  }
244  }
245  return AllTopMothers;
246 }
247 
248 // finds all mothers of "particle", that are not tops or protons
249 
250 void ttHFGenFilter::FindAllTopMothers(const reco::Candidate* particle, std::vector<const reco::Candidate*>& AllTopMothers){
251  // std::cout << "particle mother: " << particle->mother(0) << std::endl;
252  for(unsigned int i=0;i<particle->numberOfMothers();i++){
253  if(abs(particle->mother(i)->pdgId())!=6 && particle->mother(i)->pdgId()!=2212){
254  AllTopMothers.push_back(particle->mother(i));
255  if(particle->mother(i)->pdgId()!=2212 || particle->mother(i)->numberOfMothers()>1){
256  //std::cout << "Size of vector in loop = " << AllTopMothers.size() << std::endl;
257  FindAllTopMothers(particle->mother(i),AllTopMothers);
258  }
259  }
260  }
261 
262 }
263 
264 
265 // ------------ method called once each stream before processing any runs, lumis or events ------------
266 void
268 {
269 }
270 
271 // ------------ method called once each stream after processing all runs, lumis and events ------------
272 void
274 }
275 
276 // ------------ method called when starting to processes a run ------------
277 /*
278 void
279 ttHFGenFilter::beginRun(edm::Run const&, edm::EventSetup const&)
280 {
281 }
282 */
283 
284 // ------------ method called when ending the processing of a run ------------
285 /*
286 void
287 ttHFGenFilter::endRun(edm::Run const&, edm::EventSetup const&)
288 {
289 }
290 */
291 
292 // ------------ method called when starting to processes a luminosity block ------------
293 /*
294 void
295 ttHFGenFilter::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
296 {
297 }
298 */
299 
300 // ------------ method called when ending the processing of a luminosity block ------------
301 /*
302 void
303 ttHFGenFilter::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
304 {
305 }
306 */
307 
308 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
309 void
311  //The following says we do not know what parameters are allowed so do no validation
312  // Please change this to state exactly what you do use, even if it is no parameters
314  desc.setUnknown();
315  descriptions.addDefault(desc);
316 }
317 //define this as a plug-in
const edm::EDGetTokenT< std::vector< reco::GenParticle > > genBHadPlusMothersToken_
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
T getParameter(std::string const &) const
int pdgId() const final
PDG identifier.
virtual void FindAllTopMothers(const reco::Candidate *particle, std::vector< const reco::Candidate * > &)
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:125
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:517
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
ttHFGenFilter(const edm::ParameterSet &)
size_t numberOfMothers() const override
number of mothers
virtual bool analyzeMothersRecursive(const reco::Candidate *, std::vector< const reco::Candidate * > &AllTopMothers)
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
void beginStream(edm::StreamID) override
void endStream() override
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
const edm::EDGetTokenT< std::vector< int > > genBHadFlavourToken_
~ttHFGenFilter() override
virtual int status() const =0
status word
int iEvent
Definition: GenABIO.cc:224
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
void addDefault(ParameterSetDescription const &psetDescription)
const edm::EDGetTokenT< reco::GenParticleCollection > genParticlesToken_
const edm::EDGetTokenT< std::vector< std::vector< int > > > genBHadPlusMothersIndicesToken_
virtual int pdgId() const =0
PDG identifier.
const edm::EDGetTokenT< std::vector< int > > genBHadFromTopWeakDecayToken_
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int k[5][pyjets_maxn]
const edm::EDGetTokenT< std::vector< int > > genBHadIndexToken_
virtual bool HasAdditionalBHadron(const std::vector< int > &, const std::vector< int > &, const std::vector< reco::GenParticle > &, std::vector< const reco::Candidate * > &)
fixed size matrix
HLT enums.
bool filter(edm::Event &, const edm::EventSetup &) override
virtual std::vector< const reco::Candidate * > GetTops(const std::vector< reco::GenParticle > &, std::vector< const reco::Candidate * > &AllTopMothers)
const Candidate * mother(size_type=0) const override
return mother at a given position, i = 0, ... numberOfMothers() - 1 (read only mode) ...
bool OnlyHardProcessBHadrons_