EventGenerator/src/HepMCValidationHelper.cc

0001 #include "Validation/EventGenerator/interface/HepMCValidationHelper.h"
0002 #include "DataFormats/Math/interface/deltaR.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 #include <iostream>
0006 #include <limits>
0007 #include "TLorentzVector.h"
0008
0009 //#define DEBUG_HepMCValidationHelper
0010
0011 namespace HepMCValidationHelper {
0012   void findFSRPhotons(const std::vector<const HepMC::GenParticle*>& leptons,
0013                       const HepMC::GenEvent* all,
0014                       double deltaRcut,
0015                       std::vector<const HepMC::GenParticle*>& fsrphotons) {
0016     std::vector<const HepMC::GenParticle*> status1;
0017     allStatus1(all, status1);
0018     findFSRPhotons(leptons, status1, deltaRcut, fsrphotons);
0019   }
0020
0021   void findFSRPhotons(const std::vector<const HepMC::GenParticle*>& leptons,
0022                       const std::vector<const HepMC::GenParticle*>& all,
0023                       double deltaRcut,
0024                       std::vector<const HepMC::GenParticle*>& fsrphotons) {
0025     //find all status 1 photons
0026     std::vector<const HepMC::GenParticle*> allphotons;
0027     for (unsigned int i = 0; i < all.size(); ++i) {
0028       if (all[i]->status() == 1 && all[i]->pdg_id() == 22)
0029         allphotons.push_back(all[i]);
0030     }
0031
0032     //loop over the photons and check the distance wrt the leptons
0033     for (unsigned int ipho = 0; ipho < allphotons.size(); ++ipho) {
0034       bool close = false;
0035       for (unsigned int ilep = 0; ilep < leptons.size(); ++ilep) {
0036         if (deltaR(allphotons[ipho]->momentum(), leptons[ilep]->momentum()) < deltaRcut) {
0037           close = true;
0038           break;
0039         }
0040       }
0041       if (close)
0042         fsrphotons.push_back(allphotons[ipho]);
0043     }
0044   }
0045
0046   //returns true if a status 3 particle is a tau or if a status 1 particle is either an electron or a neutrino
0047   bool isChargedLepton(const HepMC::GenParticle* part) {
0048     int status = part->status();
0049     unsigned int pdg_id = abs(part->pdg_id());
0050     if (status == 2)
0051       return pdg_id == 15;
0052     else
0053       return status == 1 && (pdg_id == 11 || pdg_id == 13);
0054   }
0055
0056   //returns true if a status 1 particle is a neutrino
0057   bool isNeutrino(const HepMC::GenParticle* part) {
0058     int status = part->status();
0059     unsigned int pdg_id = abs(part->pdg_id());
0060     return status == 1 && (pdg_id == 12 || pdg_id == 14 || pdg_id == 16);
0061   }
0062
0063   //returns true is status 3 particle is tau
0064   bool isTau(const HepMC::GenParticle* part) { return part->status() == 2 && abs(part->pdg_id()) == 15; }
0065   /*
0066   void getTaus(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& taus){
0067     for (HepMC::GenEvent::particle_const_iterator iter = all->particles_begin(); iter != all->particles_end(); ++iter){
0068       if (abs((*iter)->pdg_id()) == 15) taus.push_back(*iter);
0069     }
0070   }
0071 */
0072   // get all status 1 particles
0073   void allStatus1(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& status1) {
0074     for (HepMC::GenEvent::particle_const_iterator iter = all->particles_begin(); iter != all->particles_end(); ++iter) {
0075       if ((*iter)->status() == 1)
0076         status1.push_back(*iter);
0077     }
0078   }
0079
0080   void allStatus2(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& status1) {
0081     for (HepMC::GenEvent::particle_const_iterator iter = all->particles_begin(); iter != all->particles_end(); ++iter) {
0082       if ((*iter)->status() == 2)
0083         status1.push_back(*iter);
0084     }
0085   }
0086
0087   void allStatus3(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& status1) {
0088     for (HepMC::GenEvent::particle_const_iterator iter = all->particles_begin(); iter != all->particles_end(); ++iter) {
0089       if ((*iter)->status() == 3)
0090         status1.push_back(*iter);
0091     }
0092   }
0093
0094   void findDescendents(const HepMC::GenParticle* a, std::vector<const HepMC::GenParticle*>& descendents) {
0095     HepMC::GenVertex* decayVertex = a->end_vertex();
0096     if (!decayVertex)
0097       return;
0098     HepMC::GenVertex::particles_out_const_iterator ipart;
0099     for (ipart = decayVertex->particles_out_const_begin(); ipart != decayVertex->particles_out_const_end(); ++ipart) {
0100       if ((*ipart)->status() == 1)
0101         descendents.push_back(*ipart);
0102       else
0103         findDescendents(*ipart, descendents);
0104     }
0105   }
0106
0107   void removeIsolatedLeptons(const HepMC::GenEvent* all,
0108                              double deltaRcut,
0109                              double sumPtCut,
0110                              std::vector<const HepMC::GenParticle*>& pruned) {
0111     //get all status 1 particles
0112     std::vector<const HepMC::GenParticle*> status1;
0113     allStatus1(all, status1);
0114     std::vector<const HepMC::GenParticle*> toRemove;
0115     //loop on all particles and find candidates to be isolated
0116     for (unsigned int i = 0; i < status1.size(); ++i) {
0117       //if it is a neutrino is a charged lepton (not a tau) this is a candidate to be isolated
0118       if (isNeutrino(status1[i]) || (isChargedLepton(status1[i]) && !isTau(status1[i]))) {
0119         //list of particles not to be considered in the isolation computation.
0120         //this includes the particle to be isolated and the fsr photons in case of charged lepton
0121         std::vector<const HepMC::GenParticle*> leptons;
0122         leptons.push_back(status1[i]);
0123         std::vector<const HepMC::GenParticle*> removedForIsolation;
0124         removedForIsolation.push_back(status1[i]);
0125         if (isChargedLepton(status1[i]))
0126           findFSRPhotons(leptons, status1, deltaRcut, removedForIsolation);
0127 #ifdef DEBUG_HepMCValidationHelper
0128           //std::cout << removedForIsolation.size() << " particles to be removed for isolation calculation " << std::endl;
0129 #endif
0130         //create vector of particles to compute isolation (removing removedForIsolation);
0131         std::vector<const HepMC::GenParticle*> forIsolation;
0132         std::vector<const HepMC::GenParticle*>::iterator iiso;
0133         for (iiso = status1.begin(); iiso != status1.end(); ++iiso) {
0134           std::vector<const HepMC::GenParticle*>::const_iterator iremove;
0135           bool marked = false;
0136           for (iremove = removedForIsolation.begin(); iremove != removedForIsolation.end(); ++iremove) {
0137             if ((*iiso)->barcode() == (*iremove)->barcode()) {
0138 #ifdef DEBUG_HepMCValidationHelper
0139               //std::cout << "removing particle " << **iiso << " from the list of particles to compute isolation" << std::endl;
0140 #endif
0141               marked = true;
0142               break;
0143             }
0144           }
0145           if (!marked)
0146             forIsolation.push_back(*iiso);
0147         }
0148         //now compute isolation
0149         double sumIso = 0;
0150         for (iiso = forIsolation.begin(); iiso < forIsolation.end(); ++iiso) {
0151           if (deltaR(leptons.front()->momentum(), (*iiso)->momentum()) < deltaRcut) {
0152             sumIso += (*iiso)->momentum().perp();
0153           }
0154         }
0155         //if isolated remove from the pruned list
0156         if (sumIso < sumPtCut) {
0157 #ifdef DEBUG_HepMCValidationHelper
0158           std::cout << "particle " << *status1[i] << " is considered isolated, with sumPt " << sumIso << std::endl;
0159 #endif
0160           toRemove.insert(toRemove.end(), removedForIsolation.begin(), removedForIsolation.end());
0161         }
0162 #ifdef DEBUG_HepMCValidationHelper
0163         else {
0164           std::cout << "NOT isolated! " << *status1[i] << " is considered not isolated, with sumPt " << sumIso
0165                     << std::endl;
0166         }
0167 #endif
0168       }
0169     }
0170     //at this point we have taken care of the electrons and muons, but pruned could  still contain the decay products of isolated taus,
0171     //we want to remove these as well
0172     std::vector<const HepMC::GenParticle*> status2;
0173     allStatus2(all, status2);
0174     std::vector<const HepMC::GenParticle*> taus;
0175     //getTaus(all, taus);
0176     for (unsigned int i = 0; i < status2.size(); ++i) {
0177       if (isTau(status2[i])) {
0178         //check the list we have already for duplicates
0179         //there use to be duplicates in some generators (sherpa)
0180         bool duplicate = false;
0181         TLorentzVector taumomentum(status2[i]->momentum().x(),
0182                                    status2[i]->momentum().y(),
0183                                    status2[i]->momentum().z(),
0184                                    status2[i]->momentum().t());
0185         for (unsigned int j = 0; j < taus.size(); ++j) {
0186           //compare momenta
0187           TLorentzVector othermomentum(
0188               taus[j]->momentum().x(), taus[j]->momentum().y(), taus[j]->momentum().z(), taus[j]->momentum().t());
0189           othermomentum -= taumomentum;
0190           if (status2[i]->pdg_id() == taus[j]->pdg_id() &&
0191               othermomentum.E() < 0.1 &&  //std::numeric_limits<float>::epsilon() &&
0192               othermomentum.P() < 0.1) {  //std::numeric_limits<float>::epsilon()){
0193             duplicate = true;
0194             break;
0195           }
0196         }
0197         if (!duplicate)
0198           taus.push_back(status2[i]);
0199       }
0200     }
0201     //loop over the taus, find the descendents, remove all these from the list of particles to compute isolation
0202     for (unsigned int i = 0; i < taus.size(); ++i) {
0203       std::vector<const HepMC::GenParticle*> taudaughters;
0204       findDescendents(taus[i], taudaughters);
0205       if (taudaughters.empty()) {
0206         std::ostringstream ss;
0207         auto vertex = taus[i]->end_vertex();
0208         if (vertex) {
0209           ss << "( " << vertex->point3d().x() << " " << vertex->point3d().y() << " " << vertex->point3d().z() << " )";
0210         } else {
0211           ss << "( did not decay )";
0212         }
0213         throw cms::Exception("NoTauDaugters")
0214             << " HepMCValidationHelper found no daughters for Tau within index " << i << " and info \n"
0215             << *taus[i] << " decay point " << ss.str()
0216             << "\n  This should not be able to happen and needs to be fixed.";
0217       }
0218       const HepMC::FourVector& taumom = taus[i]->momentum();
0219       //remove the daughters from the list of particles to compute isolation
0220       std::vector<const HepMC::GenParticle*> forIsolation;
0221       std::vector<const HepMC::GenParticle*>::iterator iiso;
0222       for (iiso = status1.begin(); iiso < status1.end(); ++iiso) {
0223         bool marked = false;
0224         std::vector<const HepMC::GenParticle*>::const_iterator iremove;
0225         for (iremove = taudaughters.begin(); iremove != taudaughters.end(); ++iremove) {
0226           if ((*iiso)->barcode() == (*iremove)->barcode()) {
0227 #ifdef DEBUG_HepMCValidationHelper
0228 //            std::cout << "removing particle " << **iiso << " from the list of particles to compute isolation because it comes from a tau" << std::endl;
0229 #endif
0230             marked = true;
0231             break;
0232           }
0233         }
0234         if (!marked)
0235           forIsolation.push_back(*iiso);
0236       }
0237       //no compute isolation wrt the status 2 tau direction
0238       double sumIso = 0;
0239       for (iiso = forIsolation.begin(); iiso < forIsolation.end(); ++iiso) {
0240         if (deltaR(taumom, (*iiso)->momentum()) < deltaRcut) {
0241           sumIso += (*iiso)->momentum().perp();
0242         }
0243       }
0244       //if isolated remove the tau daughters from the pruned list
0245       if (sumIso < sumPtCut) {
0246 #ifdef DEBUG_HepMCValidationHelper
0247         std::cout << "particle " << *taus[i] << " is considered isolated, with sumPt " << sumIso << std::endl;
0248 #endif
0249         toRemove.insert(toRemove.end(), taudaughters.begin(), taudaughters.end());
0250       }
0251     }
0252
0253     //now actually remove
0254     pruned.clear();
0255     for (unsigned int i = 0; i < status1.size(); ++i) {
0256       bool marked = false;
0257       std::vector<const HepMC::GenParticle*>::const_iterator iremove;
0258       for (iremove = toRemove.begin(); iremove != toRemove.end(); ++iremove) {
0259         if (status1[i]->barcode() == (*iremove)->barcode()) {
0260           marked = true;
0261           break;
0262         }
0263       }
0264       if (!marked)
0265         pruned.push_back(status1[i]);
0266     }
0267
0268 #ifdef DEBUG_HepMCValidationHelper
0269     std::cout << "list of remaining particles:" << std::endl;
0270     for (unsigned int i = 0; i < pruned.size(); ++i) {
0271       std::cout << *pruned[i] << std::endl;
0272     }
0273 #endif
0274   }
0275
0276   //get all visible status1 particles
0277   void allVisibleParticles(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& visible) {
0278     std::vector<const HepMC::GenParticle*> status1;
0279     visible.clear();
0280     allStatus1(all, status1);
0281     for (unsigned int i = 0; i < status1.size(); ++i) {
0282       if (!isNeutrino(status1[i]))
0283         visible.push_back(status1[i]);
0284     }
0285   }
0286
0287   //compute generated met
0288   TLorentzVector genMet(const HepMC::GenEvent* all, double etamin, double etamax) {
0289     std::vector<const HepMC::GenParticle*> visible;
0290     allVisibleParticles(all, visible);
0291     TLorentzVector momsum(0., 0., 0., 0.);
0292     for (unsigned int i = 0; i < visible.size(); ++i) {
0293       if (visible[i]->momentum().eta() > etamin && visible[i]->momentum().eta() < etamax) {
0294         TLorentzVector mom(visible[i]->momentum().x(),
0295                            visible[i]->momentum().y(),
0296                            visible[i]->momentum().z(),
0297                            visible[i]->momentum().t());
0298         momsum += mom;
0299       }
0300     }
0301     TLorentzVector met(-momsum.Px(), -momsum.Py(), 0., momsum.Pt());
0302     return met;
0303   }
0304 }  // namespace HepMCValidationHelper