CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
GenParticlesFromZsSelectorForMCEmbedding.cc
Go to the documentation of this file.
2 
5 
8 
9 #include <TMath.h>
10 
11 #include <vector>
12 
14 {
15  src_ = cfg.getParameter<edm::InputTag>("src");
16 
17  pdgIdsMothers_ = cfg.getParameter<vint>("pdgIdsMothers");
18  pdgIdsDaughters_ = cfg.getParameter<vint>("pdgIdsDaughters");
19 
20  maxDaughters_ = cfg.getParameter<int>("maxDaughters");
21  minDaughters_ = cfg.getParameter<int>("minDaughters");
22 
23  std::string before_or_afterFSR_string = cfg.getParameter<std::string>("before_or_afterFSR");
24  if ( before_or_afterFSR_string == "beforeFSR" ) before_or_afterFSR_ = kBeforeFSR;
25  else if ( before_or_afterFSR_string == "afterFSR" ) before_or_afterFSR_ = kAfterFSR;
26  else throw cms::Exception("Configuration")
27  << " Invalid Configuration Parameter 'before_or_afterFSR' = " << before_or_afterFSR_string << " !!\n";
28 
29  verbosity_ = ( cfg.exists("verbosity") ) ?
30  cfg.getParameter<int>("verbosity") : 0;
31 
32  produces<reco::GenParticleCollection>("");
33 }
34 
36 {
37 // nothing to be done yet...
38 }
39 
40 namespace
41 {
42  void findGenParticles(const reco::GenParticleCollection& genParticles,
43  int pdgIdMother, int pdgIdDaughter, std::vector<const reco::GenParticle*>& genParticlesFromZs)
44  {
45  for ( reco::GenParticleCollection::const_iterator genParticle = genParticles.begin();
46  genParticle != genParticles.end(); ++genParticle ) {
47  if ( TMath::Abs(genParticle->pdgId()) == pdgIdMother ) {
48  unsigned numDaughters = genParticle->numberOfDaughters();
49  for ( unsigned iDaughter = 0; iDaughter < numDaughters; ++iDaughter ) {
50  const reco::GenParticle* daughter = genParticle->daughterRef(iDaughter).get();
51  if ( TMath::Abs(daughter->pdgId()) == pdgIdDaughter ) genParticlesFromZs.push_back(daughter);
52  }
53  }
54  }
55  }
56 
57  void findGenParticles(const reco::GenParticleCollection& genParticles,
58  int pdgIdDaughter, int minDaughters, bool requireSubsequentEntries, std::vector<const reco::GenParticle*>& genParticlesFromZs)
59  {
60  int indexDaughterPlus = -1;
61  int indexDaughterMinus = -1;
62 
63  unsigned numGenParticles = genParticles.size();
64  for ( unsigned iGenParticle = 0; iGenParticle < numGenParticles; ++iGenParticle ) {
65  const reco::GenParticle& genParticle = genParticles[iGenParticle];
66 
67  if ( TMath::Abs(genParticle.pdgId()) == pdgIdDaughter ) {
68  if ( genParticle.charge() > 0 ) indexDaughterPlus = (int)iGenParticle;
69  else if ( genParticle.charge() < 0 ) indexDaughterMinus = (int)iGenParticle;
70 
71  if ( indexDaughterPlus != -1 && indexDaughterMinus != -1 ) {
72  if ( TMath::Abs(indexDaughterPlus - indexDaughterMinus) == 1 || (!requireSubsequentEntries) ) {
73  genParticlesFromZs.push_back(&genParticles[indexDaughterPlus]);
74  genParticlesFromZs.push_back(&genParticles[indexDaughterMinus]);
75  indexDaughterPlus = -1;
76  indexDaughterMinus = -1;
77  }
78  }
79  }
80 
81  if ( (int)genParticlesFromZs.size() >= minDaughters ) break;
82  }
83  }
84 
85  void findDaughters(const reco::GenParticle* mother, std::vector<const reco::GenParticle*>& daughters, int status)
86  {
87  unsigned numDaughters = mother->numberOfDaughters();
88  for ( unsigned iDaughter = 0; iDaughter < numDaughters; ++iDaughter ) {
89  const reco::GenParticle* daughter = mother->daughterRef(iDaughter).get();
90 
91  if ( status == -1 || daughter->status() == status ) daughters.push_back(daughter);
92 
93  findDaughters(daughter, daughters, status);
94  }
95  }
96 }
97 
99 {
100  if ( verbosity_ ) {
101  std::cout << "<GenParticlesFromZsSelectorForMCEmbedding::produce>:" << std::endl;
102  std::cout << " src = " << src_ << std::endl;
103  }
104 
106  evt.getByLabel(src_, genParticles);
107 
108  std::vector<const reco::GenParticle*> genParticlesFromZs_tmp;
109 
110 //--- check if HepEVT record contains any Z/gamma* --> l+ l- entries
111 //
112 // NOTE: iteration over mothers in the outer loop
113 // gives Z --> l+ l- decays priority over gamma --> e+ e-
114 //
115  for ( vint::const_iterator pdgIdMother = pdgIdsMothers_.begin();
116  pdgIdMother != pdgIdsMothers_.end(); ++pdgIdMother ) {
117  for ( vint::const_iterator pdgIdDaughter = pdgIdsDaughters_.begin();
118  pdgIdDaughter != pdgIdsDaughters_.end(); ++pdgIdDaughter ) {
119  if ( (int)genParticlesFromZs_tmp.size() < maxDaughters_ || maxDaughters_ == -1 )
120  findGenParticles(*genParticles, *pdgIdMother, *pdgIdDaughter, genParticlesFromZs_tmp);
121  }
122  }
123 
124 //--- check if HepEVT record contains l+ l- entries without Z/gamma* parent
125 //
126 // NOTE: in order to avoid ambiguities, give preference to l+ l- pairs
127 // which are subsequent in HepEVT record
128 //
129  if ( !((int)genParticlesFromZs_tmp.size() >= minDaughters_) ) {
130  for ( vint::const_iterator pdgIdDaughter = pdgIdsDaughters_.begin();
131  pdgIdDaughter != pdgIdsDaughters_.end(); ++pdgIdDaughter ) {
132  findGenParticles(*genParticles, *pdgIdDaughter, minDaughters_, true, genParticlesFromZs_tmp);
133  }
134  }
135  if ( !((int)genParticlesFromZs_tmp.size() >= minDaughters_) ) {
136  for ( vint::const_iterator pdgIdDaughter = pdgIdsDaughters_.begin();
137  pdgIdDaughter != pdgIdsDaughters_.end(); ++pdgIdDaughter ) {
138  findGenParticles(*genParticles, *pdgIdDaughter, minDaughters_, false, genParticlesFromZs_tmp);
139  }
140  }
141 
142  std::auto_ptr<reco::GenParticleCollection> genParticlesFromZs(new reco::GenParticleCollection());
143 
144  int idx = 0;
145  for ( std::vector<const reco::GenParticle*>::const_iterator genParticleFromZ_beforeFSR = genParticlesFromZs_tmp.begin();
146  genParticleFromZ_beforeFSR != genParticlesFromZs_tmp.end(); ++genParticleFromZ_beforeFSR ) {
147  if ( before_or_afterFSR_ == kBeforeFSR ) {
148  genParticlesFromZs->push_back(**genParticleFromZ_beforeFSR);
149  } else if ( before_or_afterFSR_ == kAfterFSR ) {
150  std::vector<const reco::GenParticle*> daughters;
151  findDaughters(*genParticleFromZ_beforeFSR, daughters, -1);
152  const reco::GenParticle* genParticleFromZ_afterFSR = (*genParticleFromZ_beforeFSR);
153  for ( std::vector<const reco::GenParticle*>::const_iterator daughter = daughters.begin();
154  daughter != daughters.end(); ++daughter ) {
155  if ( (*daughter)->pdgId() == (*genParticleFromZ_beforeFSR)->pdgId() &&
156  (*daughter)->energy() < genParticleFromZ_afterFSR->energy() ) genParticleFromZ_afterFSR = (*daughter);
157  }
158  if ( verbosity_ ) {
159  std::cout << "genParticleFromZ #" << idx << " (beforeFSR): Pt = " << (*genParticleFromZ_beforeFSR)->pt() << ","
160  << " eta = " << (*genParticleFromZ_beforeFSR)->eta() << ", phi = " << (*genParticleFromZ_beforeFSR)->phi() << std::endl;
161  std::cout << "genParticleFromZ #" << idx << " (afterFSR): Pt = " << genParticleFromZ_afterFSR->pt() << ","
162  << " eta = " << genParticleFromZ_afterFSR->eta() << ", phi = " << genParticleFromZ_afterFSR->phi() << std::endl;
163  }
164  genParticlesFromZs->push_back(*genParticleFromZ_afterFSR);
165  } else assert(0);
166  ++idx;
167  }
168 
169  evt.put(genParticlesFromZs);
170 }
171 
173 
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
T getParameter(std::string const &) const
tuple cfg
Definition: looper.py:293
daughters::value_type daughterRef(size_type i) const
reference to daughter at given position
virtual double energy() const final
energy
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
assert(m_qm.get())
bool exists(std::string const &parameterName) const
checks if a parameter exists
virtual double phi() const final
momentum azimuthal angle
virtual int status() const final
status word
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:121
const_iterator end() const
last daughter const_iterator
Definition: Candidate.h:146
virtual size_t numberOfDaughters() const
number of daughters
virtual int charge() const final
electric charge
Definition: LeafCandidate.h:91
T Abs(T a)
Definition: MathUtil.h:49
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:418
tuple idx
DEBUGGING if hasattr(process,&quot;trackMonIterativeTracking2012&quot;): print &quot;trackMonIterativeTracking2012 D...
virtual int pdgId() const final
PDG identifier.
const_iterator begin() const
first daughter const_iterator
Definition: Candidate.h:144
tuple cout
Definition: gather_cfg.py:145
virtual double eta() const final
momentum pseudorapidity
tuple status
Definition: mps_update.py:57
virtual double pt() const final
transverse momentum