IsolatorByNominalEfficiency.cc

Go to the documentation of this file.

#include "RecoMuon/MuonIsolation/interface/IsolatorByNominalEfficiency.h"
#include "RecoMuon/MuonIsolation/src/NominalEfficiencyThresholds.h"
#include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"

using namespace muonisolation;
using namespace std;

double IsolatorByNominalEfficiency::ConeSizes::size(int i) const {
  if (i >= 0 && i < DIM)
    return cone_dr[i];
  else
    return 0.;
}

int IsolatorByNominalEfficiency::ConeSizes::index(float dr) const {
  for (int i = DIM; i >= 1; i--)
    if (cone_dr[i - 1] < dr)
      return i;
  return 0;
}

const float IsolatorByNominalEfficiency::ConeSizes::cone_dr[] = {
    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};

IsolatorByNominalEfficiency::IsolatorByNominalEfficiency(const string& thrFile,
                                                         const vector<string>& ceff,
                                                         const vector<double>& weights)
    : theWeights(weights) {
  thresholds = new NominalEfficiencyThresholds(findPath(thrFile));
  coneForEfficiency = cones(ceff);
  theDepThresholds = std::vector<double>(weights.size(), -1e12);
}

IsolatorByNominalEfficiency::IsolatorByNominalEfficiency(const string& thrFile,
                                                         const vector<string>& ceff,
                                                         const vector<double>& weights,
                                                         const vector<double>& thresh)
    : theWeights(weights), theDepThresholds(thresh) {
  thresholds = new NominalEfficiencyThresholds(findPath(thrFile));
  coneForEfficiency = cones(ceff);
}

IsolatorByNominalEfficiency::~IsolatorByNominalEfficiency() { delete thresholds; }

IsolatorByNominalEfficiency::mapNomEff_Cone IsolatorByNominalEfficiency::cones(const vector<string>& usrVec) {
  mapNomEff_Cone result;
  for (vector<string>::const_iterator is = usrVec.begin(); is != usrVec.end(); is++) {
    char* evp = nullptr;
    int cone = strtol((*is).c_str(), &evp, 10);
    float effic = strtod(evp + 1, &evp);
    result.insert(make_pair(effic, cone));
  }
  return result;
}

string IsolatorByNominalEfficiency::findPath(const string& fileName) {
  // FIXME
  edm::FileInPath f(fileName);
  return f.fullPath();
}

MuIsoBaseIsolator::Result IsolatorByNominalEfficiency::result(const DepositContainer& deposits,
                                                              const edm::Event*) const {
  if (deposits.empty()) {
    cout << "IsolatorByNominalEfficiency: no deposit" << endl;
    return Result(resultType());  //FIXME
  }

  // To determine the threshold, the direction of the cone of the first
  // set of deposits is used.
  // For algorithms where different cone axis definitions are used
  // for different types deposits (eg. HCAL and ECAL deposits for
  // calorimeter isolation), the first one is used to determine the threshold
  // value!
  float theEta = deposits.back().dep->eta();

  // Try descending efficiency values to find the point where the candidate
  // becomes non isolated

  float nominalEfficiency = 1.;
  const float deltaeff = 0.005;
  const float mineff = deltaeff;
  for (float eff = .995; eff > mineff; eff -= deltaeff) {
    int cone = bestConeForEfficiencyIndex(eff);
    float coneSize = theConesInfo.size(cone);
    NominalEfficiencyThresholds::ThresholdLocation location = {theEta, cone};
    float thres = thresholds->thresholdValueForEfficiency(location, eff);
    float sumDep = weightedSum(deposits, coneSize);
    //      cout << "   Eff=" << eff
    //         << " eta=" << theEta
    //         << " cone=" << cone
    //         << " dR=" << coneSize
    //         << " thres=" << thres
    //         << " deposit=" << sumDep
    //         << " isolated=" << (sumDep < thres)
    //         << endl;
    if (sumDep > thres)
      break;
    nominalEfficiency = eff;
  }
  Result res(resultType());
  res.valFloat = nominalEfficiency;
  return res;
}

int IsolatorByNominalEfficiency::bestConeForEfficiencyIndex(float eff_thr) const {
  //FIXME use upper_bound
  int best_cone;
  if (!coneForEfficiency.empty()) {
    best_cone = (--(coneForEfficiency.end()))->second;
  } else
    return 0;

  mapNomEff_Cone::const_reverse_iterator it;
  for (it = coneForEfficiency.rbegin(); it != coneForEfficiency.rend(); it++) {
    if (eff_thr <= (*it).first)
      best_cone = (*it).second;
  }
  return best_cone;
}

double IsolatorByNominalEfficiency::weightedSum(const DepositContainer& deposits, float dRcone) const {
  double sumDep = 0;

  assert(deposits.size() == theWeights.size());

  vector<double>::const_iterator w = theWeights.begin();
  vector<double>::const_iterator dThresh = theDepThresholds.begin();
  for (DepositContainer::const_iterator dep = deposits.begin(); dep != deposits.end(); dep++) {
    if (dep->vetos != nullptr) {
      sumDep += dep->dep->depositAndCountWithin(dRcone, *dep->vetos, *dThresh).first * (*w);
    } else {
      sumDep += dep->dep->depositAndCountWithin(dRcone, Vetos(), *dThresh).first * (*w);
    }
    if (sumDep < 0.)
      sumDep = 0.;
    w++;
    dThresh++;
  }
  return sumDep;
}

Cuts IsolatorByNominalEfficiency::cuts(float nominalEfficiency) const {
  vector<double> etaBounds = thresholds->bins();
  vector<double> coneSizes;
  vector<double> cutvalues;
  for (vector<double>::const_iterator it = etaBounds.begin(), itEnd = etaBounds.end(); it < itEnd; ++it) {
    float eta = (*it);
    int icone = bestConeForEfficiencyIndex(nominalEfficiency);
    coneSizes.push_back(theConesInfo.size(icone));
    NominalEfficiencyThresholds::ThresholdLocation location = {eta - 1.e-3f, icone};
    cutvalues.push_back(thresholds->thresholdValueForEfficiency(location, nominalEfficiency));
  }
  return Cuts(etaBounds, coneSizes, cutvalues);
}