EgammaIsolationAlgos/src/EgammaL1TkIsolation.cc

0001 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaL1TkIsolation.h"
0002 #include "DataFormats/TrackReco/interface/Track.h"
0003 #include "DataFormats/Math/interface/deltaR.h"
0004
0005 #include <algorithm>
0006
0007 EgammaL1TkIsolation::EgammaL1TkIsolation(const edm::ParameterSet& para)
0008     : useAbsEta_(para.getParameter<bool>("useAbsEta")),
0009       etaBoundaries_(para.getParameter<std::vector<double>>("etaBoundaries")) {
0010   const auto& trkCutParams = para.getParameter<std::vector<edm::ParameterSet>>("trkCuts");
0011   for (const auto& params : trkCutParams) {
0012     trkCuts_.emplace_back(TrkCuts(params));
0013   }
0014   if (etaBoundaries_.size() + 1 != trkCuts_.size()) {
0015     throw cms::Exception("ConfigError") << "EgammaL1TkIsolation: etaBoundaries parameters size ("
0016                                         << etaBoundaries_.size()
0017                                         << ") should be one less than the size of trkCuts VPSet (" << trkCuts_.size()
0018                                         << ")";
0019   }
0020   if (!std::is_sorted(etaBoundaries_.begin(), etaBoundaries_.end())) {
0021     throw cms::Exception("ConfigError")
0022         << "EgammaL1TkIsolation: etaBoundaries parameter's entries should be in increasing value";
0023   }
0024 }
0025
0026 void EgammaL1TkIsolation::fillPSetDescription(edm::ParameterSetDescription& desc) {
0027   desc.add("useAbsEta", true);
0028   desc.add("etaBoundaries", std::vector{1.5});
0029   desc.addVPSet("trkCuts", TrkCuts::makePSetDescription(), {edm::ParameterSet(), edm::ParameterSet()});
0030 }
0031
0032 std::pair<int, double> EgammaL1TkIsolation::calIsol(const reco::TrackBase& trk, const L1TrackCollection& tracks) const {
0033   return calIsol(trk.eta(), trk.phi(), trk.vz(), tracks);
0034 }
0035
0036 std::pair<int, double> EgammaL1TkIsolation::calIsol(const double objEta,
0037                                                     const double objPhi,
0038                                                     const double objZ,
0039                                                     const L1TrackCollection& tracks) const {
0040   double ptSum = 0.;
0041   int nrTrks = 0;
0042
0043   const TrkCuts& cuts = trkCuts_[etaBinNr(objEta)];
0044
0045   for (const auto& trk : tracks) {
0046     const float trkPt = trk.momentum().perp();
0047     if (passTrkSel(trk, trkPt, cuts, objEta, objPhi, objZ)) {
0048       ptSum += trkPt;
0049       nrTrks++;
0050     }
0051   }
0052   return {nrTrks, ptSum};
0053 }
0054
0055 EgammaL1TkIsolation::TrkCuts::TrkCuts(const edm::ParameterSet& para) {
0056   minPt = para.getParameter<double>("minPt");
0057   auto sq = [](double val) { return val * val; };
0058   minDR2 = sq(para.getParameter<double>("minDR"));
0059   maxDR2 = sq(para.getParameter<double>("maxDR"));
0060   minDEta = para.getParameter<double>("minDEta");
0061   maxDZ = para.getParameter<double>("maxDZ");
0062 }
0063
0064 edm::ParameterSetDescription EgammaL1TkIsolation::TrkCuts::makePSetDescription() {
0065   edm::ParameterSetDescription desc;
0066   desc.add<double>("minPt", 2.0);
0067   desc.add<double>("maxDR", 0.3);
0068   desc.add<double>("minDR", 0.01);
0069   desc.add<double>("minDEta", 0.003);
0070   desc.add<double>("maxDZ", 0.7);
0071   return desc;
0072 }
0073
0074 //as we have verfied that trkCuts_ size is etaBoundaries_ size +1
0075 //then this is always a valid binnr for trkCuts_
0076 size_t EgammaL1TkIsolation::etaBinNr(double eta) const {
0077   if (useAbsEta_) {
0078     eta = std::abs(eta);
0079   }
0080   auto res = std::upper_bound(etaBoundaries_.begin(), etaBoundaries_.end(), eta);
0081   size_t binNr = std::distance(etaBoundaries_.begin(), res);
0082   return binNr;
0083 }
0084
0085 bool EgammaL1TkIsolation::passTrkSel(const L1Track& trk,
0086                                      const double trkPt,
0087                                      const TrkCuts& cuts,
0088                                      const double objEta,
0089                                      const double objPhi,
0090                                      const double objZ) {
0091   if (trkPt > cuts.minPt && std::abs(objZ - trk.z0()) < cuts.maxDZ) {
0092     const float trkEta = trk.eta();
0093     const float dEta = trkEta - objEta;
0094     const float dR2 = reco::deltaR2(objEta, objPhi, trkEta, trk.phi());
0095     return dR2 >= cuts.minDR2 && dR2 <= cuts.maxDR2 && std::abs(dEta) >= cuts.minDEta;
0096   }
0097
0098   return false;
0099 }