Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-10 22:24:46

0001 #include <limits>
0002 #include "FWCore/ParameterSet/interface/FileInPath.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "L1Trigger/L1THGCal/interface/backend/HGCalTriggerClusterIdentificationBase.h"
0005 #include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h"
0006 #include "CommonTools/MVAUtils/interface/TMVAEvaluator.h"
0007 
0008 class HGCalTriggerClusterIdentificationBDT : public HGCalTriggerClusterIdentificationBase {
0009 public:
0010   class Category {
0011   public:
0012     Category(float pt_min, float pt_max, float eta_min, float eta_max)
0013         : pt_min_(pt_min), pt_max_(pt_max), eta_min_(eta_min), eta_max_(eta_max) {}
0014     ~Category() {}
0015     bool contains(float pt, float eta) const {
0016       bool output = true;
0017       if (pt < pt_min_ || pt >= pt_max_)
0018         output = false;
0019       if (std::abs(eta) < eta_min_ || std::abs(eta) >= eta_max_)
0020         output = false;
0021       return output;
0022     }
0023 
0024     float pt_min() const { return pt_min_; }
0025     float pt_max() const { return pt_max_; }
0026     float eta_min() const { return eta_min_; }
0027     float eta_max() const { return eta_max_; }
0028 
0029   private:
0030     float pt_min_ = 0.;
0031     float pt_max_ = std::numeric_limits<float>::max();
0032     float eta_min_ = 1.5;
0033     float eta_max_ = 3.;
0034   };
0035 
0036   class WorkingPoint {
0037   public:
0038     WorkingPoint(const std::string& name, double wp) : name_(name), wp_(wp) {}
0039     ~WorkingPoint() = default;
0040 
0041     const std::string& name() const { return name_; }
0042     double working_point() const { return wp_; }
0043 
0044   private:
0045     std::string name_;
0046     double wp_;
0047   };
0048 
0049 public:
0050   HGCalTriggerClusterIdentificationBDT();
0051   ~HGCalTriggerClusterIdentificationBDT() override{};
0052   void initialize(const edm::ParameterSet& conf) final;
0053   float value(const l1t::HGCalMulticluster& cluster) const final;
0054   bool decision(const l1t::HGCalMulticluster& cluster, unsigned wp = 0) const final;
0055   const std::vector<std::string>& working_points() const final { return working_points_names_; }
0056 
0057 private:
0058   enum class ClusterVariable {
0059     cl3d_showerlength,
0060     cl3d_coreshowerlength,
0061     cl3d_firstlayer,
0062     cl3d_maxlayer,
0063     cl3d_seetot,
0064     cl3d_seemax,
0065     cl3d_spptot,
0066     cl3d_sppmax,
0067     cl3d_szz,
0068     cl3d_srrtot,
0069     cl3d_srrmax,
0070     cl3d_srrmean
0071   };
0072   std::vector<Category> categories_;
0073   std::vector<std::unique_ptr<TMVAEvaluator>> bdts_;
0074   std::vector<std::vector<WorkingPoint>> working_points_;
0075   std::vector<std::string> input_variables_;
0076   std::vector<ClusterVariable> input_variables_id_;
0077   std::vector<std::string> working_points_names_;
0078 
0079   float clusterVariable(ClusterVariable, const l1t::HGCalMulticluster&) const;
0080   int category(float pt, float eta) const;
0081 };
0082 
0083 DEFINE_HGC_TPG_CLUSTER_ID(HGCalTriggerClusterIdentificationBDT, "HGCalTriggerClusterIdentificationBDT");
0084 
0085 HGCalTriggerClusterIdentificationBDT::HGCalTriggerClusterIdentificationBDT()
0086     : HGCalTriggerClusterIdentificationBase() {}
0087 
0088 void HGCalTriggerClusterIdentificationBDT::initialize(const edm::ParameterSet& conf) {
0089   if (!bdts_.empty()) {
0090     edm::LogWarning("HGCalTriggerClusterIdentificationBDT|Initialization") << "BDTs already initialized.";
0091     return;
0092   }
0093   input_variables_ = conf.getParameter<std::vector<std::string>>("Inputs");
0094   std::vector<std::string> bdt_files = conf.getParameter<std::vector<std::string>>("Weights");
0095   std::vector<double> categories_etamin = conf.getParameter<std::vector<double>>("CategoriesEtaMin");
0096   std::vector<double> categories_etamax = conf.getParameter<std::vector<double>>("CategoriesEtaMax");
0097   std::vector<double> categories_ptmin = conf.getParameter<std::vector<double>>("CategoriesPtMin");
0098   std::vector<double> categories_ptmax = conf.getParameter<std::vector<double>>("CategoriesPtMax");
0099 
0100   if (bdt_files.size() != categories_etamin.size() || categories_etamin.size() != categories_etamax.size() ||
0101       categories_etamax.size() != categories_ptmin.size() || categories_ptmin.size() != categories_ptmax.size()) {
0102     throw cms::Exception("HGCalTriggerClusterIdentificationBDT|BadInitialization")
0103         << "Inconsistent numbers of categories, BDT weight files and working points";
0104   }
0105   size_t categories_size = categories_etamin.size();
0106 
0107   const auto wps_conf = conf.getParameter<std::vector<edm::ParameterSet>>("WorkingPoints");
0108   working_points_.resize(categories_size);
0109   for (const auto& wp_conf : wps_conf) {
0110     std::string wp_name = wp_conf.getParameter<std::string>("Name");
0111     std::vector<double> wps = wp_conf.getParameter<std::vector<double>>("WorkingPoint");
0112     working_points_names_.emplace_back(wp_name);
0113     if (wps.size() != categories_size) {
0114       throw cms::Exception("HGCalTriggerClusterIdentificationBDT|BadInitialization")
0115           << "Inconsistent number of categories in working point '" << wp_name << "'";
0116     }
0117     for (size_t cat = 0; cat < categories_size; cat++) {
0118       working_points_[cat].emplace_back(wp_name, wps[cat]);
0119     }
0120   }
0121 
0122   categories_.reserve(categories_size);
0123   bdts_.reserve(categories_size);
0124   for (size_t cat = 0; cat < categories_size; cat++) {
0125     categories_.emplace_back(
0126         categories_ptmin[cat], categories_ptmax[cat], categories_etamin[cat], categories_etamax[cat]);
0127   }
0128   std::vector<std::string> spectators = {};
0129   for (const auto& file : bdt_files) {
0130     bdts_.emplace_back(new TMVAEvaluator());
0131     bdts_.back()->initialize("!Color:Silent:!Error",
0132                              "BDT::bdt",
0133                              edm::FileInPath(file).fullPath(),
0134                              input_variables_,
0135                              spectators,
0136                              false,
0137                              false);
0138   }
0139 
0140   // Transform input variable strings to enum values for later comparisons
0141   input_variables_id_.reserve(input_variables_.size());
0142   for (const auto& variable : input_variables_) {
0143     if (variable == "cl3d_showerlength")
0144       input_variables_id_.push_back(ClusterVariable::cl3d_showerlength);
0145     else if (variable == "cl3d_coreshowerlength")
0146       input_variables_id_.push_back(ClusterVariable::cl3d_coreshowerlength);
0147     else if (variable == "cl3d_firstlayer")
0148       input_variables_id_.push_back(ClusterVariable::cl3d_firstlayer);
0149     else if (variable == "cl3d_maxlayer")
0150       input_variables_id_.push_back(ClusterVariable::cl3d_maxlayer);
0151     else if (variable == "cl3d_seetot")
0152       input_variables_id_.push_back(ClusterVariable::cl3d_seetot);
0153     else if (variable == "cl3d_seemax")
0154       input_variables_id_.push_back(ClusterVariable::cl3d_seemax);
0155     else if (variable == "cl3d_spptot")
0156       input_variables_id_.push_back(ClusterVariable::cl3d_spptot);
0157     else if (variable == "cl3d_sppmax")
0158       input_variables_id_.push_back(ClusterVariable::cl3d_sppmax);
0159     else if (variable == "cl3d_szz")
0160       input_variables_id_.push_back(ClusterVariable::cl3d_szz);
0161     else if (variable == "cl3d_srrtot")
0162       input_variables_id_.push_back(ClusterVariable::cl3d_srrtot);
0163     else if (variable == "cl3d_srrmax")
0164       input_variables_id_.push_back(ClusterVariable::cl3d_srrmax);
0165     else if (variable == "cl3d_srrmean")
0166       input_variables_id_.push_back(ClusterVariable::cl3d_srrmean);
0167   }
0168 }
0169 
0170 float HGCalTriggerClusterIdentificationBDT::value(const l1t::HGCalMulticluster& cluster) const {
0171   std::map<std::string, float> inputs;
0172   for (unsigned i = 0; i < input_variables_.size(); i++) {
0173     inputs[input_variables_[i]] = clusterVariable(input_variables_id_[i], cluster);
0174   }
0175   float pt = cluster.pt();
0176   float eta = cluster.eta();
0177   int cat = category(pt, eta);
0178   return (cat != -1 ? bdts_.at(cat)->evaluate(inputs) : -999.);
0179 }
0180 
0181 bool HGCalTriggerClusterIdentificationBDT::decision(const l1t::HGCalMulticluster& cluster, unsigned wp) const {
0182   float bdt_output = value(cluster);
0183   float pt = cluster.pt();
0184   float eta = cluster.eta();
0185   int cat = category(pt, eta);
0186   return (cat != -1 ? bdt_output > working_points_.at(cat).at(wp).working_point() : true);
0187 }
0188 
0189 int HGCalTriggerClusterIdentificationBDT::category(float pt, float eta) const {
0190   for (unsigned cat = 0; cat < categories_.size(); cat++) {
0191     if (categories_[cat].contains(pt, eta))
0192       return static_cast<int>(cat);
0193   }
0194   return -1;
0195 }
0196 
0197 float HGCalTriggerClusterIdentificationBDT::clusterVariable(ClusterVariable variable,
0198                                                             const l1t::HGCalMulticluster& cluster) const {
0199   switch (variable) {
0200     case ClusterVariable::cl3d_showerlength:
0201       return cluster.showerLength();
0202     case ClusterVariable::cl3d_coreshowerlength:
0203       return cluster.coreShowerLength();
0204     case ClusterVariable::cl3d_firstlayer:
0205       return cluster.firstLayer();
0206     case ClusterVariable::cl3d_maxlayer:
0207       return cluster.maxLayer();
0208     case ClusterVariable::cl3d_seetot:
0209       return cluster.sigmaEtaEtaTot();
0210     case ClusterVariable::cl3d_seemax:
0211       return cluster.sigmaEtaEtaMax();
0212     case ClusterVariable::cl3d_spptot:
0213       return cluster.sigmaPhiPhiTot();
0214     case ClusterVariable::cl3d_sppmax:
0215       return cluster.sigmaPhiPhiMax();
0216     case ClusterVariable::cl3d_szz:
0217       return cluster.sigmaZZ();
0218     case ClusterVariable::cl3d_srrtot:
0219       return cluster.sigmaRRTot();
0220     case ClusterVariable::cl3d_srrmax:
0221       return cluster.sigmaRRMax();
0222     case ClusterVariable::cl3d_srrmean:
0223       return cluster.sigmaRRMean();
0224     default:
0225       break;
0226   }
0227   return 0.;
0228 }