MuonIsolation/src/IsolatorByNominalEfficiency.cc

0001 #include "RecoMuon/MuonIsolation/interface/IsolatorByNominalEfficiency.h"
0002 #include "RecoMuon/MuonIsolation/src/NominalEfficiencyThresholds.h"
0003 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0004 #include "FWCore/ParameterSet/interface/FileInPath.h"
0005
0006 using namespace muonisolation;
0007 using namespace std;
0008
0009 double IsolatorByNominalEfficiency::ConeSizes::size(int i) const {
0010   if (i >= 0 && i < DIM)
0011     return cone_dr[i];
0012   else
0013     return 0.;
0014 }
0015
0016 int IsolatorByNominalEfficiency::ConeSizes::index(float dr) const {
0017   for (int i = DIM; i >= 1; i--)
0018     if (cone_dr[i - 1] < dr)
0019       return i;
0020   return 0;
0021 }
0022
0023 const float IsolatorByNominalEfficiency::ConeSizes::cone_dr[] = {
0024     0.001, 0.02, 0.045, 0.09, 0.13, 0.17, 0.20, 0.24, 0.28, 0.32, 0.38, 0.45, 0.5, 0.6, 0.7};
0025
0026 IsolatorByNominalEfficiency::IsolatorByNominalEfficiency(const string& thrFile,
0027                                                          const vector<string>& ceff,
0028                                                          const vector<double>& weights)
0029     : theWeights(weights) {
0030   thresholds = new NominalEfficiencyThresholds(findPath(thrFile));
0031   coneForEfficiency = cones(ceff);
0032   theDepThresholds = std::vector<double>(weights.size(), -1e12);
0033 }
0034
0035 IsolatorByNominalEfficiency::IsolatorByNominalEfficiency(const string& thrFile,
0036                                                          const vector<string>& ceff,
0037                                                          const vector<double>& weights,
0038                                                          const vector<double>& thresh)
0039     : theWeights(weights), theDepThresholds(thresh) {
0040   thresholds = new NominalEfficiencyThresholds(findPath(thrFile));
0041   coneForEfficiency = cones(ceff);
0042 }
0043
0044 IsolatorByNominalEfficiency::~IsolatorByNominalEfficiency() { delete thresholds; }
0045
0046 IsolatorByNominalEfficiency::mapNomEff_Cone IsolatorByNominalEfficiency::cones(const vector<string>& usrVec) {
0047   mapNomEff_Cone result;
0048   for (vector<string>::const_iterator is = usrVec.begin(); is != usrVec.end(); is++) {
0049     char* evp = nullptr;
0050     int cone = strtol((*is).c_str(), &evp, 10);
0051     float effic = strtod(evp + 1, &evp);
0052     result.insert(make_pair(effic, cone));
0053   }
0054   return result;
0055 }
0056
0057 string IsolatorByNominalEfficiency::findPath(const string& fileName) {
0058   // FIXME
0059   edm::FileInPath f(fileName);
0060   return f.fullPath();
0061 }
0062
0063 MuIsoBaseIsolator::Result IsolatorByNominalEfficiency::result(const DepositContainer& deposits,
0064                                                               const edm::Event*) const {
0065   if (deposits.empty()) {
0066     cout << "IsolatorByNominalEfficiency: no deposit" << endl;
0067     return Result(resultType());  //FIXME
0068   }
0069
0070   // To determine the threshold, the direction of the cone of the first
0071   // set of deposits is used.
0072   // For algorithms where different cone axis definitions are used
0073   // for different types deposits (eg. HCAL and ECAL deposits for
0074   // calorimeter isolation), the first one is used to determine the threshold
0075   // value!
0076   float theEta = deposits.back().dep->eta();
0077
0078   // Try descending efficiency values to find the point where the candidate
0079   // becomes non isolated
0080
0081   float nominalEfficiency = 1.;
0082   const float deltaeff = 0.005;
0083   const float mineff = deltaeff;
0084   for (float eff = .995; eff > mineff; eff -= deltaeff) {
0085     int cone = bestConeForEfficiencyIndex(eff);
0086     float coneSize = theConesInfo.size(cone);
0087     NominalEfficiencyThresholds::ThresholdLocation location = {theEta, cone};
0088     float thres = thresholds->thresholdValueForEfficiency(location, eff);
0089     float sumDep = weightedSum(deposits, coneSize);
0090     //      cout << "   Eff=" << eff
0091     //         << " eta=" << theEta
0092     //         << " cone=" << cone
0093     //         << " dR=" << coneSize
0094     //         << " thres=" << thres
0095     //         << " deposit=" << sumDep
0096     //         << " isolated=" << (sumDep < thres)
0097     //         << endl;
0098     if (sumDep > thres)
0099       break;
0100     nominalEfficiency = eff;
0101   }
0102   Result res(resultType());
0103   res.valFloat = nominalEfficiency;
0104   return res;
0105 }
0106
0107 int IsolatorByNominalEfficiency::bestConeForEfficiencyIndex(float eff_thr) const {
0108   //FIXME use upper_bound
0109   int best_cone;
0110   if (!coneForEfficiency.empty()) {
0111     best_cone = (--(coneForEfficiency.end()))->second;
0112   } else
0113     return 0;
0114
0115   mapNomEff_Cone::const_reverse_iterator it;
0116   for (it = coneForEfficiency.rbegin(); it != coneForEfficiency.rend(); it++) {
0117     if (eff_thr <= (*it).first)
0118       best_cone = (*it).second;
0119   }
0120   return best_cone;
0121 }
0122
0123 double IsolatorByNominalEfficiency::weightedSum(const DepositContainer& deposits, float dRcone) const {
0124   double sumDep = 0;
0125
0126   assert(deposits.size() == theWeights.size());
0127
0128   vector<double>::const_iterator w = theWeights.begin();
0129   vector<double>::const_iterator dThresh = theDepThresholds.begin();
0130   for (DepositContainer::const_iterator dep = deposits.begin(); dep != deposits.end(); dep++) {
0131     if (dep->vetos != nullptr) {
0132       sumDep += dep->dep->depositAndCountWithin(dRcone, *dep->vetos, *dThresh).first * (*w);
0133     } else {
0134       sumDep += dep->dep->depositAndCountWithin(dRcone, Vetos(), *dThresh).first * (*w);
0135     }
0136     if (sumDep < 0.)
0137       sumDep = 0.;
0138     w++;
0139     dThresh++;
0140   }
0141   return sumDep;
0142 }
0143
0144 Cuts IsolatorByNominalEfficiency::cuts(float nominalEfficiency) const {
0145   vector<double> etaBounds = thresholds->bins();
0146   vector<double> coneSizes;
0147   vector<double> cutvalues;
0148   for (vector<double>::const_iterator it = etaBounds.begin(), itEnd = etaBounds.end(); it < itEnd; ++it) {
0149     float eta = (*it);
0150     int icone = bestConeForEfficiencyIndex(nominalEfficiency);
0151     coneSizes.push_back(theConesInfo.size(icone));
0152     NominalEfficiencyThresholds::ThresholdLocation location = {eta - 1.e-3f, icone};
0153     cutvalues.push_back(thresholds->thresholdValueForEfficiency(location, nominalEfficiency));
0154   }
0155   return Cuts(etaBounds, coneSizes, cutvalues);
0156 }