Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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   }  // namespace
0012 
0013   // configuration parameters
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     //get the references to the gsf and ctf tracks that are made
0040     //by the electron
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     //the electron's CTF track must share at least 45% of the inner hits
0055     //with the electron's GSF track
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     //these vectors are for those candidate partner tracks that pass our cuts
0067     std::vector<ConversionInfo> v_candidatePartners;
0068     //track indices required to make references
0069     int ctftk_i = 0;
0070     int gsftk_i = 0;
0071 
0072     //loop over the CTF tracks and try to find the partner track
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       //apply quality cuts to remove bad tracks
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       //use the electron's CTF track, if not null, to search for the partner track
0088       //look only in a cone of 0.5 to save time, and require that the track is opp. sign
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         //need to add the track reference information for completeness
0094         //because the overloaded fnc above does not make a trackRef
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       }  //using the electron's CTF track
0101 
0102       //now we check using the electron's gsf track
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       }  //using the electron's GSF track
0113 
0114     }  //loop over the CTF track collection
0115 
0116     //------------------------------------------------------ Loop over GSF collection ----------------------------------//
0117     for (auto gsftkItr = gsfTable.begin(); gsftkItr != gsfTable.end(); ++gsftkItr, gsftk_i++) {
0118       //reject the electron's own gsfTrack
0119       if (gsfidx == gsftk_i)
0120         continue;
0121 
0122       auto gsftk = *gsftkItr;
0123 
0124       //apply quality cuts to remove bad tracks
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       //try using the electron's CTF track first if it exists
0133       //look only in a cone of 0.5 around the electron's track
0134       //require opposite sign
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         //fill the Ref info
0141         v_candidatePartners.push_back(
0142             {convInfo.dist, convInfo.dcot, convInfo.radiusOfConversion, std::nullopt, gsftk_i, deltaMissingHits, 2});
0143       }
0144 
0145       //use the electron's gsf track
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         //fill the Ref info
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     }  //loop over the gsf track collection
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     //now calculate the conversion related information
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     // this is equivalent to `1/tan(theta_1) - 1/tan(theta_2)` but requires less trigonometry
0181     float dcot = ele.pz() / ele.pt() - track.get<Pz>() / track.get<Pt>();
0182 
0183     //get the point of conversion
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     // The z-position of the conversion is unused, but here is how it could be computed if needed:
0193     // float z = ele.dz() + rEl * ele.pz() * std::acos(1 - pow(rconv, 2) / (2. * pow(rEl, 2))) / ele.pt();
0194 
0195     //now assign a sign to the radius of conversion
0196     float tempsign = ele.px() * x + ele.py() * y;
0197     tempsign = tempsign / std::abs(tempsign);
0198     rconv = tempsign * rconv;
0199 
0200     //return an instance of ConversionInfo, but with a NULL track refs
0201     return {dist, dcot, rconv};
0202   }
0203 
0204   //------------------------------------------------------------------------------------
0205 
0206   //takes in a vector of candidate conversion partners
0207   //and arbitrates between them returning the one with the
0208   //smallest R=sqrt(dist*dist + dcot*dcot)
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     //loop over the candidates
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     }  //candidate conversion loop
0269 
0270     //now do some arbitration
0271 
0272     //give preference to conversion partners found in the CTF collection
0273     //using the electron's CTF track
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     //if we get here, we didn't find a candidate conversion partner that
0287     //satisfied even the loose selections
0288     //return the the closest partner by R
0289     return arbitrateConversionPartnersbyR(v_convCandidates);
0290   }
0291 
0292 }  // namespace egamma::conv