EgammaIsolationAlgos/src/EgammaHLTTrackIsolation.cc

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 ntrack = 0;
0153   float ptSum = 0.;
0154   std::vector<float> etaele;
0155   std::vector<float> phiele;
0156
0157   // std::cout << "allEle.size() = " << allEle->size() << std::endl;
0158
0159   // Store ALL electrons eta and phi
0160   for (reco::ElectronCollection::const_iterator iElectron = allEle->begin(); iElectron != allEle->end(); iElectron++) {
0161     reco::TrackRef anothereletrackref = iElectron->track();
0162     etaele.push_back(anothereletrackref->momentum().eta());
0163     phiele.push_back(anothereletrackref->momentum().phi());
0164   }
0165
0166   for (reco::TrackCollection::const_iterator trItr = isoTracks->begin(); trItr != isoTracks->end(); ++trItr) {
0167     GlobalPoint ivtx(trItr->vertex().x(), trItr->vertex().y(), trItr->vertex().z());
0168     reco::Track::Vector ip = trItr->momentum();
0169     GlobalVector imom(ip.x(), ip.y(), ip.z());
0170
0171     float pt = imom.perp();
0172     float dperp = ivtx.perp() - vtx.perp();
0173     float dz = ivtx.z() - vtx.z();
0174     float deta = imom.eta() - mom.eta();
0175     float dphi = imom.phi() - mom.phi();
0176
0177     // Correct dmom_phi's from [-2pi,2pi] to [-pi,pi]
0178     if (dphi > M_PI)
0179       dphi = dphi - 2 * M_PI;
0180     else if (dphi < -M_PI)
0181       dphi = dphi + 2 * M_PI;
0182
0183     float R = sqrt(dphi * dphi + deta * deta);
0184
0185     // Apply boundary cut
0186     bool selected = false;
0187     bool passedconeveto = true;
0188
0189     //hmm how do I figure out if this is barrel or endcap?
0190     //abs(mom.eta())<1.5 is the obvious choice but that will be electron not detector eta for electrons
0191     //well lets leave it as that for now, its what reco does (well with eta=1.479)
0192     double innerStrip = fabs(mom.eta()) < 1.479 ? stripBarrel : stripEndcap;
0193
0194     if (pt > ptMin && R < conesize && fabs(dperp) < rspan && fabs(dz) < zspan && fabs(deta) >= innerStrip)
0195       selected = true;
0196
0197     // Check that NO electron is counted in the isolation
0198     for (unsigned int eleItr = 0; eleItr < etaele.size(); ++eleItr) {
0199       deta = etaele[eleItr] - imom.eta();
0200       dphi = phiele[eleItr] - imom.phi();
0201
0202       // Correct dmom_phi's from [-2pi,2pi] to [-pi,pi]
0203       if (dphi > M_PI)
0204         dphi = dphi - 2 * M_PI;
0205       else if (dphi < -M_PI)
0206         dphi = dphi + 2 * M_PI;
0207
0208       R = sqrt(dphi * dphi + deta * deta);
0209       if (R < vetoConesize)
0210         passedconeveto = false;
0211     }
0212
0213     if (selected && passedconeveto) {
0214       ntrack++;
0215       ptSum += pt;  //to exclude electron tracks
0216     }
0217   }
0218
0219   // ntrack-=1; //to exclude electron track
0220
0221   return (std::pair<int, float>(ntrack, ptSum));
0222 }