File indexing completed on 2024-04-06 12:13:58
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
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
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
