CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HepMCValidationHelper.cc
Go to the documentation of this file.
3 #include <iostream>
4 #include <cassert>
5 #include <limits>
6 #include "TLorentzVector.h"
7 
8 //#define DEBUG_HepMCValidationHelper
9 
10 namespace HepMCValidationHelper {
11  void findFSRPhotons(const std::vector<const HepMC::GenParticle*>& leptons,
12  const HepMC::GenEvent* all,
13  double deltaRcut,
14  std::vector<const HepMC::GenParticle*>& fsrphotons){
15 
16  std::vector<const HepMC::GenParticle*> status1;
17  allStatus1(all, status1);
18  findFSRPhotons(leptons, status1, deltaRcut, fsrphotons);
19  }
20 
21  void findFSRPhotons(const std::vector<const HepMC::GenParticle*>& leptons,
22  const std::vector<const HepMC::GenParticle*>& all,
23  double deltaRcut,
24  std::vector<const HepMC::GenParticle*>& fsrphotons){
25 
26  //find all status 1 photons
27  std::vector<const HepMC::GenParticle*> allphotons;
28  for (unsigned int i = 0; i < all.size(); ++i){
29  if (all[i]->status()==1 && all[i]->pdg_id()==22) allphotons.push_back(all[i]);
30  }
31 
32  //loop over the photons and check the distance wrt the leptons
33  for (unsigned int ipho = 0; ipho < allphotons.size(); ++ipho){
34  bool close = false;
35  for (unsigned int ilep = 0; ilep < leptons.size(); ++ilep){
36  if (deltaR(allphotons[ipho]->momentum(), leptons[ilep]->momentum()) < deltaRcut){
37  close = true;
38  break;
39  }
40  }
41  if (close) fsrphotons.push_back(allphotons[ipho]);
42  }
43  }
44 
45  //returns true if a status 3 particle is a tau or if a status 1 particle is either an electron or a neutrino
47  int status = part->status();
48  unsigned int pdg_id = abs(part->pdg_id());
49  if (status == 2) return pdg_id == 15;
50  else return status == 1 && (pdg_id == 11 || pdg_id == 13 ) ;
51  }
52 
53  //returns true if a status 1 particle is a neutrino
55  int status = part->status();
56  unsigned int pdg_id = abs(part->pdg_id());
57  return status == 1 && (pdg_id == 12 || pdg_id == 14 || pdg_id == 16) ;
58  }
59 
60  //returns true is status 3 particle is tau
62  return part->status() == 2 && abs(part->pdg_id()) == 15;
63  }
64 /*
65  void getTaus(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& taus){
66  for (HepMC::GenEvent::particle_const_iterator iter = all->particles_begin(); iter != all->particles_end(); ++iter){
67  if (abs((*iter)->pdg_id()) == 15) taus.push_back(*iter);
68  }
69  }
70 */
71  // get all status 1 particles
72  void allStatus1(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& status1){
73  for (HepMC::GenEvent::particle_const_iterator iter = all->particles_begin(); iter != all->particles_end(); ++iter){
74  if ((*iter)->status() == 1) status1.push_back(*iter);
75  }
76  }
77 
78  void allStatus2(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& status1){
79  for (HepMC::GenEvent::particle_const_iterator iter = all->particles_begin(); iter != all->particles_end(); ++iter){
80  if ((*iter)->status() == 2) status1.push_back(*iter);
81  }
82  }
83 
84  void allStatus3(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& status1){
85  for (HepMC::GenEvent::particle_const_iterator iter = all->particles_begin(); iter != all->particles_end(); ++iter){
86  if ((*iter)->status() == 3) status1.push_back(*iter);
87  }
88  }
89 
90  void findDescendents(const HepMC::GenParticle* a, std::vector<const HepMC::GenParticle*>& descendents){
91  HepMC::GenVertex* decayVertex = a->end_vertex();
92  if (!decayVertex) return;
93  HepMC::GenVertex::particles_out_const_iterator ipart;
94  for (ipart = decayVertex->particles_out_const_begin(); ipart != decayVertex->particles_out_const_end(); ++ipart ){
95  if ((*ipart)->status() == 1) descendents.push_back(*ipart);
96  else findDescendents(*ipart, descendents);
97  }
98  }
99 
100  void removeIsolatedLeptons(const HepMC::GenEvent* all, double deltaRcut, double sumPtCut, std::vector<const HepMC::GenParticle*>& pruned) {
101  //get all status 1 particles
102  std::vector<const HepMC::GenParticle*> status1;
103  allStatus1(all, status1);
104  std::vector<const HepMC::GenParticle*> toRemove;
105  //loop on all particles and find candidates to be isolated
106  for (unsigned int i = 0; i < status1.size(); ++i){
107  //if it is a neutrino is a charged lepton (not a tau) this is a candidate to be isolated
108  if (isNeutrino(status1[i]) || (isChargedLepton(status1[i]) && !isTau(status1[i]))){
109  //list of particles not to be considered in the isolation computation.
110  //this includes the particle to be isolated and the fsr photons in case of charged lepton
111  std::vector<const HepMC::GenParticle*> leptons;
112  leptons.push_back(status1[i]);
113  std::vector<const HepMC::GenParticle*> removedForIsolation;
114  removedForIsolation.push_back(status1[i]);
115  if (isChargedLepton(status1[i])) findFSRPhotons(leptons, status1, deltaRcut, removedForIsolation);
116 #ifdef DEBUG_HepMCValidationHelper
117  //std::cout << removedForIsolation.size() << " particles to be removed for isolation calculation " << std::endl;
118 #endif
119  //create vector of particles to compute isolation (removing removedForIsolation);
120  std::vector<const HepMC::GenParticle*> forIsolation;
121  std::vector<const HepMC::GenParticle*>::iterator iiso;
122  for (iiso = status1.begin(); iiso != status1.end(); ++iiso){
123  std::vector<const HepMC::GenParticle*>::const_iterator iremove;
124  bool marked = false;
125  for ( iremove = removedForIsolation.begin(); iremove != removedForIsolation.end(); ++iremove){
126  if ((*iiso)->barcode() == (*iremove)->barcode()) {
127 #ifdef DEBUG_HepMCValidationHelper
128  //std::cout << "removing particle " << **iiso << " from the list of particles to compute isolation" << std::endl;
129 #endif
130  marked = true;
131  break;
132  }
133  }
134  if (!marked) forIsolation.push_back(*iiso);
135  }
136  //now compute isolation
137  double sumIso = 0;
138  for (iiso = forIsolation.begin(); iiso < forIsolation.end(); ++iiso){
139  if (deltaR(leptons.front()->momentum(), (*iiso)->momentum()) < deltaRcut){
140  sumIso += (*iiso)->momentum().perp();
141  }
142  }
143  //if isolated remove from the pruned list
144  if (sumIso < sumPtCut) {
145 #ifdef DEBUG_HepMCValidationHelper
146  std::cout << "particle " << *status1[i] << " is considered isolated, with sumPt " << sumIso << std::endl;
147 #endif
148  toRemove.insert(toRemove.end(), removedForIsolation.begin(), removedForIsolation.end());
149  }
150 #ifdef DEBUG_HepMCValidationHelper
151  else {
152  std::cout << "NOT isolated! " << *status1[i] << " is considered not isolated, with sumPt " << sumIso << std::endl;
153  }
154 #endif
155  }
156  }
157  //at this point we have taken care of the electrons and muons, but pruned could still contain the decay products of isolated taus,
158  //we want to remove these as well
159  std::vector<const HepMC::GenParticle*> status2;
160  allStatus2(all, status2);
161  std::vector<const HepMC::GenParticle*> taus;
162  //getTaus(all, taus);
163  for (unsigned int i = 0; i < status2.size(); ++i){
164  if (isTau(status2[i])) {
165  //check the list we have already for duplicates
166  //there use to be duplicates in some generators (sherpa)
167  bool duplicate = false;
168  TLorentzVector taumomentum(status2[i]->momentum().x(), status2[i]->momentum().y(), status2[i]->momentum().z(), status2[i]->momentum().t());
169  for (unsigned int j = 0; j < taus.size(); ++j){
170  //compare momenta
171  TLorentzVector othermomentum(taus[j]->momentum().x(), taus[j]->momentum().y(), taus[j]->momentum().z(), taus[j]->momentum().t());
172  othermomentum-=taumomentum;
173  if (status2[i]->pdg_id() == taus[j]->pdg_id() &&
174  othermomentum.E() < 0.1 &&//std::numeric_limits<float>::epsilon() &&
175  othermomentum.P() < 0.1 ){ //std::numeric_limits<float>::epsilon()){
176  duplicate = true;
177  break;
178  }
179  }
180  if (!duplicate) taus.push_back(status2[i]);
181  }
182  }
183  //loop over the taus, find the descendents, remove all these from the list of particles to compute isolation
184  for (unsigned int i = 0; i < taus.size(); ++i){
185  std::vector<const HepMC::GenParticle*> taudaughters;
186  findDescendents(taus[i], taudaughters);
187  assert(taudaughters.size()>0);
188  const HepMC::FourVector& taumom = taus[i]->momentum();
189  //remove the daughters from the list of particles to compute isolation
190  std::vector<const HepMC::GenParticle*> forIsolation;
191  std::vector<const HepMC::GenParticle*>::iterator iiso;
192  for (iiso = status1.begin(); iiso < status1.end(); ++iiso){
193  bool marked = false;
194  std::vector<const HepMC::GenParticle*>::const_iterator iremove;
195  for ( iremove = taudaughters.begin(); iremove != taudaughters.end(); ++iremove){
196  if ((*iiso)->barcode() == (*iremove)->barcode()) {
197 #ifdef DEBUG_HepMCValidationHelper
198 // std::cout << "removing particle " << **iiso << " from the list of particles to compute isolation because it comes from a tau" << std::endl;
199 #endif
200  marked = true;
201  break;
202  }
203  if (!marked) forIsolation.push_back(*iiso);
204  }
205  }
206  //no compute isolation wrt the status 2 tau direction
207  double sumIso = 0;
208  for (iiso = forIsolation.begin(); iiso < forIsolation.end(); ++iiso){
209  if (deltaR(taumom, (*iiso)->momentum()) < deltaRcut){
210  sumIso += (*iiso)->momentum().perp();
211  }
212  }
213  //if isolated remove the tau daughters from the pruned list
214  if (sumIso < sumPtCut) {
215 #ifdef DEBUG_HepMCValidationHelper
216  std::cout << "particle " << *taus[i] << " is considered isolated, with sumPt " << sumIso << std::endl;
217 #endif
218  toRemove.insert(toRemove.end(), taudaughters.begin(), taudaughters.end());
219  }
220  }
221 
222  //now actually remove
223  pruned.clear();
224  for (unsigned int i = 0; i < status1.size(); ++i){
225  bool marked = false;
226  std::vector<const HepMC::GenParticle*>::const_iterator iremove;
227  for ( iremove = toRemove.begin(); iremove != toRemove.end(); ++iremove){
228  if (status1[i]->barcode() == (*iremove)->barcode()){
229  marked = true;
230  break;
231  }
232  }
233  if (!marked) pruned.push_back(status1[i]);
234  }
235 
236 #ifdef DEBUG_HepMCValidationHelper
237  std::cout<< "list of remaining particles:" <<std::endl;
238  for (unsigned int i = 0; i < pruned.size(); ++i){
239  std::cout << *pruned[i] << std::endl;
240  }
241 #endif
242  }
243 
244  //get all visible status1 particles
245  void allVisibleParticles(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& visible) {
246  std::vector<const HepMC::GenParticle*> status1;
247  visible.clear();
248  allStatus1(all, status1);
249  for (unsigned int i = 0; i < status1.size(); ++i){
250  if (!isNeutrino(status1[i])) visible.push_back(status1[i]);
251  }
252  }
253 
254  //compute generated met
255  TLorentzVector genMet(const HepMC::GenEvent* all, double etamin, double etamax){
256  std::vector<const HepMC::GenParticle*> visible;
257  allVisibleParticles(all, visible);
258  TLorentzVector momsum(0., 0., 0., 0.);
259  for (unsigned int i = 0; i < visible.size(); ++i){
260  if (visible[i]->momentum().eta() > etamin && visible[i]->momentum().eta() < etamax) {
261  TLorentzVector mom(visible[i]->momentum().x(), visible[i]->momentum().y(), visible[i]->momentum().z(), visible[i]->momentum().t());
262  momsum += mom;
263  }
264  }
265  TLorentzVector met(-momsum.Px(), -momsum.Py(), 0., momsum.Pt());
266  return met;
267  }
268 }
int i
Definition: DBlmapReader.cc:9
bool isTau(const HepMC::GenParticle *part)
void allStatus3(const HepMC::GenEvent *all, std::vector< const HepMC::GenParticle * > &status3)
void findFSRPhotons(const std::vector< const HepMC::GenParticle * > &leptons, const std::vector< const HepMC::GenParticle * > &all, double deltaR, std::vector< const HepMC::GenParticle * > &photons)
#define abs(x)
Definition: mlp_lapack.h:159
bool isChargedLepton(const HepMC::GenParticle *part)
void allStatus1(const HepMC::GenEvent *all, std::vector< const HepMC::GenParticle * > &status1)
float float float z
void removeIsolatedLeptons(const HepMC::GenEvent *all, double deltaR, double sumPt, std::vector< const HepMC::GenParticle * > &pruned)
void allStatus2(const HepMC::GenEvent *all, std::vector< const HepMC::GenParticle * > &status2)
int j
Definition: DBlmapReader.cc:9
void allVisibleParticles(const HepMC::GenEvent *all, std::vector< const HepMC::GenParticle * > &visible)
void findDescendents(const HepMC::GenParticle *a, std::vector< const HepMC::GenParticle * > &descendents)
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
part
Definition: HCALResponse.h:20
TLorentzVector genMet(const HepMC::GenEvent *all, double etamin=-9999., double etamax=9999.)
bool isNeutrino(const HepMC::GenParticle *part)
double a
Definition: hdecay.h:121
tuple cout
Definition: gather_cfg.py:121
Definition: DDAxes.h:10
tuple status
Definition: ntuplemaker.py:245