Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:13

0001 // -*- C++ -*-
0002 //
0003 // Package:     EgammaHLTAlgos
0004 // Class  :     EgammaHLTTrackIsolation
0005 //
0006 // Implementation:
0007 //     <Notes on implementation>
0008 //
0009 // Original Author:  Monica Vazquez Acosta
0010 //         Created:  Tue Jun 13 12:17:19 CEST 2006
0011 //
0012 
0013 // system include files
0014 
0015 // user include files
0016 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHLTTrackIsolation.h"
0017 
0018 std::pair<int, float> EgammaHLTTrackIsolation::electronIsolation(const reco::Track* const tr,
0019                                                                  const reco::TrackCollection* isoTracks) {
0020   GlobalPoint vtx(0, 0, tr->vertex().z());
0021   const reco::Track::Vector& p = tr->momentum();
0022   GlobalVector mom(p.x(), p.y(), p.z());
0023   return findIsoTracks(mom, vtx, isoTracks, true);
0024 }
0025 
0026 std::pair<int, float> EgammaHLTTrackIsolation::electronIsolation(const reco::Track* const tr,
0027                                                                  const reco::TrackCollection* isoTracks,
0028                                                                  GlobalPoint zvtx) {
0029   // Just to insure consistency with no-vertex-code
0030   GlobalPoint vtx(0, 0, zvtx.z());
0031   const reco::Track::Vector& p = tr->momentum();
0032   GlobalVector mom(p.x(), p.y(), p.z());
0033   return findIsoTracks(mom, vtx, isoTracks, true);
0034 }
0035 
0036 std::pair<int, float> EgammaHLTTrackIsolation::electronIsolation(const reco::Track* const tr,
0037                                                                  const reco::ElectronCollection* allEle,
0038                                                                  const reco::TrackCollection* isoTracks) {
0039   GlobalPoint vtx(0, 0, tr->vertex().z());
0040   const reco::Track::Vector& p = tr->momentum();
0041   GlobalVector mom(p.x(), p.y(), p.z());
0042   return findIsoTracksWithoutEle(mom, vtx, allEle, isoTracks);
0043 }
0044 
0045 std::pair<int, float> EgammaHLTTrackIsolation::photonIsolation(const reco::RecoCandidate* const recocandidate,
0046                                                                const reco::TrackCollection* isoTracks,
0047                                                                bool useVertex) {
0048   if (useVertex) {
0049     GlobalPoint vtx(0, 0, recocandidate->vertex().z());
0050     return photonIsolation(recocandidate, isoTracks, vtx);
0051   } else {
0052     reco::RecoCandidate::Point pos = recocandidate->superCluster()->position();
0053     GlobalVector mom(pos.x(), pos.y(), pos.z());
0054     return findIsoTracks(mom, GlobalPoint(), isoTracks, false, false);
0055   }
0056 }
0057 
0058 std::pair<int, float> EgammaHLTTrackIsolation::photonIsolation(const reco::RecoCandidate* const recocandidate,
0059                                                                const reco::TrackCollection* isoTracks,
0060                                                                GlobalPoint zvtx) {
0061   // to insure consistency with no-free-vertex-code
0062   GlobalPoint vtx(0, 0, zvtx.z());
0063 
0064   reco::RecoCandidate::Point pos = recocandidate->superCluster()->position();
0065   GlobalVector mom(pos.x() - vtx.x(), pos.y() - vtx.y(), pos.z() - vtx.z());
0066 
0067   return findIsoTracks(mom, vtx, isoTracks, false);
0068 }
0069 
0070 std::pair<int, float> EgammaHLTTrackIsolation::photonIsolation(const reco::RecoCandidate* const recocandidate,
0071                                                                const reco::ElectronCollection* allEle,
0072                                                                const reco::TrackCollection* isoTracks) {
0073   reco::RecoCandidate::Point pos = recocandidate->superCluster()->position();
0074   GlobalVector mom(pos.x(), pos.y(), pos.z());
0075   return findIsoTracksWithoutEle(mom, GlobalPoint(), allEle, isoTracks);
0076 }
0077 
0078 std::pair<int, float> EgammaHLTTrackIsolation::findIsoTracks(
0079     GlobalVector mom, GlobalPoint vtx, const reco::TrackCollection* isoTracks, bool isElectron, bool useVertex) {
0080   // Check that reconstructed tracks fit within cone boundaries,
0081   // (Note: tracks will not always stay within boundaries)
0082   int ntrack = 0;
0083   float ptSum = 0.;
0084 
0085   for (reco::TrackCollection::const_iterator trItr = isoTracks->begin(); trItr != isoTracks->end(); ++trItr) {
0086     GlobalPoint ivtx(trItr->vertex().x(), trItr->vertex().y(), trItr->vertex().z());
0087     reco::Track::Vector ip = trItr->momentum();
0088     GlobalVector imom(ip.x(), ip.y(), ip.z());
0089 
0090     float pt = imom.perp();
0091     float dperp = 0.;
0092     float dz = 0.;
0093     float deta = 0.;
0094     float dphi = 0.;
0095     if (useVertex) {
0096       dperp = ivtx.perp() - vtx.perp();
0097       dz = ivtx.z() - vtx.z();
0098       deta = imom.eta() - mom.eta();
0099       dphi = imom.phi() - mom.phi();
0100     } else {
0101       //in case of unkown photon vertex, modify direction of photon to point from
0102       //current track vertex to sc instead of from (0.,0.,0.) to sc.  In this
0103       //way, all tracks are considered based on direction alone.
0104       GlobalVector mom_temp = mom - GlobalVector(ivtx.x(), ivtx.y(), ivtx.z());
0105       deta = imom.eta() - mom_temp.eta();
0106       dphi = imom.phi() - mom_temp.phi();
0107     }
0108     // Correct dmom_phi's from [-2pi,2pi] to [-pi,pi]
0109     if (dphi > M_PI)
0110       dphi = dphi - 2 * M_PI;
0111     else if (dphi < -M_PI)
0112       dphi = dphi + 2 * M_PI;
0113 
0114     float R = sqrt(dphi * dphi + deta * deta);
0115 
0116     // Apply boundary cut
0117     // bool selected=false;
0118 
0119     // if (pt > ptMin && R < conesize &&
0120     //  fabs(dperp) < rspan && fabs(dz) < zspan) selected=true;
0121 
0122     // if (selected) {
0123     //  ntrack++;
0124     //  if (!isElectron || R > vetoConesize) ptSum+=pt; //to exclude electron track
0125     // }
0126     // float theVetoVar = R;
0127     // if (isElectron) theVetoVar = R;
0128 
0129     //hmm how do I figure out if this is barrel or endcap?
0130     //abs(mom.eta())<1.5 is the obvious choice but that will be electron not detector eta for electrons
0131     //well lets leave it as that for now, its what reco does (well with eta=1.479)
0132     double innerStrip = fabs(mom.eta()) < 1.479 ? stripBarrel : stripEndcap;
0133 
0134     if (pt > ptMin && R < conesize && R > vetoConesize && fabs(dperp) < rspan && fabs(dz) < zspan &&
0135         fabs(deta) >= innerStrip) {
0136       ntrack++;
0137       ptSum += pt;
0138     }
0139   }
0140 
0141   // if (isElectron) ntrack-=1; //to exclude electron track
0142 
0143   return (std::pair<int, float>(ntrack, ptSum));
0144 }
0145 
0146 std::pair<int, float> EgammaHLTTrackIsolation::findIsoTracksWithoutEle(GlobalVector mom,
0147                                                                        GlobalPoint vtx,
0148                                                                        const reco::ElectronCollection* allEle,
0149                                                                        const reco::TrackCollection* isoTracks) {
0150   // Check that reconstructed tracks fit within cone boundaries,
0151   // (Note: tracks will not always stay within boundaries)
0152   int iele = 0;
0153   int ntrack = 0;
0154   float ptSum = 0.;
0155   std::vector<float> etaele;
0156   std::vector<float> phiele;
0157 
0158   // std::cout << "allEle.size() = " << allEle->size() << std::endl;
0159 
0160   // Store ALL electrons eta and phi
0161   for (reco::ElectronCollection::const_iterator iElectron = allEle->begin(); iElectron != allEle->end(); iElectron++) {
0162     iele++;
0163     reco::TrackRef anothereletrackref = iElectron->track();
0164     etaele.push_back(anothereletrackref->momentum().eta());
0165     phiele.push_back(anothereletrackref->momentum().phi());
0166     // std::cout << "Electron " << iele << ": phi = " << anothereletrackref->momentum().phi() << ", eta = " << anothereletrackref->momentum().eta() << std::endl;
0167   }
0168 
0169   for (reco::TrackCollection::const_iterator trItr = isoTracks->begin(); trItr != isoTracks->end(); ++trItr) {
0170     GlobalPoint ivtx(trItr->vertex().x(), trItr->vertex().y(), trItr->vertex().z());
0171     reco::Track::Vector ip = trItr->momentum();
0172     GlobalVector imom(ip.x(), ip.y(), ip.z());
0173 
0174     float pt = imom.perp();
0175     float dperp = ivtx.perp() - vtx.perp();
0176     float dz = ivtx.z() - vtx.z();
0177     float deta = imom.eta() - mom.eta();
0178     float dphi = imom.phi() - mom.phi();
0179 
0180     // Correct dmom_phi's from [-2pi,2pi] to [-pi,pi]
0181     if (dphi > M_PI)
0182       dphi = dphi - 2 * M_PI;
0183     else if (dphi < -M_PI)
0184       dphi = dphi + 2 * M_PI;
0185 
0186     float R = sqrt(dphi * dphi + deta * deta);
0187 
0188     // Apply boundary cut
0189     bool selected = false;
0190     bool passedconeveto = true;
0191 
0192     //hmm how do I figure out if this is barrel or endcap?
0193     //abs(mom.eta())<1.5 is the obvious choice but that will be electron not detector eta for electrons
0194     //well lets leave it as that for now, its what reco does (well with eta=1.479)
0195     double innerStrip = fabs(mom.eta()) < 1.479 ? stripBarrel : stripEndcap;
0196 
0197     if (pt > ptMin && R < conesize && fabs(dperp) < rspan && fabs(dz) < zspan && fabs(deta) >= innerStrip)
0198       selected = true;
0199 
0200     // Check that NO electron is counted in the isolation
0201     for (unsigned int eleItr = 0; eleItr < etaele.size(); ++eleItr) {
0202       deta = etaele[eleItr] - imom.eta();
0203       dphi = phiele[eleItr] - imom.phi();
0204 
0205       // Correct dmom_phi's from [-2pi,2pi] to [-pi,pi]
0206       if (dphi > M_PI)
0207         dphi = dphi - 2 * M_PI;
0208       else if (dphi < -M_PI)
0209         dphi = dphi + 2 * M_PI;
0210 
0211       R = sqrt(dphi * dphi + deta * deta);
0212       if (R < vetoConesize)
0213         passedconeveto = false;
0214     }
0215 
0216     if (selected && passedconeveto) {
0217       ntrack++;
0218       ptSum += pt;  //to exclude electron tracks
0219     }
0220   }
0221 
0222   // ntrack-=1; //to exclude electron track
0223 
0224   return (std::pair<int, float>(ntrack, ptSum));
0225 }