File indexing completed on 2024-04-06 12:24:47
0001 #include "RecoEgamma/EgammaElectronAlgos/interface/ConversionFinder.h"
0002 #include "DataFormats/Math/interface/deltaR.h"
0003 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0004
0005 namespace egamma::conv {
0006
0007 namespace {
0008
0009 constexpr float square(float x) { return x * x; };
0010
0011 }
0012
0013
0014 constexpr float maxRefPtErrorForKfConv = 0.05f;
0015 constexpr float maxRelGsfPtErrorForKfConv = 0.25f;
0016 constexpr float maxRelPtDiffForKfConv = 0.2f;
0017 constexpr float maxRelPtDiffForGsfConv = 0.25f;
0018 constexpr float dR2Max = square(0.5f);
0019 constexpr int minNumberOfValidHits = 5;
0020 constexpr float maxDist2Dcot2 = square(0.05f);
0021 constexpr int maxDeltaMissingHits = 2;
0022 constexpr float maxRelGsfPtError = 0.5f;
0023 constexpr int maxDeltaMissingHitsForKFtoKF = 3;
0024 constexpr float maxDistOrCotForKFtoKF = 0.02f;
0025
0026 ConversionInfo getConversionInfo(reco::Track const& el_track, TrackRowView const& track, float bFieldAtOrigin);
0027
0028
0029 std::vector<ConversionInfo> findConversions(const reco::GsfElectronCore& gsfElectron,
0030 TrackTableView ctfTable,
0031 TrackTableView gsfTable,
0032 float bFieldAtOrigin,
0033 float minFracSharedHits) {
0034 using namespace reco;
0035 using namespace std;
0036 using namespace edm;
0037 using namespace edm::soa::col;
0038
0039
0040
0041 const reco::TrackRef eleCtfTk = gsfElectron.ctfTrack();
0042 const reco::GsfTrackRef& eleGsfTk = gsfElectron.gsfTrack();
0043
0044 float eleGsfPt = eleGsfTk->pt();
0045 float eleGsfEta = eleGsfTk->eta();
0046 float eleGsfPhi = eleGsfTk->phi();
0047
0048 const bool useEleCtfTrack = eleCtfTk.isNonnull() && gsfElectron.ctfGsfOverlap() > minFracSharedHits;
0049
0050 std::optional<float> eleCtfPt = std::nullopt;
0051 std::optional<float> eleCtfEta = std::nullopt;
0052 std::optional<float> eleCtfPhi = std::nullopt;
0053
0054
0055
0056 std::optional<int> ctfidx = std::nullopt;
0057 int gsfidx = static_cast<int>(eleGsfTk.key());
0058
0059 if (useEleCtfTrack) {
0060 eleCtfPt = eleCtfTk->pt();
0061 eleCtfEta = eleCtfTk->eta();
0062 eleCtfPhi = eleCtfTk->phi();
0063 ctfidx = static_cast<int>(eleCtfTk.key());
0064 }
0065
0066
0067 std::vector<ConversionInfo> v_candidatePartners;
0068
0069 int ctftk_i = 0;
0070 int gsftk_i = 0;
0071
0072
0073 for (auto ctftkItr = ctfTable.begin(); ctftkItr != ctfTable.end(); ++ctftkItr, ctftk_i++) {
0074 if (useEleCtfTrack && ctftk_i == ctfidx.value())
0075 continue;
0076
0077 auto ctftk = *ctftkItr;
0078
0079
0080 if (ctftk.get<PtError>() > maxRefPtErrorForKfConv * ctftk.get<Pt>() ||
0081 ctftk.get<NumberOfValidHits>() < minNumberOfValidHits)
0082 continue;
0083
0084 if (useEleCtfTrack && std::abs(ctftk.get<Pt>() - eleCtfPt.value()) < maxRelPtDiffForKfConv * eleCtfPt.value())
0085 continue;
0086
0087
0088
0089 if (useEleCtfTrack && eleCtfTk->charge() + ctftk.get<Charge>() == 0 &&
0090 deltaR2(eleCtfEta.value(), eleCtfPhi.value(), ctftk.get<Eta>(), ctftk.get<Phi>()) < dR2Max) {
0091 ConversionInfo convInfo = getConversionInfo(*eleCtfTk, ctftk, bFieldAtOrigin);
0092
0093
0094
0095 int deltaMissingHits = ctftk.get<MissingInnerHits>() - eleCtfTk->missingInnerHits();
0096
0097 v_candidatePartners.push_back(
0098 {convInfo.dist, convInfo.dcot, convInfo.radiusOfConversion, ctftk_i, std::nullopt, deltaMissingHits, 0});
0099
0100 }
0101
0102
0103 if (eleGsfTk->charge() + ctftk.get<Charge>() == 0 &&
0104 deltaR2(eleGsfEta, eleGsfPhi, ctftk.get<Eta>(), ctftk.get<Phi>()) < dR2Max &&
0105 eleGsfTk->ptError() < maxRelGsfPtErrorForKfConv * eleGsfPt) {
0106 int deltaMissingHits = ctftk.get<MissingInnerHits>() - eleGsfTk->missingInnerHits();
0107
0108 ConversionInfo convInfo = getConversionInfo(*eleGsfTk, ctftk, bFieldAtOrigin);
0109
0110 v_candidatePartners.push_back(
0111 {convInfo.dist, convInfo.dcot, convInfo.radiusOfConversion, ctftk_i, std::nullopt, deltaMissingHits, 1});
0112 }
0113
0114 }
0115
0116
0117 for (auto gsftkItr = gsfTable.begin(); gsftkItr != gsfTable.end(); ++gsftkItr, gsftk_i++) {
0118
0119 if (gsfidx == gsftk_i)
0120 continue;
0121
0122 auto gsftk = *gsftkItr;
0123
0124
0125 if (gsftk.get<PtError>() > maxRelGsfPtError * gsftk.get<Pt>() ||
0126 gsftk.get<NumberOfValidHits>() < minNumberOfValidHits)
0127 continue;
0128
0129 if (std::abs(gsftk.get<Pt>() - eleGsfPt) < maxRelPtDiffForGsfConv * eleGsfPt)
0130 continue;
0131
0132
0133
0134
0135 if (useEleCtfTrack && eleCtfTk->charge() + gsftk.get<Charge>() == 0 &&
0136 deltaR2(eleCtfEta.value(), eleCtfPhi.value(), gsftk.get<Eta>(), gsftk.get<Phi>()) < dR2Max) {
0137 int deltaMissingHits = gsftk.get<MissingInnerHits>() - eleCtfTk->missingInnerHits();
0138
0139 ConversionInfo convInfo = getConversionInfo(*eleCtfTk, gsftk, bFieldAtOrigin);
0140
0141 v_candidatePartners.push_back(
0142 {convInfo.dist, convInfo.dcot, convInfo.radiusOfConversion, std::nullopt, gsftk_i, deltaMissingHits, 2});
0143 }
0144
0145
0146 if (eleGsfTk->charge() + gsftk.get<Charge>() == 0 &&
0147 deltaR2(eleGsfEta, eleGsfPhi, gsftk.get<Eta>(), gsftk.get<Phi>()) < dR2Max &&
0148 (eleGsfTk->ptError() < maxRelGsfPtError * eleGsfPt)) {
0149 ConversionInfo convInfo = getConversionInfo(*eleGsfTk, gsftk, bFieldAtOrigin);
0150
0151
0152 int deltaMissingHits = gsftk.get<MissingInnerHits>() - eleGsfTk->missingInnerHits();
0153
0154 v_candidatePartners.push_back(
0155 {convInfo.dist, convInfo.dcot, convInfo.radiusOfConversion, std::nullopt, gsftk_i, deltaMissingHits, 3});
0156 }
0157 }
0158
0159 return v_candidatePartners;
0160 }
0161
0162
0163 ConversionInfo getConversionInfo(reco::Track const& ele, TrackRowView const& track, float bFieldAtOrigin) {
0164 using namespace edm::soa::col;
0165
0166
0167 float rEl = 100.f * ele.pt() / (-0.3f * bFieldAtOrigin * ele.charge());
0168 float xEl = -1.f * (rEl - ele.d0()) * sin(ele.phi());
0169 float yEl = (rEl - ele.d0()) * cos(ele.phi());
0170 rEl = std::abs(rEl);
0171
0172 float rCand = 100.f * track.get<Pt>() / (-0.3f * bFieldAtOrigin * track.get<Charge>());
0173 float xCand = -1.f * (rCand - track.get<D0>()) * sin(track.get<Phi>());
0174 float yCand = (rCand - track.get<D0>()) * cos(track.get<Phi>());
0175 rCand = std::abs(rCand);
0176
0177 float d = sqrt(pow(xEl - xCand, 2) + pow(yEl - yCand, 2));
0178 float dist = d - (rEl + rCand);
0179
0180
0181 float dcot = ele.pz() / ele.pt() - track.get<Pz>() / track.get<Pt>();
0182
0183
0184 float xa1 = xEl + (xCand - xEl) * rEl / d;
0185 float xa2 = xCand + (xEl - xCand) * rCand / d;
0186 float ya1 = yEl + (yCand - yEl) * rEl / d;
0187 float ya2 = yCand + (yEl - yCand) * rCand / d;
0188
0189 float x = .5f * (xa1 + xa2);
0190 float y = .5f * (ya1 + ya2);
0191 float rconv = sqrt(pow(x, 2) + pow(y, 2));
0192
0193
0194
0195
0196 float tempsign = ele.px() * x + ele.py() * y;
0197 tempsign = tempsign / std::abs(tempsign);
0198 rconv = tempsign * rconv;
0199
0200
0201 return {dist, dcot, rconv};
0202 }
0203
0204
0205
0206
0207
0208
0209 ConversionInfo const& arbitrateConversionPartnersbyR(const std::vector<ConversionInfo>& convCandidates) {
0210 ConversionInfo const* closestConversion = &convCandidates.front();
0211
0212 if (convCandidates.size() == 1)
0213 return *closestConversion;
0214
0215 float R = pow(closestConversion->dist, 2) + pow(closestConversion->dcot, 2);
0216
0217 for (auto const& temp : convCandidates) {
0218 float temp_R = pow(temp.dist, 2) + pow(temp.dcot, 2);
0219 if (temp_R < R) {
0220 R = temp_R;
0221 closestConversion = &temp;
0222 }
0223 }
0224
0225 return *closestConversion;
0226 }
0227
0228
0229 ConversionInfo findBestConversionMatch(const std::vector<ConversionInfo>& v_convCandidates) {
0230 using namespace std;
0231
0232 if (v_convCandidates.empty())
0233 return {};
0234
0235 if (v_convCandidates.size() == 1)
0236 return v_convCandidates.at(0);
0237
0238 vector<ConversionInfo> v_0;
0239 vector<ConversionInfo> v_1;
0240 vector<ConversionInfo> v_2;
0241 vector<ConversionInfo> v_3;
0242
0243
0244 for (auto const& temp : v_convCandidates) {
0245 if (temp.radiusOfConversion <= -2)
0246 continue;
0247
0248 if (temp.flag == 0) {
0249 if ((std::abs(temp.dist) < maxDistOrCotForKFtoKF && std::abs(temp.dcot) < maxDistOrCotForKFtoKF &&
0250 temp.deltaMissingHits < maxDeltaMissingHitsForKFtoKF) ||
0251 (pow(temp.dist, 2) + pow(temp.dcot, 2) < maxDist2Dcot2 && temp.deltaMissingHits < maxDeltaMissingHits))
0252 v_0.push_back(temp);
0253 }
0254
0255 if (temp.flag == 1) {
0256 if (square(temp.dist) + square(temp.dcot) < maxDist2Dcot2 && temp.deltaMissingHits < maxDeltaMissingHits)
0257 v_1.push_back(temp);
0258 }
0259 if (temp.flag == 2) {
0260 if (square(temp.dist) + square(temp.dcot) < maxDist2Dcot2 && temp.deltaMissingHits < maxDeltaMissingHits)
0261 v_2.push_back(temp);
0262 }
0263 if (temp.flag == 3) {
0264 if (square(temp.dist) + square(temp.dcot) < maxDist2Dcot2 && temp.deltaMissingHits < maxDeltaMissingHits)
0265 v_3.push_back(temp);
0266 }
0267
0268 }
0269
0270
0271
0272
0273
0274 if (!v_0.empty())
0275 return arbitrateConversionPartnersbyR(v_0);
0276
0277 if (!v_1.empty())
0278 return arbitrateConversionPartnersbyR(v_1);
0279
0280 if (!v_2.empty())
0281 return arbitrateConversionPartnersbyR(v_2);
0282
0283 if (!v_3.empty())
0284 return arbitrateConversionPartnersbyR(v_3);
0285
0286
0287
0288
0289 return arbitrateConversionPartnersbyR(v_convCandidates);
0290 }
0291
0292 }