Phase2L1ParticleFlow/src/ParametricResolution.cc

0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/ParametricResolution.h"
0002
0003 #ifdef CMSSW_GIT_HASH
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/Utilities/interface/Exception.h"
0006
0007 std::vector<float> l1tpf::ParametricResolution::getVFloat(const edm::ParameterSet &cpset, const std::string &name) {
0008   std::vector<double> vd = cpset.getParameter<std::vector<double>>(name);
0009   return std::vector<float>(vd.begin(), vd.end());
0010 }
0011
0012 l1tpf::ParametricResolution::ParametricResolution(const edm::ParameterSet &cpset)
0013     : etas_(getVFloat(cpset, "etaBins")), offsets_(getVFloat(cpset, "offset")), scales_(getVFloat(cpset, "scale")) {
0014   if (cpset.existsAs<std::vector<double>>("ptMin")) {
0015     ptMins_ = getVFloat(cpset, "ptMin");
0016   } else {
0017     float ptMin = cpset.existsAs<double>("ptMin") ? cpset.getParameter<double>("ptMin") : 0;
0018     ptMins_ = std::vector<float>(etas_.size(), ptMin);
0019   }
0020   if (cpset.existsAs<std::vector<double>>("ptMax")) {
0021     ptMaxs_ = getVFloat(cpset, "ptMax");
0022   } else {
0023     ptMaxs_ = std::vector<float>(etas_.size(), 1e6);
0024   }
0025
0026   std::string skind = cpset.getParameter<std::string>("kind");
0027   if (skind == "track")
0028     kind_ = Kind::Track;
0029   else if (skind == "calo")
0030     kind_ = Kind::Calo;
0031   else
0032     throw cms::Exception("Configuration", "Bad kind of resolution: " + skind);
0033 }
0034 #endif
0035
0036 float l1tpf::ParametricResolution::operator()(const float pt, const float abseta) const {
0037   for (unsigned int i = 0, n = etas_.size(); i < n; ++i) {
0038     if (pt > ptMaxs_[i])
0039       continue;
0040     if (abseta < etas_[i]) {
0041       switch (kind_) {
0042         case Kind::Track:
0043           return pt * std::min<float>(1.f, std::hypot(pt * scales_[i] * 0.001, offsets_[i]));
0044         case Kind::Calo:
0045           return std::min<float>(pt, pt * scales_[i] + offsets_[i]);
0046           if (pt < ptMins_[i])
0047             return pt * std::min<float>(1, scales_[i] + offsets_[i] / ptMins_[i]);
0048           return std::min<float>(pt, pt * scales_[i] + offsets_[i]);
0049       }
0050     }
0051   }
0052   return std::min<float>(pt, 0.3 * pt + 7);  // saturate to 100% at 10 GeV, and to 30% at high pt
0053 }