File indexing completed on 2024-04-06 12:20:39
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
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 }
