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