File indexing completed on 2023-05-08 23:19:12
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
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
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
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
0081
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
0102
0103
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
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
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
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
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
0151
0152 int ntrack = 0;
0153 float ptSum = 0.;
0154 std::vector<float> etaele;
0155 std::vector<float> phiele;
0156
0157
0158
0159
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
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
0186 bool selected = false;
0187 bool passedconeveto = true;
0188
0189
0190
0191
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
0198 for (unsigned int eleItr = 0; eleItr < etaele.size(); ++eleItr) {
0199 deta = etaele[eleItr] - imom.eta();
0200 dphi = phiele[eleItr] - imom.phi();
0201
0202
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;
0216 }
0217 }
0218
0219
0220
0221 return (std::pair<int, float>(ntrack, ptSum));
0222 }