Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-06-20 01:53:35 

0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/hgcalinput_ref.h"
0002 #include "L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/nn_activation.h"
0003 
0004 #ifdef CMSSW_GIT_HASH
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0007 
0008 l1ct::HgcalClusterDecoderEmulator::HgcalClusterDecoderEmulator(const edm::ParameterSet &pset)
0009     : slim_(pset.getParameter<bool>("slim")),
0010       multiclass_id_(pset.getParameterSet("multiclass_id")),
0011       corrector_(pset.getParameter<std::string>("corrector"),
0012                  pset.getParameter<double>("correctorEmfMax"),
0013                  false,
0014                  pset.getParameter<bool>("emulateCorrections"),
0015                  l1tpf::corrector::EmulationMode::Correction),
0016       emInterpScenario_(setEmInterpScenario(pset.getParameter<std::string>("emInterpScenario"))) {}
0017 
0018 edm::ParameterSetDescription l1ct::HgcalClusterDecoderEmulator::getParameterSetDescription() {
0019   edm::ParameterSetDescription description;
0020   description.add<bool>("slim", false);
0021   description.add<std::string>("corrector", "");
0022   description.add<double>("correctorEmfMax", -1);
0023   description.add<bool>("emulateCorrections", false);
0024   description.add<edm::ParameterSetDescription>("multiclass_id", MultiClassID::getParameterSetDescription());
0025   description.add<std::string>("emInterpScenario", "No");
0026   return description;
0027 }
0028 
0029 l1ct::HgcalClusterDecoderEmulator::MultiClassID::MultiClassID(const edm::ParameterSet &pset)
0030     : l1ct::HgcalClusterDecoderEmulator::MultiClassID::MultiClassID(
0031           pset.getParameter<std::string>("model"),
0032           pset.getParameter<std::vector<double>>("wp_pt"),
0033           pset.getParameter<std::vector<double>>("wp_PU"),
0034           pset.getParameter<std::vector<double>>("wp_Pi"),
0035           pset.getParameter<std::vector<double>>("wp_PFEm"),
0036           pset.getParameter<std::vector<double>>("wp_EgEm"),
0037           pset.getParameter<std::vector<double>>("wp_EgEm_tight")) {}
0038 
0039 edm::ParameterSetDescription l1ct::HgcalClusterDecoderEmulator::MultiClassID::getParameterSetDescription() {
0040   edm::ParameterSetDescription description;
0041   description.add<std::string>("model");
0042   description.add<std::vector<double>>("wp_pt");
0043   description.add<std::vector<double>>("wp_PU");
0044   description.add<std::vector<double>>("wp_Pi");
0045   description.add<std::vector<double>>("wp_EgEm");
0046   description.add<std::vector<double>>("wp_EgEm_tight");
0047   description.add<std::vector<double>>("wp_PFEm");
0048   return description;
0049 }
0050 
0051 #endif
0052 
0053 l1ct::HgcalClusterDecoderEmulator::HgcalClusterDecoderEmulator(const std::string &model,
0054                                                                const std::vector<double> &wp_pt,
0055                                                                const std::vector<double> &wp_PU,
0056                                                                const std::vector<double> &wp_Pi,
0057                                                                const std::vector<double> &wp_PFEm,
0058                                                                const std::vector<double> &wp_EgEm,
0059                                                                const std::vector<double> &wp_EgEm_tight,
0060                                                                bool slim,
0061                                                                const std::string &corrector,
0062                                                                float correctorEmfMax,
0063                                                                bool emulateCorrections,
0064                                                                const std::string &emInterpScenario)
0065     : slim_{slim},
0066       multiclass_id_(model, wp_pt, wp_PU, wp_Pi, wp_PFEm, wp_EgEm, wp_EgEm_tight),
0067       corrector_(corrector, correctorEmfMax, false, emulateCorrections, l1tpf::corrector::EmulationMode::Correction),
0068       emInterpScenario_(setEmInterpScenario(emInterpScenario)) {}
0069 
0070 l1ct::HgcalClusterDecoderEmulator::UseEmInterp l1ct::HgcalClusterDecoderEmulator::setEmInterpScenario(
0071     const std::string &emInterpScenario) {
0072   if (emInterpScenario == "no")
0073     return UseEmInterp::No;
0074   if (emInterpScenario == "emOnly")
0075     return UseEmInterp::EmOnly;
0076   if (emInterpScenario == "allKeepHad")
0077     return UseEmInterp::AllKeepHad;
0078   if (emInterpScenario == "allKeepTot")
0079     return UseEmInterp::AllKeepTot;
0080   throw std::runtime_error("Unknown emInterpScenario: " + emInterpScenario);
0081 }
0082 
0083 l1ct::HadCaloObjEmu l1ct::HgcalClusterDecoderEmulator::decode(const l1ct::PFRegionEmu &sector,
0084                                                               const ap_uint<256> &in,
0085                                                               bool &valid) const {
0086   // Word 0
0087   ap_uint<14> w_pt = in(13, 0);       // 14 bits: 13-0
0088   ap_uint<14> w_empt = in(27, 14);    // 14 bits: 27-14
0089   ap_uint<4> w_gctqual = in(31, 28);  //  4 bits: 31-28
0090   ap_uint<8> w_emf_tot = in(39, 32);  //  8 bits: 39-32
0091   ap_uint<8> w_emf = in(47, 40);      //  8 bits: 47-40
0092 
0093   // Word 1
0094   ap_uint<10> w_abseta = in(64 + 9, 64 + 0);   // 10 bits: 9-0
0095   ap_int<9> w_phi = in(64 + 18, 64 + 10);      //  9 bits: 18-10
0096   ap_uint<12> w_meanz = in(64 + 30, 64 + 19);  // 12 bits: 30-19
0097 
0098   // Word 2
0099   ap_uint<6> w_showerlength = in(128 + 18, 128 + 13);      //  6 bits: 18-13
0100   ap_uint<7> w_sigmazz = in(128 + 38, 128 + 32);           //  7 bits: 38-32
0101   ap_uint<7> w_sigmaphiphi = in(128 + 45, 128 + 39);       //  7 bits: 45-39
0102   ap_uint<6> w_coreshowerlength = in(128 + 51, 128 + 46);  //  6 bits: 51-46
0103   ap_uint<5> w_sigmaetaeta = in(128 + 56, 128 + 52);       //  5 bits: 56-52
0104 
0105   // Word 3
0106   ap_uint<13> w_sigmarrtot = in(213, 201);  // 13 bits: 213-201 // FIXME: use word3 spare bits
0107 
0108   // Conversion to local (input sector) coordinates
0109   ap_int<9> w_eta = l1ct::glbeta_t(w_abseta.to_int() * (sector.floatEtaCenter() > 0 ? +1 : -1)) - sector.hwEtaCenter;
0110 
0111   l1ct::HadCaloObjEmu out;
0112   out.clear();
0113   if (w_pt == 0)
0114     return out;
0115   // if (w_pt == 0 || w_phi > sector.hwPhiHalfWidth || w_phi <= -sector.hwPhiHalfWidth)
0116   //   return out;
0117   out.hwPt = w_pt * l1ct::pt_t(l1ct::Scales::INTPT_LSB);
0118   out.hwEta = w_eta;
0119   out.hwPhi = w_phi;  // relative to the region center, at calo
0120   out.hwEmPt = w_empt * l1ct::pt_t(l1ct::Scales::INTPT_LSB);
0121 
0122   if (!slim_) {
0123     // FIXME: the scaling here is added to the encoded word.
0124     out.hwSrrTot = w_sigmarrtot * l1ct::srrtot_t(l1ct::Scales::SRRTOT_LSB);
0125     // We just downscale precision and round to the nearest integer
0126     out.hwMeanZ = l1ct::meanz_t(std::min(w_meanz.to_int() + 1, (1 << 12) - 1) >> 1);
0127     // Compute an H/E value: 1/emf - 1 as needed by Composite ID
0128     // NOTE: this uses the total cluster energy, which is not the case for the eot shower shape!
0129     // FIXME: could drop once we move the model to the eot fraction
0130     ap_ufixed<10, 5, AP_RND_CONV, AP_SAT> w_hoe = 256.0 / (w_emf_tot.to_int() + 0.5) - 1;
0131     out.hwHoe = w_hoe;
0132   }
0133   std::vector<MultiClassID::bdt_feature_t> inputs = {w_showerlength,
0134                                                      w_coreshowerlength,
0135                                                      w_emf,
0136                                                      w_abseta - 256,
0137                                                      w_meanz,  // We use the full resolution here
0138                                                      w_sigmaetaeta,
0139                                                      w_sigmaphiphi,
0140                                                      w_sigmazz};
0141 
0142   // Apply EM interpretation scenario
0143   if (emInterpScenario_ == UseEmInterp::No) {  // we do not use EM interpretation
0144     out.hwEmPt = w_emf_tot * out.hwPt / 256;
0145     // NOTE: only case where hoe consisten with hwEmPt
0146   } else if (emInterpScenario_ == UseEmInterp::EmOnly) {  // for emID objs, use EM interp as pT and set H = 0
0147     if (out.hwEmID) {
0148       out.hwPt = out.hwEmPt;
0149       out.hwHoe = 0;
0150     }
0151   } else if (emInterpScenario_ ==
0152              UseEmInterp::AllKeepHad) {  // for all objs, replace EM part with EM interp, preserve H
0153     l1ct::pt_t had_pt = out.hwPt - w_emf_tot * out.hwPt / 256;
0154     out.hwPt = had_pt + out.hwEmPt;
0155     // FIXME: we do not recompute hoe for now...
0156   } else if (emInterpScenario_ ==
0157              UseEmInterp::AllKeepTot) {  // for all objs, replace EM part with EM interp, preserve pT
0158     // FIXME: we do not recompute hoe for now...
0159   }
0160 
0161   bool notPU = multiclass_id_.evaluate(out, inputs);
0162 
0163   // Calibrate pt and set error
0164   if (corrector_.valid()) {
0165     float newpt = corrector_.correctedPt(out.floatPt(), out.floatEmPt(), sector.floatGlbEta(out.hwEta));
0166     out.hwPt = l1ct::Scales::makePtFromFloat(newpt);
0167     // NOTE: hoe/emfrac are not updated
0168   }
0169 
0170   // evaluate multiclass model
0171   valid = notPU && out.hwPt > 0;
0172 
0173   return out;
0174 }
0175 
0176 l1ct::HgcalClusterDecoderEmulator::MultiClassID::MultiClassID(const std::string &model,
0177                                                               const std::vector<double> &wp_pt,
0178                                                               const std::vector<double> &wp_PU,
0179                                                               const std::vector<double> &wp_Pi,
0180                                                               const std::vector<double> &wp_PFEm,
0181                                                               const std::vector<double> &wp_EgEm,
0182                                                               const std::vector<double> &wp_EgEm_tight) {
0183   for (auto pt : wp_pt)
0184     wp_pt_.emplace_back(pt);
0185   for (auto pu : wp_PU)
0186     wp_PU_.emplace_back(pu);
0187   for (auto pi : wp_Pi)
0188     wp_Pi_.emplace_back(pi);
0189   for (auto egem : wp_EgEm)
0190     wp_EgEm_.emplace_back(egem);
0191   for (auto pfem : wp_PFEm)
0192     wp_PFEm_.emplace_back(pfem);
0193   for (auto egem : wp_EgEm_tight)
0194     wp_EgEm_tight_.emplace_back(egem);
0195 
0196 #ifdef CMSSW_GIT_HASH
0197   auto resolvedFileName = edm::FileInPath(model).fullPath();
0198 #else
0199   auto resolvedFileName = model;
0200 #endif
0201   multiclass_bdt_ = std::make_unique<conifer::BDT<bdt_feature_t, bdt_score_t, false>>(resolvedFileName);
0202 }
0203 
0204 bool l1ct::HgcalClusterDecoderEmulator::MultiClassID::evaluate(l1ct::HadCaloObjEmu &cl,
0205                                                                const std::vector<bdt_feature_t> &inputs) const {
0206   auto bdt_score = multiclass_bdt_->decision_function(inputs);  //0 is pu, 1 is pi, 2 is eg
0207   bdt_score_t raw_scores[3] = {bdt_score[0], bdt_score[1], bdt_score[2]};
0208   l1ct::id_prob_t sm_scores[3];
0209   nnet::softmax_stable<bdt_score_t, l1ct::id_prob_t, softmax_config>(raw_scores, sm_scores);
0210 
0211   // softmax_stable<>
0212 
0213   unsigned int pt_bin = 0;
0214   for (size_t i = wp_pt_.size(); i > 0; --i) {
0215     if (cl.hwPt >= wp_pt_[i - 1]) {
0216       pt_bin = i;
0217       break;
0218     }
0219   }
0220   bool passPu = (sm_scores[0] >= wp_PU_[pt_bin]);
0221   // bool passPi = (sm_scores[1] >= wp_Pi_[pt_bin]);  // FIXME: where do we store this?
0222   bool passPFEm = (sm_scores[2] >= wp_PFEm_[pt_bin]);
0223   bool passEgEm = (sm_scores[2] >= wp_EgEm_[pt_bin]);
0224   bool passEgEm_tight = (sm_scores[2] >= wp_EgEm_tight_[pt_bin]);
0225 
0226   // bit 0: PF EM ID
0227   // bit 1: EG EM ID
0228   // bit 2: EG Loose ID
0229   cl.hwEmID = passPFEm | (passEgEm_tight << 1) | (passEgEm << 2);
0230 
0231   cl.hwPiProb = sm_scores[1];
0232   cl.hwEmProb = sm_scores[2];
0233   return !passPu;
0234 }
0235 
0236 void l1ct::HgcalClusterDecoderEmulator::MultiClassID::softmax(const float rawScores[3], float scores[3]) const {
0237   // softmax (for now, let's compute the softmax in this code; this needs to be changed to implement on firmware)
0238   // Softmax implemented in conifer (standalone) is to be integrated here soon; for now, just do "offline" softmax :(
0239   float denom = exp(rawScores[0]) + exp(rawScores[1]) + exp(rawScores[2]);
0240   scores[0] = exp(rawScores[0]) / denom;
0241   scores[1] = exp(rawScores[1]) / denom;
0242   scores[2] = exp(rawScores[2]) / denom;
0243 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!