File indexing completed on 2023-03-17 11:20:39
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
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());
0068 }
0069
0070
0071
0072
0073
0074
0075
0076 float theEta = deposits.back().dep->eta();
0077
0078
0079
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
0091
0092
0093
0094
0095
0096
0097
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
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 }