EgammaIsolationAlgos/src/EleTkIsolFromCands.cc

0001 #include "RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h"
0002 #include "DataFormats/TrackReco/interface/Track.h"
0003 #include "DataFormats/Math/interface/deltaR.h"
0004
0005 EleTkIsolFromCands::TrkCuts::TrkCuts(const edm::ParameterSet& para) {
0006   minPt = para.getParameter<double>("minPt");
0007   auto sq = [](double val) { return val * val; };
0008   minDR2 = sq(para.getParameter<double>("minDR"));
0009   maxDR2 = sq(para.getParameter<double>("maxDR"));
0010   minDEta = para.getParameter<double>("minDEta");
0011   maxDZ = para.getParameter<double>("maxDZ");
0012   minHits = para.getParameter<int>("minHits");
0013   minPixelHits = para.getParameter<int>("minPixelHits");
0014   maxDPtPt = para.getParameter<double>("maxDPtPt");
0015
0016   auto qualNames = para.getParameter<std::vector<std::string> >("allowedQualities");
0017   auto algoNames = para.getParameter<std::vector<std::string> >("algosToReject");
0018
0019   for (auto& qualName : qualNames) {
0020     allowedQualities.push_back(reco::TrackBase::qualityByName(qualName));
0021   }
0022   for (auto& algoName : algoNames) {
0023     algosToReject.push_back(reco::TrackBase::algoByName(algoName));
0024   }
0025   std::sort(algosToReject.begin(), algosToReject.end());
0026 }
0027
0028 edm::ParameterSetDescription EleTkIsolFromCands::TrkCuts::pSetDescript() {
0029   edm::ParameterSetDescription desc;
0030   desc.add<double>("minPt", 1.0);
0031   desc.add<double>("maxDR", 0.3);
0032   desc.add<double>("minDR", 0.000);
0033   desc.add<double>("minDEta", 0.005);
0034   desc.add<double>("maxDZ", 0.1);
0035   desc.add<double>("maxDPtPt", -1);
0036   desc.add<int>("minHits", 8);
0037   desc.add<int>("minPixelHits", 1);
0038   desc.add<std::vector<std::string> >("allowedQualities");
0039   desc.add<std::vector<std::string> >("algosToReject");
0040   return desc;
0041 }
0042
0043 edm::ParameterSetDescription EleTkIsolFromCands::pSetDescript() {
0044   edm::ParameterSetDescription desc;
0045   desc.add("barrelCuts", TrkCuts::pSetDescript());
0046   desc.add("endcapCuts", TrkCuts::pSetDescript());
0047   return desc;
0048 }
0049
0050 EleTkIsolFromCands::TrackTable EleTkIsolFromCands::preselectTracks(reco::TrackCollection const& tracks,
0051                                                                    TrkCuts const& cuts) {
0052   std::vector<float> pt;
0053   std::vector<float> eta;
0054   std::vector<float> phi;
0055   std::vector<float> vz;
0056   pt.reserve(tracks.size());
0057   eta.reserve(tracks.size());
0058   phi.reserve(tracks.size());
0059   vz.reserve(tracks.size());
0060
0061   for (auto const& trk : tracks) {
0062     auto trkPt = trk.pt();
0063     if (passTrackPreselection(trk, trkPt, cuts)) {
0064       pt.emplace_back(trkPt);
0065       eta.emplace_back(trk.eta());
0066       phi.emplace_back(trk.phi());
0067       vz.emplace_back(trk.vz());
0068     }
0069   }
0070
0071   return {pt, std::move(eta), std::move(phi), std::move(vz)};
0072 }
0073
0074 EleTkIsolFromCands::TrackTable EleTkIsolFromCands::preselectTracksFromCands(pat::PackedCandidateCollection const& cands,
0075                                                                             TrkCuts const& cuts,
0076                                                                             PIDVeto pidVeto) {
0077   std::vector<float> pt;
0078   std::vector<float> eta;
0079   std::vector<float> phi;
0080   std::vector<float> vz;
0081   pt.reserve(cands.size());
0082   eta.reserve(cands.size());
0083   phi.reserve(cands.size());
0084   vz.reserve(cands.size());
0085
0086   for (auto const& cand : cands) {
0087     if (cand.hasTrackDetails() && cand.charge() != 0 && passPIDVeto(cand.pdgId(), pidVeto)) {
0088       const reco::Track& trk = cand.pseudoTrack();
0089       float trkPt = trk.pt();
0090       if (passTrackPreselection(trk, trkPt, cuts)) {
0091         pt.emplace_back(trkPt);
0092         eta.emplace_back(trk.eta());
0093         phi.emplace_back(trk.phi());
0094         vz.emplace_back(trk.vz());
0095       }
0096     }
0097   }
0098
0099   return {pt, std::move(eta), std::move(phi), std::move(vz)};
0100 }
0101
0102 EleTkIsolFromCands::Output EleTkIsolFromCands::operator()(const reco::TrackBase& eleTrk) {
0103   using namespace edm::soa::col;
0104
0105   float ptSum = 0.;
0106   int nrTrks = 0;
0107
0108   const float eleEta = eleTrk.eta();
0109   const float elePhi = eleTrk.phi();
0110   const float eleVz = eleTrk.vz();
0111
0112   const bool isBarrelElectron = std::abs(eleEta) < 1.5;
0113
0114   auto const& preselectedTracks = getPreselectedTracks(isBarrelElectron);
0115   auto const& cuts = isBarrelElectron ? cfg_.barrelCuts : cfg_.endcapCuts;
0116
0117   for (auto const& trk : preselectedTracks) {
0118     const float dR2 = reco::deltaR2(eleEta, elePhi, trk.get<Eta>(), trk.get<Phi>());
0119     const float dEta = trk.get<Eta>() - eleEta;
0120     const float dZ = eleVz - trk.get<Vz>();
0121
0122     if (dR2 >= cuts.minDR2 && dR2 <= cuts.maxDR2 && std::abs(dEta) >= cuts.minDEta && std::abs(dZ) < cuts.maxDZ) {
0123       ptSum += trk.get<Pt>();
0124       nrTrks++;
0125     }
0126   }
0127   return {nrTrks, ptSum};
0128 }
0129
0130 bool EleTkIsolFromCands::passPIDVeto(const int pdgId, const EleTkIsolFromCands::PIDVeto veto) {
0131   int pidAbs = std::abs(pdgId);
0132   switch (veto) {
0133     case PIDVeto::NONE:
0134       return true;
0135     case PIDVeto::ELES:
0136       if (pidAbs == 11)
0137         return false;
0138       else
0139         return true;
0140     case PIDVeto::NONELES:
0141       if (pidAbs == 11)
0142         return true;
0143       else
0144         return false;
0145   }
0146   throw cms::Exception("CodeError") << "invalid PIDVeto " << static_cast<int>(veto) << ", "
0147                                     << "this is likely due to some static casting of invalid ints somewhere";
0148 }
0149
0150 EleTkIsolFromCands::PIDVeto EleTkIsolFromCands::pidVetoFromStr(const std::string& vetoStr) {
0151   if (vetoStr == "NONE")
0152     return PIDVeto::NONE;
0153   else if (vetoStr == "ELES")
0154     return PIDVeto::ELES;
0155   else if (vetoStr == "NONELES")
0156     return PIDVeto::NONELES;
0157   else {
0158     throw cms::Exception("CodeError") << "unrecognised string " << vetoStr
0159                                       << ", either a typo or this function needs to be updated";
0160   }
0161 }
0162
0163 bool EleTkIsolFromCands::passTrackPreselection(const reco::TrackBase& trk, float trackPt, const TrkCuts& cuts) {
0164   return trackPt > cuts.minPt && trk.numberOfValidHits() >= cuts.minHits &&
0165          trk.hitPattern().numberOfValidPixelHits() >= cuts.minPixelHits && passQual(trk, cuts.allowedQualities) &&
0166          passAlgo(trk, cuts.algosToReject) && (cuts.maxDPtPt < 0 || trk.ptError() / trackPt < cuts.maxDPtPt);
0167 }
0168
0169 bool EleTkIsolFromCands::passQual(const reco::TrackBase& trk, const std::vector<reco::TrackBase::TrackQuality>& quals) {
0170   if (quals.empty())
0171     return true;
0172
0173   for (auto qual : quals) {
0174     if (trk.quality(qual))
0175       return true;
0176   }
0177
0178   return false;
0179 }
0180
0181 bool EleTkIsolFromCands::passAlgo(const reco::TrackBase& trk,
0182                                   const std::vector<reco::TrackBase::TrackAlgorithm>& algosToRej) {
0183   return algosToRej.empty() ||
0184          //check also the originalAlgo in case the track is reconstructed by more than one
0185          //reject only if both algo and originalAlgo are not acceptable
0186          !(std::binary_search(algosToRej.begin(), algosToRej.end(), trk.algo()) &&
0187            std::binary_search(algosToRej.begin(), algosToRej.end(), trk.originalAlgo()));
0188 }
0189
0190 EleTkIsolFromCands::TrackTable const& EleTkIsolFromCands::getPreselectedTracks(bool isBarrel) {
0191   auto const& cuts = isBarrel ? cfg_.barrelCuts : cfg_.endcapCuts;
0192   auto& preselectedTracks = isBarrel ? preselectedTracksWithBarrelCuts_ : preselectedTracksWithEndcapCuts_;
0193   bool& tracksCached = isBarrel ? tracksCachedForBarrelCuts_ : tracksCachedForEndcapCuts_;
0194
0195   if (!tracksCached) {
0196     preselectedTracks = tracks_ ? preselectTracks(*tracks_, cuts) : preselectTracksFromCands(*cands_, cuts, pidVeto_);
0197     tracksCached = true;
0198   }
0199
0200   return preselectedTracks;
0201 }