Egamma/src/ConversionTools.cc

0001 #include "CommonTools/Egamma/interface/ConversionTools.h"
0002 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
0003 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0004 #include "DataFormats/TrackReco/interface/Track.h"
0005 #include "DataFormats/Common/interface/RefToPtr.h"
0006 #include "DataFormats/Math/interface/deltaPhi.h"
0007 #include "DataFormats/Math/interface/deltaR.h"
0008
0009 #include <TMath.h>
0010
0011 using namespace edm;
0012 using namespace reco;
0013
0014 //--------------------------------------------------------------------------------------------------
0015 bool ConversionTools::isGoodConversion(const Conversion &conv,
0016                                        const math::XYZPoint &beamspot,
0017                                        float lxyMin,
0018                                        float probMin,
0019                                        unsigned int nHitsBeforeVtxMax) {
0020   //Check if a given conversion candidate passes the conversion selection cuts
0021
0022   const reco::Vertex &vtx = conv.conversionVertex();
0023
0024   //vertex validity
0025   if (!vtx.isValid())
0026     return false;
0027
0028   //fit probability
0029   if (TMath::Prob(vtx.chi2(), vtx.ndof()) < probMin)
0030     return false;
0031
0032   //compute transverse decay length
0033   math::XYZVector mom(conv.refittedPairMomentum());
0034   double dbsx = vtx.x() - beamspot.x();
0035   double dbsy = vtx.y() - beamspot.y();
0036   double lxy = (mom.x() * dbsx + mom.y() * dbsy) / mom.rho();
0037
0038   //transverse decay length
0039   if (lxy < lxyMin)
0040     return false;
0041
0042   //loop through daughters to check nhitsbeforevtx
0043   for (std::vector<uint8_t>::const_iterator it = conv.nHitsBeforeVtx().begin(); it != conv.nHitsBeforeVtx().end();
0044        ++it) {
0045     if ((*it) > nHitsBeforeVtxMax)
0046       return false;
0047   }
0048
0049   return true;
0050 }
0051
0052 //--------------------------------------------------------------------------------------------------
0053 bool ConversionTools::matchesConversion(const reco::GsfElectron &ele,
0054                                         const reco::Conversion &conv,
0055                                         bool allowCkfMatch,
0056                                         bool allowAmbiguousGsfMatch) {
0057   //check if a given GsfElectron matches a given conversion (no quality cuts applied)
0058   //matching is always attempted through the gsf track ref, and optionally attempted through the
0059   //closest ctf track ref
0060
0061   const std::vector<edm::RefToBase<reco::Track> > &convTracks = conv.tracks();
0062   for (std::vector<edm::RefToBase<reco::Track> >::const_iterator it = convTracks.begin(); it != convTracks.end();
0063        ++it) {
0064     if (ele.reco::GsfElectron::gsfTrack().isNonnull() && ele.reco::GsfElectron::gsfTrack().id() == it->id() &&
0065         ele.reco::GsfElectron::gsfTrack().key() == it->key())
0066       return true;
0067     else if (allowCkfMatch && ele.reco::GsfElectron::closestCtfTrackRef().isNonnull() &&
0068              ele.reco::GsfElectron::closestCtfTrackRef().id() == it->id() &&
0069              ele.reco::GsfElectron::closestCtfTrackRef().key() == it->key())
0070       return true;
0071     if (allowAmbiguousGsfMatch) {
0072       for (auto const &tk : ele.ambiguousGsfTracks()) {
0073         if (tk.isNonnull() && tk.id() == it->id() && tk.key() == it->key())
0074           return true;
0075       }
0076     }
0077   }
0078
0079   return false;
0080 }
0081
0082 //--------------------------------------------------------------------------------------------------
0083 bool ConversionTools::matchesConversion(const reco::GsfElectronCore &eleCore,
0084                                         const reco::Conversion &conv,
0085                                         bool allowCkfMatch) {
0086   //check if a given GsfElectronCore matches a given conversion (no quality cuts applied)
0087   //matching is always attempted through the gsf track ref, and optionally attempted through the
0088   //closest ctf track ref
0089
0090   for (const auto &trkRef : conv.tracks()) {
0091     if (eleCore.gsfTrack().isNonnull() && eleCore.gsfTrack().id() == trkRef.id() &&
0092         eleCore.gsfTrack().key() == trkRef.key())
0093       return true;
0094     else if (allowCkfMatch && eleCore.ctfTrack().isNonnull() && eleCore.ctfTrack().id() == trkRef.id() &&
0095              eleCore.ctfTrack().key() == trkRef.key())
0096       return true;
0097   }
0098
0099   return false;
0100 }
0101
0102 //--------------------------------------------------------------------------------------------------
0103 bool ConversionTools::matchesConversion(
0104     const reco::SuperCluster &sc, const reco::Conversion &conv, float dRMax, float dEtaMax, float dPhiMax) {
0105   //check if a given SuperCluster matches a given conversion (no quality cuts applied)
0106   //matching is geometric between conversion momentum and vector joining conversion vertex
0107   //to supercluster position
0108
0109   math::XYZVector mom(conv.refittedPairMomentum());
0110
0111   const math::XYZPoint &scpos(sc.position());
0112   math::XYZPoint cvtx(conv.conversionVertex().position());
0113
0114   math::XYZVector cscvector = scpos - cvtx;
0115   float dR = reco::deltaR(mom, cscvector);
0116   float dEta = mom.eta() - cscvector.eta();
0117   float dPhi = reco::deltaPhi(mom.phi(), cscvector.phi());
0118
0119   if (dR > dRMax)
0120     return false;
0121   if (dEta > dEtaMax)
0122     return false;
0123   if (dPhi > dPhiMax)
0124     return false;
0125
0126   return true;
0127 }
0128
0129 //--------------------------------------------------------------------------------------------------
0130 bool ConversionTools::matchesConversion(const edm::RefToBase<reco::Track> &trk, const reco::Conversion &conv) {
0131   //check if given track matches given conversion (matching by ref)
0132
0133   if (trk.isNull())
0134     return false;
0135
0136   const std::vector<edm::RefToBase<reco::Track> > &convTracks = conv.tracks();
0137   for (std::vector<edm::RefToBase<reco::Track> >::const_iterator it = convTracks.begin(); it != convTracks.end();
0138        ++it) {
0139     if (trk.id() == it->id() && trk.key() == it->key())
0140       return true;
0141   }
0142
0143   return false;
0144 }
0145
0146 //--------------------------------------------------------------------------------------------------
0147 bool ConversionTools::matchesConversion(const reco::TrackRef &trk, const reco::Conversion &conv) {
0148   //check if given track matches given conversion (matching by ref)
0149
0150   if (trk.isNull())
0151     return false;
0152
0153   const std::vector<edm::RefToBase<reco::Track> > &convTracks = conv.tracks();
0154   for (std::vector<edm::RefToBase<reco::Track> >::const_iterator it = convTracks.begin(); it != convTracks.end();
0155        ++it) {
0156     if (trk.id() == it->id() && trk.key() == it->key())
0157       return true;
0158   }
0159
0160   return false;
0161 }
0162
0163 //--------------------------------------------------------------------------------------------------
0164 bool ConversionTools::matchesConversion(const reco::GsfTrackRef &trk, const reco::Conversion &conv) {
0165   //check if given track matches given conversion (matching by ref)
0166
0167   if (trk.isNull())
0168     return false;
0169
0170   const std::vector<edm::RefToBase<reco::Track> > &convTracks = conv.tracks();
0171   for (std::vector<edm::RefToBase<reco::Track> >::const_iterator it = convTracks.begin(); it != convTracks.end();
0172        ++it) {
0173     if (trk.id() == it->id() && trk.key() == it->key())
0174       return true;
0175   }
0176
0177   return false;
0178 }
0179
0180 //--------------------------------------------------------------------------------------------------
0181 bool ConversionTools::hasMatchedConversion(const reco::GsfElectron &ele,
0182                                            const reco::ConversionCollection &convCol,
0183                                            const math::XYZPoint &beamspot,
0184                                            bool allowCkfMatch,
0185                                            float lxyMin,
0186                                            float probMin,
0187                                            unsigned int nHitsBeforeVtxMax) {
0188   //check if a given electron candidate matches to at least one conversion candidate in the
0189   //collection which also passes the selection cuts, optionally match with the closestckf track in
0190   //in addition to just the gsf track (enabled in default arguments)
0191
0192   for (auto const &it : convCol) {
0193     if (!matchesConversion(ele, it, allowCkfMatch))
0194       continue;
0195     if (!isGoodConversion(it, beamspot, lxyMin, probMin, nHitsBeforeVtxMax))
0196       continue;
0197
0198     return true;
0199   }
0200
0201   return false;
0202 }
0203
0204 //--------------------------------------------------------------------------------------------------
0205 bool ConversionTools::hasMatchedConversion(const reco::TrackRef &trk,
0206                                            const reco::ConversionCollection &convCol,
0207                                            const math::XYZPoint &beamspot,
0208                                            float lxyMin,
0209                                            float probMin,
0210                                            unsigned int nHitsBeforeVtxMax) {
0211   //check if a given track matches to at least one conversion candidate in the
0212   //collection which also passes the selection cuts
0213
0214   if (trk.isNull())
0215     return false;
0216
0217   for (auto const &it : convCol) {
0218     if (!matchesConversion(trk, it))
0219       continue;
0220     if (!isGoodConversion(it, beamspot, lxyMin, probMin, nHitsBeforeVtxMax))
0221       continue;
0222
0223     return true;
0224   }
0225
0226   return false;
0227 }
0228
0229 //--------------------------------------------------------------------------------------------------
0230 bool ConversionTools::hasMatchedConversion(const reco::SuperCluster &sc,
0231                                            const reco::ConversionCollection &convCol,
0232                                            const math::XYZPoint &beamspot,
0233                                            float dRMax,
0234                                            float dEtaMax,
0235                                            float dPhiMax,
0236                                            float lxyMin,
0237                                            float probMin,
0238                                            unsigned int nHitsBeforeVtxMax) {
0239   //check if a given SuperCluster matches to at least one conversion candidate in the
0240   //collection which also passes the selection cuts
0241
0242   for (auto const &it : convCol) {
0243     if (!matchesConversion(sc, it))
0244       continue;
0245     if (!isGoodConversion(it, beamspot, lxyMin, probMin, nHitsBeforeVtxMax))
0246       continue;
0247
0248     return true;
0249   }
0250
0251   return false;
0252 }
0253
0254 //--------------------------------------------------------------------------------------------------
0255 reco::Conversion const *ConversionTools::matchedConversion(const reco::GsfElectron &ele,
0256                                                            const reco::ConversionCollection &convCol,
0257                                                            const math::XYZPoint &beamspot,
0258                                                            bool allowCkfMatch,
0259                                                            float lxyMin,
0260                                                            float probMin,
0261                                                            unsigned int nHitsBeforeVtxMax) {
0262   //check if a given electron candidate matches to at least one conversion candidate in the
0263   //collection which also passes the selection cuts, optionally match with the closestckf track in
0264   //in addition to just the gsf track (enabled in default arguments)
0265   //If multiple conversions are found, returned reference corresponds to minimum
0266   //conversion radius
0267
0268   reco::Conversion const *match = nullptr;
0269
0270   double minRho = 999.;
0271   for (auto const &it : convCol) {
0272     float rho = it.conversionVertex().position().rho();
0273     if (rho > minRho)
0274       continue;
0275     if (!matchesConversion(ele, it, allowCkfMatch))
0276       continue;
0277     if (!isGoodConversion(it, beamspot, lxyMin, probMin, nHitsBeforeVtxMax))
0278       continue;
0279
0280     minRho = rho;
0281     match = &it;
0282   }
0283
0284   return match;
0285 }
0286
0287 //--------------------------------------------------------------------------------------------------
0288 reco::Conversion const *ConversionTools::matchedConversion(const reco::GsfElectronCore &eleCore,
0289                                                            const reco::ConversionCollection &convCol,
0290                                                            const math::XYZPoint &beamspot,
0291                                                            bool allowCkfMatch,
0292                                                            float lxyMin,
0293                                                            float probMin,
0294                                                            unsigned int nHitsBeforeVtxMax) {
0295   //check if a given electron candidate matches to at least one conversion candidate in the
0296   //collection which also passes the selection cuts, optionally match with the closestckf track in
0297   //in addition to just the gsf track (enabled in default arguments)
0298   //If multiple conversions are found, returned reference corresponds to minimum
0299   //conversion radius
0300
0301   reco::Conversion const *match = nullptr;
0302
0303   double minRho = 999.;
0304   for (auto const &it : convCol) {
0305     float rho = it.conversionVertex().position().rho();
0306     if (rho > minRho)
0307       continue;
0308     if (!matchesConversion(eleCore, it, allowCkfMatch))
0309       continue;
0310     if (!isGoodConversion(it, beamspot, lxyMin, probMin, nHitsBeforeVtxMax))
0311       continue;
0312
0313     minRho = rho;
0314     match = &it;
0315   }
0316
0317   return match;
0318 }
0319
0320 //--------------------------------------------------------------------------------------------------
0321 reco::Conversion const *ConversionTools::matchedConversion(const reco::TrackRef &trk,
0322                                                            const reco::ConversionCollection &convCol,
0323                                                            const math::XYZPoint &beamspot,
0324                                                            float lxyMin,
0325                                                            float probMin,
0326                                                            unsigned int nHitsBeforeVtxMax) {
0327   //check if a given track matches to at least one conversion candidate in the
0328   //collection which also passes the selection cuts
0329   //If multiple conversions are found, returned reference corresponds to minimum
0330   //conversion radius
0331
0332   reco::Conversion const *match = nullptr;
0333
0334   if (trk.isNull())
0335     return match;
0336
0337   double minRho = 999.;
0338   for (auto const &it : convCol) {
0339     float rho = it.conversionVertex().position().rho();
0340     if (rho > minRho)
0341       continue;
0342     if (!matchesConversion(trk, it))
0343       continue;
0344     if (!isGoodConversion(it, beamspot, lxyMin, probMin, nHitsBeforeVtxMax))
0345       continue;
0346
0347     minRho = rho;
0348     match = &it;
0349   }
0350
0351   return match;
0352 }
0353
0354 //--------------------------------------------------------------------------------------------------
0355 reco::Conversion const *ConversionTools::matchedConversion(const reco::SuperCluster &sc,
0356                                                            const reco::ConversionCollection &convCol,
0357                                                            const math::XYZPoint &beamspot,
0358                                                            float dRMax,
0359                                                            float dEtaMax,
0360                                                            float dPhiMax,
0361                                                            float lxyMin,
0362                                                            float probMin,
0363                                                            unsigned int nHitsBeforeVtxMax) {
0364   //check if a given SuperCluster matches to at least one conversion candidate in the
0365   //collection which also passes the selection cuts
0366   //If multiple conversions are found, returned reference corresponds to minimum
0367   //conversion radius
0368
0369   reco::Conversion const *match = nullptr;
0370
0371   double minRho = 999.;
0372   for (auto const &it : convCol) {
0373     float rho = it.conversionVertex().position().rho();
0374     if (rho > minRho)
0375       continue;
0376     if (!matchesConversion(sc, it, dRMax, dEtaMax, dPhiMax))
0377       continue;
0378     if (!isGoodConversion(it, beamspot, lxyMin, probMin, nHitsBeforeVtxMax))
0379       continue;
0380
0381     minRho = rho;
0382     match = &it;
0383   }
0384
0385   return match;
0386 }
0387
0388 //--------------------------------------------------------------------------------------------------
0389 bool ConversionTools::hasMatchedPromptElectron(const reco::SuperClusterRef &sc,
0390                                                const reco::GsfElectronCollection &eleCol,
0391                                                const reco::ConversionCollection &convCol,
0392                                                const math::XYZPoint &beamspot,
0393                                                bool allowCkfMatch,
0394                                                float lxyMin,
0395                                                float probMin,
0396                                                unsigned int nHitsBeforeVtxMax) {
0397   return !(matchedPromptElectron(sc, eleCol, convCol, beamspot, allowCkfMatch, lxyMin, probMin, nHitsBeforeVtxMax) ==
0398            nullptr);
0399 }
0400
0401 //--------------------------------------------------------------------------------------------------
0402 reco::GsfElectron const *ConversionTools::matchedPromptElectron(const reco::SuperClusterRef &sc,
0403                                                                 const reco::GsfElectronCollection &eleCol,
0404                                                                 const reco::ConversionCollection &convCol,
0405                                                                 const math::XYZPoint &beamspot,
0406                                                                 bool allowCkfMatch,
0407                                                                 float lxyMin,
0408                                                                 float probMin,
0409                                                                 unsigned int nHitsBeforeVtxMax) {
0410   //check if a given SuperCluster matches to at least one GsfElectron having zero expected inner hits
0411   //and not matching any conversion in the collection passing the quality cuts
0412
0413   reco::GsfElectron const *match = nullptr;
0414
0415   if (sc.isNull())
0416     return match;
0417
0418   for (auto const &it : eleCol) {
0419     //match electron to supercluster
0420     if (it.superCluster() != sc)
0421       continue;
0422
0423     //check expected inner hits
0424     if (it.gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) > 0)
0425       continue;
0426
0427     //check if electron is matching to a conversion
0428     if (hasMatchedConversion(it, convCol, beamspot, allowCkfMatch, lxyMin, probMin, nHitsBeforeVtxMax))
0429       continue;
0430
0431     match = &it;
0432   }
0433
0434   return match;
0435 }
0436
0437 //--------------------------------------------------------------------------------------------------
0438 float ConversionTools::getVtxFitProb(const reco::Conversion *conv) {
0439   if (conv != nullptr) {
0440     const reco::Vertex &vtx = conv->conversionVertex();
0441     if (vtx.isValid()) {
0442       return TMath::Prob(vtx.chi2(), vtx.ndof());
0443     }
0444   }
0445   return -1;
0446 }