CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
JetFlavourIdentifier.cc
Go to the documentation of this file.
1 //
2 // Translation of BTag MCJetFlavour tool to identify real flavour of a jet
3 // work with CaloJet objects
4 // Store Infos by Values in JetFlavour.h
5 // Author: Attilio
6 // Date: 05.10.2007
7 //
8 //
9 // \class JetFlavourIdentifier
10 //
11 // \brief Interface to pull out the proper flavour identifier from a
12 // jet->parton matching collection.
13 //
14 // In detail, the definitions are as follows:
15 //
16 // Definitions:
17 // The default behavior is that the "definition" is NULL_DEF,
18 // so the software will fall back to either "physics" or "algorithmic" definition
19 // as per the "physDefinition" switch.
20 //
21 // If the user specifies "definition", then that definition is taken. However,
22 // if the requested definition is not defined, the software reverts back to
23 // either "physics" or "algorithmic" based on the "physDefinition" switch.
24 // For example, if the user specifies "heaviest" as a flavor ID, and there
25 // are no bottom, charm, or top quarks in the event that match to the jet,
26 // then the software will fall back to the "physics" or "algorithmic" definition.
27 //
28 // Modifications:
29 //
30 // 09.03.2008: Sal Rappoccio.
31 // Added capability for all methods of identification, not just
32 // "physics" or "algorithmic". If the requested method does not exist
33 // (i.e. is unphysical), the "physics" or "algorithmic" definition
34 // is defaulted to.
35 
36 //=======================================================================
37 
38 // user include files
43 
50 
51 //#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
52 
53 
61 
65 
68 
69 #include <memory>
70 #include <string>
71 #include <iostream>
72 #include <vector>
73 #include <Math/VectorUtil.h>
74 #include <TMath.h>
75 
76 using namespace std;
77 using namespace reco;
78 using namespace edm;
79 using namespace ROOT::Math::VectorUtil;
80 
81 namespace reco { namespace modules {
82 
83 //--------------------------------------------------------------------------
84 //
85 //--------------------------------------------------------------------------
87 {
88  public:
89  enum DEFINITION_T { PHYSICS=0, ALGO, NEAREST_STATUS2, NEAREST_STATUS3, HEAVIEST,
91  NULL_DEF};
92 
95 
96  private:
97  virtual void produce(edm::Event&, const edm::EventSetup& ) override;
98 
99  JetFlavour::Leptons findLeptons(const GenParticleRef &);
100  std::vector<const reco::Candidate*> findCandidates(const reco::Candidate*, int);
101  void fillLeptons(const std::vector<const reco::Candidate*> &, JetFlavour::Leptons &, int, int);
102  static int heaviestFlavour(int);
103 
110 
111 };
112 } }
114 
115 //=========================================================================
116 
117 JetFlavourIdentifier::JetFlavourIdentifier( const edm::ParameterSet& iConfig )
118 {
119  produces<JetFlavourMatchingCollection>();
120  sourceByReferToken_ = consumes<JetMatchedPartonsCollection>(iConfig.getParameter<InputTag>("srcByReference"));
121  physDefinition = iConfig.getParameter<bool>("physicsDefinition");
122  leptonInfo_ = iConfig.exists("leptonInfo") ? iConfig.getParameter<bool>("leptonInfo") : false;
123  // If we have a definition of which parton to identify, use it,
124  // otherwise we default to the "old" behavior of either "physics" or "algorithmic".
125  // Furthermore, if the specified definition is not sensible for the given jet,
126  // then the "physDefinition" switch is used to identify the flavour of the jet.
127  if ( iConfig.exists("definition") ) {
128  definition = static_cast<DEFINITION_T>( iConfig.getParameter<int>("definition") );
129  } else {
130  definition = NULL_DEF;
131  }
132 }
133 
134 //=========================================================================
135 
136 JetFlavourIdentifier::~JetFlavourIdentifier()
137 {
138 }
139 
140 // ------------ method called to produce the data ------------
141 
142 void JetFlavourIdentifier::produce( Event& iEvent, const EventSetup& iEs )
143 {
144  // Get the JetMatchedPartons
145  iEvent.getByToken (sourceByReferToken_, theTagByRef);
146 
147  // Create a JetFlavourMatchingCollection
149  if (!theTagByRef->empty()) {
150  RefToBase<Jet> jj = theTagByRef->begin()->first;
152  } else {
153  jfmc = new JetFlavourMatchingCollection();
154  }
155  auto_ptr<reco::JetFlavourMatchingCollection> jetFlavMatching(jfmc);
156 
157  // Loop over the matched partons and see which match.
158  for ( JetMatchedPartonsCollection::const_iterator j = theTagByRef->begin();
159  j != theTagByRef->end();
160  j ++ ) {
161 
162 
163  // Consider this match.
164  const MatchedPartons aMatch = (*j).second;
165 
166  // This will hold the 4-vector, vertex, flavour and the leptonian decays (0: no lepton, xyz: x leptons in layer 2, y in layer 1 and z in layer 0) of the requested object.
167  math::XYZTLorentzVector thePartonLorentzVector(0,0,0,0);
168  math::XYZPoint thePartonVertex(0,0,0);
169  int thePartonFlavour = 0;
170  JetFlavour::Leptons theLeptons;
171 
172  // get the partons based on which definition to use.
173  switch (definition) {
174  case PHYSICS: {
175  const GenParticleRef aPartPhy = aMatch.physicsDefinitionParton();
176  if (aPartPhy.isNonnull()) {
177  thePartonLorentzVector = aPartPhy.get()->p4();
178  thePartonVertex = aPartPhy.get()->vertex();
179  thePartonFlavour = aPartPhy.get()->pdgId();
180  if (leptonInfo_) theLeptons = findLeptons(aPartPhy);
181  }
182  break;
183  }
184  case ALGO: {
185  const GenParticleRef aPartAlg = aMatch.algoDefinitionParton();
186  if (aPartAlg.isNonnull()) {
187  thePartonLorentzVector = aPartAlg.get()->p4();
188  thePartonVertex = aPartAlg.get()->vertex();
189  thePartonFlavour = aPartAlg.get()->pdgId();
190  if (leptonInfo_) theLeptons = findLeptons(aPartAlg);
191  }
192  break;
193  }
194  case NEAREST_STATUS2 : {
195  const GenParticleRef aPartN2 = aMatch.nearest_status2();
196  if (aPartN2.isNonnull()) {
197  thePartonLorentzVector = aPartN2.get()->p4();
198  thePartonVertex = aPartN2.get()->vertex();
199  thePartonFlavour = aPartN2.get()->pdgId();
200  if (leptonInfo_) theLeptons = findLeptons(aPartN2);
201  }
202  break;
203  }
204  case NEAREST_STATUS3: {
205  const GenParticleRef aPartN3 = aMatch.nearest_status3();
206  if (aPartN3.isNonnull()) {
207  thePartonLorentzVector = aPartN3.get()->p4();
208  thePartonVertex = aPartN3.get()->vertex();
209  thePartonFlavour = aPartN3.get()->pdgId();
210  if (leptonInfo_) theLeptons = findLeptons(aPartN3);
211  }
212  break;
213  }
214  case HEAVIEST: {
215  const GenParticleRef aPartHeaviest = aMatch.heaviest();
216  if (aPartHeaviest.isNonnull()) {
217  thePartonLorentzVector = aPartHeaviest.get()->p4();
218  thePartonVertex = aPartHeaviest.get()->vertex();
219  thePartonFlavour = aPartHeaviest.get()->pdgId();
220  if (leptonInfo_) theLeptons = findLeptons(aPartHeaviest);
221  }
222  break;
223  }
224  // Default case is backwards-compatible
225  default:{
226  if (physDefinition) {
227  const GenParticleRef aPartPhy = aMatch.physicsDefinitionParton();
228  if (aPartPhy.isNonnull()) {
229  thePartonLorentzVector = aPartPhy.get()->p4();
230  thePartonVertex = aPartPhy.get()->vertex();
231  thePartonFlavour = aPartPhy.get()->pdgId();
232  if (leptonInfo_) theLeptons = findLeptons(aPartPhy);
233  }
234  } else {
235  const GenParticleRef aPartAlg = aMatch.algoDefinitionParton();
236  if (aPartAlg.isNonnull()) {
237  thePartonLorentzVector = aPartAlg.get()->p4();
238  thePartonVertex = aPartAlg.get()->vertex();
239  thePartonFlavour = aPartAlg.get()->pdgId();
240  if (leptonInfo_) theLeptons = findLeptons(aPartAlg);
241  }
242  }
243  } break;
244  } // end switch on definition
245 
246  // Now make sure we have a match. If the user specified "heaviest", for instance,
247  // and there is no b- or c-quarks, then fall back to the "physDefinition" switch.
248 
249  if (thePartonFlavour == 0) {
250  if (physDefinition) {
251  const GenParticleRef aPartPhy = aMatch.physicsDefinitionParton();
252  if (aPartPhy.isNonnull()) {
253  thePartonLorentzVector = aPartPhy.get()->p4();
254  thePartonVertex = aPartPhy.get()->vertex();
255  thePartonFlavour = aPartPhy.get()->pdgId();
256  if (leptonInfo_) theLeptons = findLeptons(aPartPhy);
257  }
258  } else {
259  const GenParticleRef aPartAlg = aMatch.algoDefinitionParton();
260  if (aPartAlg.isNonnull()) {
261  thePartonLorentzVector = aPartAlg.get()->p4();
262  thePartonVertex = aPartAlg.get()->vertex();
263  thePartonFlavour = aPartAlg.get()->pdgId();
264  if (leptonInfo_) theLeptons = findLeptons(aPartAlg);
265  }
266  }
267  }
268 
269 /*
270  std::cout << "Leptons of " <<thePartonFlavour << " Jet: " << std::endl;
271  std::cout << " electrons: " <<theLeptons.electron << std::endl;
272  std::cout << " muons : " <<theLeptons.muon << std::endl;
273  std::cout << " tau : " <<theLeptons.tau << std::endl;
274 */
275  // Add the jet->flavour match to the map.
276  (*jetFlavMatching)[(*j).first] = JetFlavour(thePartonLorentzVector, thePartonVertex, thePartonFlavour, theLeptons);
277  }// end loop over jets
278 
279 
280  // Put the object into the event.
281  iEvent.put( jetFlavMatching );
282 
283 }
284 
285 JetFlavour::Leptons JetFlavourIdentifier::findLeptons(const GenParticleRef &parton)
286 {
287  JetFlavour::Leptons theLeptons;
288 
289  thePartonLV = parton->p4();
290 
292  const reco::Candidate *mcstring = parton->daughter(0);
293  int partonFlavour = std::abs(parton->pdgId());
294 // std::cout << "parton DeltaR: " << DeltaR(thePartonLV, parton->p4()) << std::endl;
295 
297  std::vector<const reco::Candidate*> candidates = findCandidates(mcstring, partonFlavour);
298 // std::cout << "Candidates are:" << std::endl;
299 // for(unsigned int j = 0; j < candidates.size(); j++) std::cout << " --> " << candidates[j]->pdgId() << std::endl;
300 
302  fillLeptons(candidates, theLeptons, 1, partonFlavour);
303 
304  return theLeptons;
305 }
306 
307 std::vector<const reco::Candidate*> JetFlavourIdentifier::findCandidates(const reco::Candidate *cand, int partonFlavour)
308 {
309  std::vector<const reco::Candidate*> cands;
310  if(!cand) return cands;
311 
312  for(unsigned int i = 0; i < cand->numberOfDaughters(); i++) {
313 /*
314  std::cout << "DeltaR - " << partonFlavour << " ";
315  if (DeltaR(thePartonLV, cand->daughter(i)->p4()) > 0.7) std::cout << "(";
316  std::cout << cand->daughter(i)->pdgId() << ": " << DeltaR(thePartonLV, cand->daughter(i)->p4());
317  if (DeltaR(thePartonLV, cand->daughter(i)->p4()) > 0.7) std::cout << ")";
318  std::cout << std::endl;
319 */
320  if (DeltaR(thePartonLV, cand->daughter(i)->p4()) < 0.7) {
321  int pdgId = std::abs(cand->daughter(i)->pdgId());
322  int flavour = heaviestFlavour(pdgId);
323  if (flavour == partonFlavour ||
324  (flavour >= 10 && partonFlavour >= 10)) {
325 // std::cout << "<------- " << std::endl;
326  std::vector<const reco::Candidate*> newcands = findCandidates(cand->daughter(i), partonFlavour);
327 // std::cout << " ------->" << std::endl;
328  std::copy(newcands.begin(), newcands.end(), std::back_inserter(cands));
329  }
330  if (partonFlavour >= 10)
331  cands.push_back(cand->daughter(i));
332  }
333  }
334 
335  if (cands.empty() && std::abs(cand->pdgId()) > 110 &&
336  !(partonFlavour >= 4 && partonFlavour < 10 &&
337  heaviestFlavour(cand->pdgId()) < 4))
338  cands.push_back(cand);
339 
340  return cands;
341 }
342 
343 void JetFlavourIdentifier::fillLeptons(const std::vector<const reco::Candidate*> &cands, JetFlavour::Leptons &leptons, int rank, int flavour)
344 {
345  for(unsigned int j = 0; j < cands.size(); j++) {
346  for(unsigned int i = 0; i < cands[j]->numberOfDaughters(); i++) {
347  int pdgId = std::abs(cands[j]->daughter(i)->pdgId());
348 
349 // for(int z = 1; z <= rank; z *= 10) std::cout << " ------ ";
350 // std::cout << pdgId << std::endl;
351 
353  if (pdgId == 12)
354  leptons.electron += rank;
355  else if (pdgId == 14)
356  leptons.muon += rank;
357  else if (pdgId == 16)
358  leptons.tau += rank;
359  else {
360  int heaviest = heaviestFlavour(pdgId);
361  int heaviest_ = heaviest < 10 ? heaviest : 0;
362  if (!heaviest || (flavour < 4 ? (heaviest_ < 4) : (heaviest >= 4))) {
363  std::vector<const reco::Candidate*> newcands = findCandidates(cands[j]->daughter(i), heaviest);
364  if (pdgId <= 110) newcands.push_back(cands[j]->daughter(i));
365  fillLeptons(newcands, leptons, rank * 10, std::max(heaviest_, flavour));
366  }
367  }
368  }
369  }
370 }
371 
372 int JetFlavourIdentifier::heaviestFlavour(int pdgId)
373 {
374  int flavour = 0;
375 
376  pdgId = std::abs(pdgId) % 100000;
377  if (pdgId > 110) {
378  while(pdgId % 10 > 0 && pdgId % 10 < 6) {
379  pdgId /= 10;
380  if (pdgId % 10 > flavour)
381  flavour = pdgId % 10;
382  }
383  } else
384  flavour = pdgId;
385 
386  return flavour;
387 }
388 
389 
390 //define this as a plug-in
T getParameter(std::string const &) const
Handle< JetMatchedPartonsCollection > theTagByRef
int i
Definition: DBlmapReader.cc:9
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:252
transient_vector_type::const_iterator const_iterator
const GenParticleRef & nearest_status3() const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:464
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
bool exists(std::string const &parameterName) const
checks if a parameter exists
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
virtual size_type numberOfDaughters() const =0
number of daughters
int iEvent
Definition: GenABIO.cc:230
const GenParticleRef & algoDefinitionParton() const
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:120
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:244
const GenParticleRef & nearest_status2() const
virtual int pdgId() const =0
PDG identifier.
lepton info struct
Definition: JetFlavour.h:25
const GenParticleRef heaviest() const
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
const GenParticleRef & physicsDefinitionParton() const
RefToBaseProd< T > makeRefToBaseProdFrom(RefToBase< T > const &iRef, Event const &iEvent)
EDGetTokenT< JetMatchedPartonsCollection > sourceByReferToken_
int flavour(const Candidate &part)
Definition: pdgIdUtils.h:31
virtual const LorentzVector & p4() const =0
four-momentum Lorentz vector