RecoTau/interface/TauWPThreshold.h

0001 #ifndef RecoTauTag_RecoTau_TauWPThreshold_h
0002 #define RecoTauTag_RecoTau_TauWPThreshold_h
0003
0004 #include "DataFormats/TauReco/interface/BaseTau.h"
0005 #include "DataFormats/PatCandidates/interface/Tau.h"
0006
0007 #include <TF1.h>
0008
0009 #include <cmath>
0010 #include <cstdlib>
0011
0012 namespace tau {
0013   class TauWPThreshold {
0014   public:
0015     explicit TauWPThreshold(const std::string& cut_str) {
0016       bool simple_value = false;
0017       const char* cut_cstr = cut_str.c_str();
0018       char* end_cstr{};
0019       value_ = std::strtod(cut_cstr, &end_cstr);
0020       // simple number if end_cstr is empty and not equal to initial string (cut_cstr), and returned value is finite
0021       simple_value = !*end_cstr && cut_cstr != end_cstr && std::isfinite(value_);
0022       if (!simple_value) {
0023         static const std::string prefix =
0024             "[&](double *x, double *p) { const int decayMode = p[0];"
0025             "const double pt = p[1]; const double eta = p[2];";
0026         static const int n_params = 3;
0027         static const auto handler = [](int, Bool_t, const char*, const char*) -> void {};
0028
0029         std::string fn_str = prefix;
0030         if (cut_str.find("return") == std::string::npos)
0031           fn_str += " return " + cut_str + ";}";
0032         else
0033           fn_str += cut_str + "}";
0034         auto old_handler = SetErrorHandler(handler);
0035         fn_ = std::make_unique<TF1>("fn_", fn_str.c_str(), 0, 1, n_params);
0036         SetErrorHandler(old_handler);
0037         if (!fn_->IsValid())
0038           throw cms::Exception("TauWPThreshold: invalid formula") << "Invalid WP cut formula = '" << cut_str << "'.";
0039       }
0040     }
0041
0042     double operator()(int dm, double pt, double eta) const {
0043       if (!fn_)
0044         return value_;
0045
0046       fn_->SetParameter(0, dm);
0047       fn_->SetParameter(1, pt);
0048       fn_->SetParameter(2, eta);
0049       return fn_->Eval(0);
0050     }
0051
0052     double operator()(const reco::BaseTau& tau, bool isPFTau) const {
0053       const int dm =
0054           isPFTau ? dynamic_cast<const reco::PFTau&>(tau).decayMode() : dynamic_cast<const pat::Tau&>(tau).decayMode();
0055       return (*this)(dm, tau.pt(), tau.eta());
0056     }
0057
0058     double operator()(const reco::Candidate& tau) const { return (*this)(-1, tau.pt(), tau.eta()); }
0059
0060   private:
0061     std::unique_ptr<TF1> fn_;
0062     double value_;
0063   };
0064 }  // namespace tau
0065
0066 #endif