L1TGlobal/src/AXOL1TLCondition.cc

0001 /**
0002  * \class AXOL1TLCondition
0003  *
0004  *
0005  * Description: evaluation of a condition for axol1tl anomaly detection algorithm
0006  *
0007  * Author: Melissa Quinnan
0008  *
0009  **/
0010
0011 // this class header
0012 #include "L1Trigger/L1TGlobal/interface/CorrCondition.h"
0013
0014 // system include files
0015 #include <iostream>
0016 #include <iomanip>
0017
0018 #include <string>
0019 #include <vector>
0020 #include <algorithm>
0021 #include "ap_fixed.h"
0022
0023 // user include files
0024 //   base classes
0025 #include "L1Trigger/L1TGlobal/interface/AXOL1TLTemplate.h"
0026 #include "L1Trigger/L1TGlobal/interface/ConditionEvaluation.h"
0027
0028 #include "L1Trigger/L1TGlobal/interface/MuCondition.h"
0029 #include "L1Trigger/L1TGlobal/interface/AXOL1TLCondition.h"
0030 #include "L1Trigger/L1TGlobal/interface/CaloCondition.h"
0031 #include "L1Trigger/L1TGlobal/interface/EnergySumCondition.h"
0032 #include "L1Trigger/L1TGlobal/interface/MuonTemplate.h"
0033 #include "L1Trigger/L1TGlobal/interface/CaloTemplate.h"
0034 #include "L1Trigger/L1TGlobal/interface/EnergySumTemplate.h"
0035 #include "L1Trigger/L1TGlobal/interface/GlobalScales.h"
0036
0037 #include "DataFormats/L1Trigger/interface/L1Candidate.h"
0038
0039 #include "L1Trigger/L1TGlobal/interface/GlobalBoard.h"
0040
0041 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0042 #include "FWCore/MessageLogger/interface/MessageDrop.h"
0043
0044 namespace {
0045   //template function for reading results
0046   template <typename ResultType, typename LossType>
0047   LossType readResult(hls4mlEmulator::Model& model) {
0048     std::pair<ResultType, LossType> ADModelResult;  //model outputs a pair of the (result vector, loss)
0049     model.read_result(&ADModelResult);
0050     return ADModelResult.second;
0051   }
0052 }  // namespace
0053
0054 l1t::AXOL1TLCondition::AXOL1TLCondition()
0055     : ConditionEvaluation(), m_gtAXOL1TLTemplate{nullptr}, m_gtGTB{nullptr}, m_model{nullptr} {}
0056
0057 l1t::AXOL1TLCondition::AXOL1TLCondition(const GlobalCondition* axol1tlTemplate, const GlobalBoard* ptrGTB)
0058     : ConditionEvaluation(),
0059       m_gtAXOL1TLTemplate(static_cast<const AXOL1TLTemplate*>(axol1tlTemplate)),
0060       m_gtGTB(ptrGTB),
0061       m_model_loader{kModelNamePrefix + m_gtAXOL1TLTemplate->modelVersion()} {
0062   loadModel();
0063 }
0064
0065 // copy constructor
0066 void l1t::AXOL1TLCondition::copy(const l1t::AXOL1TLCondition& cp) {
0067   m_gtAXOL1TLTemplate = cp.gtAXOL1TLTemplate();
0068   m_gtGTB = cp.gtGTB();
0069
0070   m_condMaxNumberObjects = cp.condMaxNumberObjects();
0071   m_condLastResult = cp.condLastResult();
0072   m_combinationsInCond = cp.getCombinationsInCond();
0073
0074   m_verbosity = cp.m_verbosity;
0075
0076   m_model_loader.reset(cp.model_loader().model_name());
0077   loadModel();
0078 }
0079
0080 l1t::AXOL1TLCondition::AXOL1TLCondition(const l1t::AXOL1TLCondition& cp) : ConditionEvaluation() { copy(cp); }
0081
0082 // destructor
0083 l1t::AXOL1TLCondition::~AXOL1TLCondition() {
0084   // empty
0085 }
0086
0087 // equal operator
0088 l1t::AXOL1TLCondition& l1t::AXOL1TLCondition::operator=(const l1t::AXOL1TLCondition& cp) {
0089   copy(cp);
0090   return *this;
0091 }
0092
0093 // methods
0094 void l1t::AXOL1TLCondition::setGtAXOL1TLTemplate(const AXOL1TLTemplate* caloTempl) { m_gtAXOL1TLTemplate = caloTempl; }
0095
0096 ///   set the pointer to uGT GlobalBoard
0097 void l1t::AXOL1TLCondition::setuGtB(const GlobalBoard* ptrGTB) { m_gtGTB = ptrGTB; }
0098
0099 /// set score for score saving
0100 void l1t::AXOL1TLCondition::setScore(const float scoreval) const { m_savedscore = scoreval; }
0101
0102 void l1t::AXOL1TLCondition::loadModel() {
0103   try {
0104     m_model = m_model_loader.load_model();
0105   } catch (std::runtime_error& e) {
0106     throw cms::Exception("ModelError") << " ERROR: failed to load AXOL1TL model version \""
0107                                        << m_model_loader.model_name()
0108                                        << "\". Model version not found in cms-hls4ml externals.";
0109   }
0110 }
0111
0112 const bool l1t::AXOL1TLCondition::evaluateCondition(const int bxEval) const {
0113   if (m_model == nullptr) {
0114     throw cms::Exception("ModelError") << " ERROR: no model was loaded for AXOL1TL model version \""
0115                                        << m_model_loader.model_name() << "\".";
0116   }
0117
0118   bool condResult = false;
0119   int useBx = bxEval + m_gtAXOL1TLTemplate->condRelativeBx();
0120
0121   // //pointers to objects
0122   const BXVector<const l1t::Muon*>* candMuVec = m_gtGTB->getCandL1Mu();
0123   const BXVector<const l1t::L1Candidate*>* candJetVec = m_gtGTB->getCandL1Jet();
0124   const BXVector<const l1t::L1Candidate*>* candEGVec = m_gtGTB->getCandL1EG();
0125   const BXVector<const l1t::EtSum*>* candEtSumVec = m_gtGTB->getCandL1EtSum();
0126
0127   const int NMuons = 4;
0128   const int NJets = 10;
0129   const int NEgammas = 4;
0130   //const int NEtSums = 1;
0131
0132   //number of indices in vector is #objects * 3 for et, eta, phi
0133   const int MuVecSize = 12;    //NMuons * 3;      //so 12
0134   const int JVecSize = 30;     //NJets * 3;        //so 30
0135   const int EGVecSize = 12;    //NEgammas * 3;    //so 12
0136   const int EtSumVecSize = 3;  //NEtSums * 3;    //so 3
0137
0138   //total # inputs in vector is (4+10+4+1)*3 = 57
0139   const int NInputs = 57;
0140
0141   //types of inputs and outputs
0142   typedef ap_fixed<18, 13> inputtype;
0143   typedef ap_ufixed<18, 14> losstype;
0144
0145   //define zero
0146   inputtype fillzero = 0.0;
0147
0148   //AD vector declaration, will fill later
0149   inputtype ADModelInput[NInputs] = {};
0150
0151   //initializing vector by type for my sanity
0152   inputtype MuInput[MuVecSize];
0153   inputtype JetInput[JVecSize];
0154   inputtype EgammaInput[EGVecSize];
0155   inputtype EtSumInput[EtSumVecSize];
0156
0157   //declare result vectors +score
0158   // resulttype result;
0159   losstype loss;
0160   // pairtype ADModelResult;  //model outputs a pair of the (result vector, loss)
0161   float score = -1.0;  //not sure what the best default is hm??
0162
0163   //check number of input objects we actually have (muons, jets etc)
0164   int NCandMu = candMuVec->size(useBx);
0165   int NCandJet = candJetVec->size(useBx);
0166   int NCandEG = candEGVec->size(useBx);
0167   int NCandEtSum = candEtSumVec->size(useBx);
0168
0169   //initialize arrays to zero (std::fill(first, last, value);)
0170   std::fill(EtSumInput, EtSumInput + EtSumVecSize, fillzero);
0171   std::fill(MuInput, MuInput + MuVecSize, fillzero);
0172   std::fill(JetInput, JetInput + JVecSize, fillzero);
0173   std::fill(EgammaInput, EgammaInput + EGVecSize, fillzero);
0174   std::fill(ADModelInput, ADModelInput + NInputs, fillzero);
0175
0176   //then fill the object vectors
0177   //NOTE assume candidates are already sorted by pt
0178   //loop over EtSums first, easy because there is max 1 of them
0179   if (NCandEtSum > 0) {  //check if not empty
0180     for (int iEtSum = 0; iEtSum < NCandEtSum; iEtSum++) {
0181       if ((candEtSumVec->at(useBx, iEtSum))->getType() == l1t::EtSum::EtSumType::kMissingEt) {
0182         EtSumInput[0] =
0183             ((candEtSumVec->at(useBx, iEtSum))->hwPt()) / 2;  //have to do hwPt/2 in order to match original et inputs
0184         // EtSumInput[1] = (candEtSumVec->at(useBx, iEtSum))->hwEta(); //this one is zero, so leave it zero
0185         EtSumInput[2] = (candEtSumVec->at(useBx, iEtSum))->hwPhi();
0186       }
0187     }
0188   }
0189
0190   //next egammas
0191   if (NCandEG > 0) {  //check if not empty
0192     for (int iEG = 0; iEG < NCandEG; iEG++) {
0193       if (iEG < NEgammas) {  //stop if fill the Nobjects we need
0194         EgammaInput[0 + (3 * iEG)] = ((candEGVec->at(useBx, iEG))->hwPt()) /
0195                                      2;  //index 0,3,6,9 //have to do hwPt/2 in order to match original et inputs
0196         EgammaInput[1 + (3 * iEG)] = (candEGVec->at(useBx, iEG))->hwEta();  //index 1,4,7,10
0197         EgammaInput[2 + (3 * iEG)] = (candEGVec->at(useBx, iEG))->hwPhi();  //index 2,5,8,11
0198       }
0199     }
0200   }
0201
0202   //next muons
0203   if (NCandMu > 0) {  //check if not empty
0204     for (int iMu = 0; iMu < NCandMu; iMu++) {
0205       if (iMu < NMuons) {  //stop if fill the Nobjects we need
0206         MuInput[0 + (3 * iMu)] = ((candMuVec->at(useBx, iMu))->hwPt()) /
0207                                  2;  //index 0,3,6,9 //have to do hwPt/2 in order to match original et inputs
0208         MuInput[1 + (3 * iMu)] = (candMuVec->at(useBx, iMu))->hwEtaAtVtx();  //index 1,4,7,10
0209         MuInput[2 + (3 * iMu)] = (candMuVec->at(useBx, iMu))->hwPhiAtVtx();  //index 2,5,8,11
0210       }
0211     }
0212   }
0213
0214   //next jets
0215   if (NCandJet > 0) {  //check if not empty
0216     for (int iJet = 0; iJet < NCandJet; iJet++) {
0217       if (iJet < NJets) {  //stop if fill the Nobjects we need
0218         JetInput[0 + (3 * iJet)] = ((candJetVec->at(useBx, iJet))->hwPt()) /
0219                                    2;  //index 0,3,6,9...27 //have to do hwPt/2 in order to match original et inputs
0220         JetInput[1 + (3 * iJet)] = (candJetVec->at(useBx, iJet))->hwEta();  //index 1,4,7,10...28
0221         JetInput[2 + (3 * iJet)] = (candJetVec->at(useBx, iJet))->hwPhi();  //index 2,5,8,11...29
0222       }
0223     }
0224   }
0225
0226   //now put it all together-> EtSum+EGamma+Muon+Jet into ADModelInput
0227   int index = 0;
0228   for (int idET = 0; idET < EtSumVecSize; idET++) {
0229     ADModelInput[index++] = EtSumInput[idET];
0230   }
0231   for (int idEG = 0; idEG < EGVecSize; idEG++) {
0232     ADModelInput[index++] = EgammaInput[idEG];
0233   }
0234   for (int idMu = 0; idMu < MuVecSize; idMu++) {
0235     ADModelInput[index++] = MuInput[idMu];
0236   }
0237   for (int idJ = 0; idJ < JVecSize; idJ++) {
0238     ADModelInput[index++] = JetInput[idJ];
0239   }
0240
0241   //now run the inference
0242   m_model->prepare_input(ADModelInput);  //scaling internal here
0243   m_model->predict();
0244   // m_model->read_result(&ADModelResult);  // this should be the square sum model result
0245   if ((m_model_loader.model_name() == "GTADModel_v3") ||
0246       (m_model_loader.model_name() == "GTADModel_v4")) {  //v3/v4 overwrite
0247     using resulttype = std::array<ap_fixed<10, 7, AP_RND_CONV, AP_SAT>, 8>;
0248     loss = readResult<resulttype, losstype>(*m_model);
0249   } else {  //v5 default
0250     using resulttype = ap_fixed<18, 14, AP_RND_CONV, AP_SAT>;
0251     loss = readResult<resulttype, losstype>(*m_model);
0252   }
0253
0254   // result = ADModelResult.first;
0255   // loss = ADModelResult.second;
0256   score = ((loss).to_float()) * 16.0;  //scaling to match threshold
0257   //save score to class variable in case score saving needed
0258   setScore(score);
0259
0260   //number of objects/thrsholds to check
0261   int iCondition = 0;  // number of conditions: there is only one
0262   int nObjInCond = m_gtAXOL1TLTemplate->nrObjects();
0263
0264   if (iCondition >= nObjInCond || iCondition < 0) {
0265     return false;
0266   }
0267
0268   const AXOL1TLTemplate::ObjectParameter objPar = (*(m_gtAXOL1TLTemplate->objectParameter()))[iCondition];
0269
0270   // condGEqVal indicates the operator used for the condition (>=, =): true for >=
0271   bool condGEqVal = m_gtAXOL1TLTemplate->condGEq();
0272   bool passCondition = false;
0273
0274   passCondition = checkCut(objPar.minAXOL1TLThreshold, score, condGEqVal);
0275
0276   condResult |= passCondition;  //condresult true if passCondition true else it is false
0277
0278   //return result
0279   return condResult;
0280 }
0281
0282 void l1t::AXOL1TLCondition::print(std::ostream& myCout) const {
0283   myCout << "Dummy Print for AXOL1TLCondition" << std::endl;
0284   m_gtAXOL1TLTemplate->print(myCout);
0285
0286   ConditionEvaluation::print(myCout);
0287 }