CMS 3D CMS Logo

FlavorHistoryProducer.cc
Go to the documentation of this file.
1 
9 // -------------------------------------------------------------
10 // Identify the ancestry of the Quark
11 //
12 //
13 // Matrix Element:
14 // Status 3 parent with precisely 2 "grandparents" that
15 // is outside of the "initial" section (0-5) that has the
16 // same ID as the status 2 parton in question.
17 //
18 // Flavor excitation:
19 // Almost the same as the matrix element classification,
20 // but has only one outgoing parton product instead of two.
21 //
22 // Gluon splitting:
23 // Parent is a quark of a different flavor than the parton
24 // in question, or a gluon. Can come from either ISR or FSR.
25 //
26 // True decay:
27 // Decays from a resonance like top, Higgs, etc.
28 // -------------------------------------------------------------
29 
31 #include <string>
32 #include <vector>
33 #include <set>
34 #include <utility>
35 #include <algorithm>
36 
49 #include <fstream>
50 
53 
55 public:
58 
59 private:
61  void produce(edm::StreamID, edm::Event &e, const edm::EventSetup &) const override;
62 
63  void getAncestors(const reco::Candidate &c, std::vector<reco::Candidate const *> &moms) const;
64 
66  reco::Candidate const &parton) const;
67 
68  const edm::EDGetTokenT<reco::CandidateView> srcToken_; // GenParticles source collection name
69  const edm::EDGetTokenT<reco::CandidateView> matchedSrcToken_; // matched particles source collection name
70  const double matchDR_; // delta r to match matched particles
71  const int pdgIdToSelect_; // pdg of hf partons to select
72  const double ptMinParticle_; // minimum pt of the partons
73  const double ptMinShower_; // minimum pt of the shower
74  const double etaMaxParticle_; // max eta of the parton
75  const double etaMaxShower_; // max eta of the shower
76  const std::string flavorHistoryName_; // name to give flavor history
77  const bool verbose_; // verbose flag
78 };
79 
80 // #include "DataFormats/Common/interface/ValueMap.h"
81 // #include <iterators>
82 
83 // -------------------------------------------------------------
84 // Identify the ancestry of the Quark
85 //
86 //
87 // Matrix Element:
88 // Status 3 parent with precisely 2 "grandparents" that
89 // is outside of the "initial" section (0-5) that has the
90 // same ID as the status 2 parton in question.
91 //
92 // Flavor excitation:
93 // If we find only one outgoing parton.
94 //
95 // Gluon splitting:
96 // Parent is a quark of a different flavor than the parton
97 // in question, or a gluon.
98 // Can come from either ISR or FSR.
99 //
100 // True decay:
101 // Decays from a resonance like top, Higgs, etc.
102 // -------------------------------------------------------------
103 
104 using namespace std;
105 using namespace reco;
106 using namespace edm;
107 
109  : srcToken_(consumes<CandidateView>(p.getParameter<InputTag>("src"))),
110  matchedSrcToken_(consumes<CandidateView>(p.getParameter<InputTag>("matchedSrc"))),
111  matchDR_(p.getParameter<double>("matchDR")),
112  pdgIdToSelect_(p.getParameter<int>("pdgIdToSelect")),
113  ptMinParticle_(p.getParameter<double>("ptMinParticle")),
114  ptMinShower_(p.getParameter<double>("ptMinShower")),
115  etaMaxParticle_(p.getParameter<double>("etaMaxParticle")),
116  etaMaxShower_(p.getParameter<double>("etaMaxShower")),
117  flavorHistoryName_(p.getParameter<string>("flavorHistoryName")),
118  verbose_(p.getUntrackedParameter<bool>("verbose")) {
119  produces<FlavorHistoryEvent>(flavorHistoryName_);
120 }
121 
123  if (verbose_)
124  cout << "---------- Hello from FlavorHistoryProducer! -----" << endl;
125 
126  // Get a handle to the particle collection (OwnVector)
127  Handle<CandidateView> particlesViewH;
128  evt.getByToken(srcToken_, particlesViewH);
129 
130  Handle<CandidateView> matchedView;
131  evt.getByToken(matchedSrcToken_, matchedView);
132 
133  // Copy the View to an vector for easier iterator manipulation convenience
134  vector<const Candidate *> particles;
135  for (CandidateView::const_iterator p = particlesViewH->begin(); p != particlesViewH->end(); ++p) {
136  particles.push_back(&*p);
137  }
138 
139  // List of indices for the partons to add
140  vector<int> partonIndices;
141  // List of progenitors for those partons
142  vector<int> progenitorIndices;
143  // List of sisters for those partons
144  vector<int> sisterIndices;
145  // Flavor sources
146  vector<FlavorHistory::FLAVOR_T> flavorSources;
147 
148  // Make a new flavor history vector
149  auto flavorHistoryEvent = std::make_unique<FlavorHistoryEvent>();
150 
151  // ------------------------------------------------------------
152  // Loop over partons
153  // ------------------------------------------------------------
156  for (j = 0; j < j_max; ++j) {
157  // Get the candidate
158  const Candidate *p = particles[j];
159  // Set up indices that we'll need for the flavor history
161  vector<Candidate const *>::size_type progenitorIndex = j_max;
162  bool foundProgenitor = false;
163  FlavorHistory::FLAVOR_T flavorSource = FlavorHistory::FLAVOR_NULL;
164 
165  int idabs = std::abs((p)->pdgId());
166  int nDa = (p)->numberOfDaughters();
167 
168  // Check if we have a status 2 or 3 particle, which is a parton before the string.
169  // Only consider quarks.
170  if (idabs <= FlavorHistory::tQuarkId && p->status() == 2 && idabs == pdgIdToSelect_) {
171  // Ensure the parton in question has daughters
172  if (nDa > 0) {
173  // Ensure the parton passes some minimum kinematic cuts
174  if ((p)->pt() > ptMinShower_ && fabs((p)->eta()) < etaMaxShower_) {
175  if (verbose_)
176  cout << "--------------------------" << endl;
177  if (verbose_)
178  cout << "Processing particle " << j << ", status = " << p->status() << ", pdg id = " << p->pdgId() << endl;
179 
180  //Main (but simple) workhorse to get all ancestors
181  vector<Candidate const *> allParents;
182  getAncestors(*p, allParents);
183 
184  if (verbose_) {
185  cout << "Parents : " << endl;
186  vector<Candidate const *>::const_iterator iprint = allParents.begin(), iprintend = allParents.end();
187  for (; iprint != iprintend; ++iprint)
188  cout << " status = " << (*iprint)->status() << ", pdg id = " << (*iprint)->pdgId()
189  << ", pt = " << (*iprint)->pt() << endl;
190  }
191 
192  // -------------------------------------------------------------
193  // Identify the ancestry of the Quark
194  //
195  //
196  // Matrix Element:
197  // Status 3 parent with precisely 2 "grandparents" that
198  // is outside of the "initial" section (0-5) that has the
199  // same ID as the status 2 parton in question.
200  //
201  // Flavor excitation:
202  // If we find only one outgoing parton.
203  //
204  // Gluon splitting:
205  // Parent is a quark of a different flavor than the parton
206  // in question, or a gluon.
207  // Can come from either ISR or FSR.
208  //
209  // True decay:
210  // Decays from a resonance like top, Higgs, etc.
211  // -------------------------------------------------------------
212  vector<Candidate const *>::size_type a_size = allParents.size();
213 
214  //
215  // Loop over all the ancestors of this parton and find the progenitor.
216  //
217  for (vector<Candidate const *>::size_type i = 0; i < a_size && !foundProgenitor; ++i) {
218  const Candidate *aParent = allParents[i];
219  vector<Candidate const *>::const_iterator found = find(particles.begin(), particles.end(), aParent);
220 
221  // Get the index of the progenitor candidate
222  progenitorIndex = found - particles.begin();
223 
224  int aParentId = std::abs(aParent->pdgId());
225 
226  // -----------------------------------------------------------------------
227  // Here we examine particles that were produced after the collision
228  // -----------------------------------------------------------------------
229  if (progenitorIndex > 5) {
230  // Here is where we have a matrix element process.
231  // The signature is to have a status 3 parent with precisely 2 parents
232  // that is outside of the "initial" section (0-5) that has the
233  // same ID as the status 2 parton in question.
234  // NOTE: If this parton has no sister, then this will be classified as
235  // a flavor excitation. Often the only difference is that the matrix element
236  // cases have a sister, while the flavor excitation cases do not.
237  // If we do not find a sister below, this will be classified as a flavor
238  // excitation.
239  if (aParent->numberOfMothers() == 2 && aParent->pdgId() == p->pdgId() && aParent->status() == 3) {
240  if (verbose_)
241  cout << "Matrix Element progenitor" << endl;
242  if (flavorSource == FlavorHistory::FLAVOR_NULL)
243  flavorSource = FlavorHistory::FLAVOR_ME;
244  foundProgenitor = true;
245  }
246  // Here we have a true decay. Parent is not a quark or a gluon.
247  else if ((aParentId > pdgIdToSelect_ && aParentId < FlavorHistory::gluonId) ||
248  aParentId > FlavorHistory::gluonId) {
249  if (verbose_)
250  cout << "Flavor decay progenitor" << endl;
251  if (flavorSource == FlavorHistory::FLAVOR_NULL)
252  flavorSource = FlavorHistory::FLAVOR_DECAY;
253  foundProgenitor = true;
254  }
255  } // end if progenitorIndex > 5
256  } // end loop over parents
257 
258  // Otherwise, classify it as gluon splitting. Take the first status 3 parton with 1 parent
259  // that you see as the progenitor
260  if (!foundProgenitor) {
261  if (flavorSource == FlavorHistory::FLAVOR_NULL)
262  flavorSource = FlavorHistory::FLAVOR_GS;
263  // Now get the parton with only one parent (the proton) and that is the progenitor
264  for (vector<Candidate const *>::size_type i = 0; i < a_size && !foundProgenitor; ++i) {
265  const Candidate *aParent = allParents[i];
266  vector<Candidate const *>::const_iterator found = find(particles.begin(), particles.end(), aParent);
267  // Get the index of the progenitor candidate
268  progenitorIndex = found - particles.begin();
269 
270  if (aParent->numberOfMothers() == 1 && aParent->status() == 3) {
271  foundProgenitor = true;
272  }
273  } // end loop over parents
274  } // end if we haven't found a progenitor, and if we haven't found a status 3 parton of the same flavor
275  // (i.e. end if examining gluon splitting)
276  } // End if this parton passes some minimal kinematic cuts
277  } // End if this particle is status 2 (has strings as daughters)
278 
279  // Make sure we've actually found a progenitor
280  if (!foundProgenitor)
281  progenitorIndex = j_max;
282 
283  // We've found the particle, add to the list
284 
285  partonIndices.push_back(partonIndex);
286  progenitorIndices.push_back(progenitorIndex);
287  flavorSources.push_back(flavorSource);
288  sisterIndices.push_back(-1); // set below
289 
290  } // End if this particle was a status==2 parton
291  } // end loop over particles
292 
293  // Now add sisters.
294  // Also if the event is preliminarily classified as "matrix element", check to
295  // make sure that they have a sister. If not, it is flavor excitation.
296 
297  if (verbose_)
298  cout << "Making sisters" << endl;
299  // First make sure nothing went terribly wrong:
300  if (partonIndices.size() == progenitorIndices.size() && !partonIndices.empty()) {
301  // Now loop over the candidates
302  for (unsigned int ii = 0; ii < partonIndices.size(); ++ii) {
303  // Get the iith particle
304  const Candidate *candi = particles[partonIndices[ii]];
305  if (verbose_)
306  cout << " Sister candidate particle 1: " << ii << ", pdgid = " << candi->pdgId()
307  << ", status = " << candi->status() << endl;
308  // Get the iith progenitor
309  // Now loop over the other flavor history candidates and
310  // attempt to find a sister to this one
311  for (unsigned int jj = 0; jj < partonIndices.size(); ++jj) {
312  if (ii != jj) {
313  const Candidate *candj = particles[partonIndices[jj]];
314  if (verbose_)
315  cout << " Sister candidate particle 2: " << jj << ", pdgid = " << candj->pdgId()
316  << ", status = " << candj->status() << endl;
317  // They should be opposite in pdgid and have the same status, and
318  // the same progenitory.
319  if (candi->pdgId() == -1 * candj->pdgId() && candi->status() == candj->status()) {
320  sisterIndices[ii] = partonIndices[jj];
321  if (verbose_)
322  cout << "Parton " << partonIndices[ii] << " has sister " << sisterIndices[ii] << endl;
323  }
324  }
325  }
326 
327  // Here, ensure that there is a sister. Otherwise this is "flavor excitation"
328  if (sisterIndices[ii] < 0) {
329  if (verbose_)
330  cout << "No sister. Classified as flavor excitation" << endl;
331  flavorSources[ii] = FlavorHistory::FLAVOR_EXC;
332  }
333 
334  if (verbose_)
335  cout << "Getting matched jet" << endl;
336  // Get the closest match in the matched object collection
337  CandidateView::const_iterator matched = getClosestJet(matchedView, *candi);
338  CandidateView::const_iterator matchedBegin = matchedView->begin();
339  CandidateView::const_iterator matchedEnd = matchedView->end();
340 
341  if (verbose_)
342  cout << "Getting sister jet" << endl;
343  // Get the closest sister in the matched object collection, if sister is found
345  (sisterIndices[ii] >= 0 && static_cast<unsigned int>(sisterIndices[ii]) < particles.size())
346  ? getClosestJet(matchedView, *particles[sisterIndices[ii]])
347  : matchedEnd;
348 
349  if (verbose_)
350  cout << "Making matched shallow clones : ";
351  ShallowClonePtrCandidate matchedCand;
352  if (matched != matchedEnd) {
353  if (verbose_)
354  cout << " found" << endl;
355  matchedCand = ShallowClonePtrCandidate(CandidatePtr(matchedView, matched - matchedBegin));
356  } else {
357  if (verbose_)
358  cout << " NOT found" << endl;
359  }
360 
361  if (verbose_)
362  cout << "Making sister shallow clones : ";
363  ShallowClonePtrCandidate sisterCand;
364  if (sister != matchedEnd) {
365  if (verbose_)
366  cout << " found" << endl;
367  sisterCand = ShallowClonePtrCandidate(CandidatePtr(matchedView, sister - matchedBegin));
368  } else {
369  if (verbose_)
370  cout << " NOT found" << endl;
371  }
372 
373  if (verbose_)
374  cout << "Making history object" << endl;
375  // Now create the flavor history object
376  FlavorHistory history(flavorSources[ii],
377  particlesViewH,
378  partonIndices[ii],
379  progenitorIndices[ii],
380  sisterIndices[ii],
381  matchedCand,
382  sisterCand);
383  if (verbose_)
384  cout << "Adding flavor history : " << history << endl;
385  flavorHistoryEvent->push_back(history);
386  }
387  }
388 
389  // If we got any partons, cache them and then write them out
390  if (flavorHistoryEvent->size() > 0) {
391  // Calculate some nice variables for FlavorHistoryEvent
392  if (verbose_)
393  cout << "Caching flavor history event" << endl;
394  flavorHistoryEvent->cache();
395 
396  if (verbose_) {
397  cout << "Outputting pdg id = " << pdgIdToSelect_ << " with nelements = " << flavorHistoryEvent->size() << endl;
398  vector<FlavorHistory>::const_iterator i = flavorHistoryEvent->begin(), iend = flavorHistoryEvent->end();
399  for (; i != iend; ++i) {
400  cout << *i << endl;
401  }
402  }
403  }
404 
405  // Now add the flavor history to the event record
406  evt.put(std::move(flavorHistoryEvent), flavorHistoryName_);
407 }
408 
409 // Helper function to get all ancestors of this candidate
410 void FlavorHistoryProducer::getAncestors(const Candidate &c, vector<Candidate const *> &moms) const {
411  if (c.numberOfMothers() == 1) {
412  const Candidate *dau = &c;
413  const Candidate *mom = c.mother();
414  while (dau->numberOfMothers() != 0) {
415  moms.push_back(dau);
416  dau = mom;
417  mom = dau->mother();
418  }
419  }
420 }
421 
423  reco::Candidate const &parton) const {
424  double dr = matchDR_;
425  CandidateView::const_iterator j = pJets->begin(), jend = pJets->end();
426  CandidateView::const_iterator bestJet = pJets->end();
427  for (; j != jend; ++j) {
428  double dri = deltaR(parton.p4(), j->p4());
429  if (dri < dr) {
430  dr = dri;
431  bestJet = j;
432  }
433  }
434  return bestJet;
435 }
436 
438 
FlavorHistoryProducer(const edm::ParameterSet &)
constructor
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:133
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
virtual int status() const =0
status word
void produce(edm::StreamID, edm::Event &e, const edm::EventSetup &) const override
process one event
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
void getAncestors(const reco::Candidate &c, std::vector< reco::Candidate const *> &moms) const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:540
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
uint16_t size_type
const edm::EDGetTokenT< reco::CandidateView > srcToken_
const edm::EDGetTokenT< reco::CandidateView > matchedSrcToken_
reco::CandidateView::const_iterator getClosestJet(edm::Handle< reco::CandidateView > const &pJets, reco::Candidate const &parton) const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
ii
Definition: cuy.py:589
virtual int pdgId() const =0
PDG identifier.
edm::Ptr< Candidate > CandidatePtr
persistent reference to an object in a collection of Candidate objects
Definition: CandidateFwd.h:25
fixed size matrix
HLT enums.
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:86
const std::string flavorHistoryName_
const_iterator begin() const
const_iterator end() const
def move(src, dest)
Definition: eostools.py:511
virtual const LorentzVector & p4() const =0
four-momentum Lorentz vector