Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:07:27

0001 #include "GeneratorInterface/Pythia8Interface/test/analyserhepmc/LeptonAnalyserHepMC.h"
0002 #include "HepMC/GenVertex.h"
0003 #include "HepMC/GenParticle.h"
0004 #include <iterator>
0005 #include <algorithm>
0006 
0007 LeptonAnalyserHepMC::LeptonAnalyserHepMC(double aMaxEta, double aThresholdEt)
0008     : MaxEta(aMaxEta), ThresholdEt(aThresholdEt), RConeIsol(0.3), MaxPtIsol(2.), RIdJet(0.3), EpsIdJet(0.6) {
0009   ;
0010 }
0011 
0012 std::vector<HepMC::GenParticle> LeptonAnalyserHepMC::isolatedLeptons(const HepMC::GenEvent* pEv) {
0013   HepMC::GenEvent::particle_const_iterator part;
0014   HepMC::GenEvent::particle_const_iterator part1;
0015 
0016   std::vector<HepMC::GenParticle> isoleptons;
0017   bool lepton = false;
0018   for (part = pEv->particles_begin(); part != pEv->particles_end(); ++part) {
0019     lepton = false;
0020     if (abs((*part)->pdg_id()) == 11)
0021       lepton = true;
0022     if (abs((*part)->pdg_id()) == 13)
0023       lepton = true;
0024     if (!(*part)->end_vertex() && (*part)->status() == 1 && lepton && fabs((*part)->momentum().eta()) < MaxEta &&
0025         (*part)->momentum().perp() > ThresholdEt) {
0026       double eta0 = (*part)->momentum().eta();
0027       double phi0 = (*part)->momentum().phi();
0028       double pti, dist, etai, phii, dphi;
0029       bool isol = true;
0030       for (part1 = pEv->particles_begin(); part1 != part && part1 != pEv->particles_end(); part1++) {
0031         if (!(*part1)->end_vertex() && (*part1)->status() == 1) {
0032           pti = (*part1)->momentum().perp();
0033           etai = (*part1)->momentum().eta();
0034           phii = (*part1)->momentum().phi();
0035           dphi = phi0 - phii;
0036           if (fabs(phi0 - phii - 6.2832) < fabs(dphi))
0037             dphi = phi0 - phii - 6.2832;
0038           if (fabs(phi0 - phii + 6.2832) < fabs(dphi))
0039             dphi = phi0 - phii + 6.2832;
0040           dist = sqrt((eta0 - etai) * (eta0 - etai) + dphi * dphi);
0041           if (dist < RConeIsol && pti > MaxPtIsol) {
0042             isol = false;
0043             break;
0044           }
0045         }
0046       }
0047       if (isol)
0048         isoleptons.push_back(HepMC::GenParticle(**part));
0049     }
0050   }
0051   return isoleptons;
0052 }
0053 
0054 int LeptonAnalyserHepMC::nIsolatedLeptons(const HepMC::GenEvent* pEv) {
0055   std::vector<HepMC::GenParticle> isoleptons = isolatedLeptons(pEv);
0056   return isoleptons.size();
0057 }
0058 
0059 double LeptonAnalyserHepMC::MinMass(const HepMC::GenEvent* pEv) {
0060   std::vector<HepMC::GenParticle> isoleptons = isolatedLeptons(pEv);
0061   if (isoleptons.size() < 2)
0062     return 0.;
0063   double MinM = 100000.;
0064   std::vector<HepMC::GenParticle>::iterator ipart, ipart1;
0065   for (ipart = isoleptons.begin(); ipart != isoleptons.end(); ipart++) {
0066     for (ipart1 = isoleptons.begin(); ipart1 != isoleptons.end(); ipart1++) {
0067       if (ipart1 == ipart)
0068         continue;
0069       double px = ipart->momentum().px() + ipart1->momentum().px();
0070       double py = ipart->momentum().py() + ipart1->momentum().py();
0071       double pz = ipart->momentum().pz() + ipart1->momentum().pz();
0072       double e = ipart->momentum().e() + ipart1->momentum().e();
0073       double mass = sqrt(e * e - px * px - py * py - pz * pz);
0074       if (mass < MinM)
0075         MinM = mass;
0076     }
0077   }
0078   return MinM;
0079 }
0080 
0081 std::vector<fastjet::PseudoJet> LeptonAnalyserHepMC::removeLeptonsFromJets(std::vector<fastjet::PseudoJet>& jets,
0082                                                                            const HepMC::GenEvent* pEv) {
0083   std::vector<HepMC::GenParticle> isoleptons = isolatedLeptons(pEv);
0084   if (isoleptons.empty())
0085     return jets;
0086   std::vector<fastjet::PseudoJet>::iterator ijet;
0087   std::vector<HepMC::GenParticle>::iterator ipart;
0088   std::vector<fastjet::PseudoJet> newjets;
0089   for (ijet = jets.begin(); ijet != jets.end(); ijet++) {
0090     if (fabs(ijet->rap()) > 5.0)
0091       continue;
0092     bool isLepton = false;
0093     for (ipart = isoleptons.begin(); ipart != isoleptons.end(); ipart++) {
0094       fastjet::PseudoJet fjLept(
0095           ipart->momentum().px(), ipart->momentum().py(), ipart->momentum().pz(), ipart->momentum().e());
0096       //cout << "lepton eta = " << fjLept.rap() << endl;
0097 
0098       if (fjLept.squared_distance(*ijet) < RIdJet * RIdJet &&
0099           fabs(ijet->e() - ipart->momentum().e()) / ijet->e() < EpsIdJet)
0100         isLepton = true;
0101     }
0102 
0103     if (!isLepton)
0104       newjets.push_back(*ijet);
0105   }
0106   return newjets;
0107 }
0108 
0109 #if 0
0110 vector<Jet> LeptonAnalyserHepMC::removeLeptonsFromJets(vector<Jet>& jets,
0111                                                        HepMC::GenEvent* pEv)
0112 {
0113   vector<HepMC::GenParticle> isoleptons = isolatedLeptons(pEv);
0114   if(isoleptons.empty()) return jets;
0115   vector<Jet>::iterator ijet;
0116   vector<HepMC::GenParticle>::iterator ipart;
0117 #if 00
0118   // this code can be used if first argument is vector<Jet>, not vector<Jet>&
0119   for ( ijet = jets.begin(); ijet < jets.end(); ijet++) {
0120     bool bad = false;
0121     for ( ipart=isoleptons.begin(); ipart != isoleptons.end(); ipart++) {
0122       JetableObject jpart(*ipart);
0123       if(ijet->dist(jpart) < RIdJet &&
0124          fabs(ijet->e()-ipart->momentum().e())/ijet->e() < EpsIdJet )
0125                                                              bad = true;
0126     }
0127     if(bad) {ijet = jets.erase(ijet); ijet--;}
0128   }
0129   return jets;
0130 #endif
0131   vector<Jet> newjets;
0132   for ( ijet = jets.begin(); ijet != jets.end(); ijet++) {
0133     bool islepton = false;
0134     for ( ipart=isoleptons.begin(); ipart != isoleptons.end(); ipart++) {
0135       JetableObject jpart(*ipart);
0136       if(ijet->dist(jpart) < RIdJet &&
0137          fabs(ijet->e()-ipart->momentum().e())/ijet->e() < EpsIdJet )
0138                                                        islepton = true;
0139     }
0140     if(!islepton) newjets.push_back(*ijet);
0141   }
0142   return newjets;
0143 }
0144 #endif