Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-02-03 10:47:56

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 const bool l1t::AXOL1TLCondition::evaluateCondition(const int bxEval) const {
0089   bool condResult = false;
0090   int useBx = bxEval + m_gtAXOL1TLTemplate->condRelativeBx();
0091 
0092   //HLS4ML stuff
0093   std::string AXOL1TLmodelversion = m_AXOL1TLmodelversion;  //config loading method
0094   hls4mlEmulator::ModelLoader loader(AXOL1TLmodelversion);
0095   std::shared_ptr<hls4mlEmulator::Model> model;
0096   model = loader.load_model();
0097 
0098   // //pointers to objects
0099   const BXVector<const l1t::Muon*>* candMuVec = m_gtGTB->getCandL1Mu();
0100   const BXVector<const l1t::L1Candidate*>* candJetVec = m_gtGTB->getCandL1Jet();
0101   const BXVector<const l1t::L1Candidate*>* candEGVec = m_gtGTB->getCandL1EG();
0102   const BXVector<const l1t::EtSum*>* candEtSumVec = m_gtGTB->getCandL1EtSum();
0103 
0104   const int NMuons = 4;
0105   const int NJets = 10;
0106   const int NEgammas = 4;
0107   //const int NEtSums = 1;
0108 
0109   //number of indices in vector is #objects * 3 for et, eta, phi
0110   const int MuVecSize = 12;    //NMuons * 3;      //so 12
0111   const int JVecSize = 30;     //NJets * 3;        //so 30
0112   const int EGVecSize = 12;    //NEgammas * 3;    //so 12
0113   const int EtSumVecSize = 3;  //NEtSums * 3;    //so 3
0114 
0115   //total # inputs in vector is (4+10+4+1)*3 = 57
0116   const int NInputs = 57;
0117 
0118   //types of inputs and outputs
0119   typedef ap_fixed<18, 13> inputtype;
0120   typedef std::array<ap_fixed<10, 7, AP_RND_CONV, AP_SAT>, 8> resulttype;  //v3
0121   typedef ap_ufixed<18, 14> losstype;
0122   typedef std::pair<resulttype, losstype> pairtype;
0123   // typedef std::array<ap_fixed<10, 7>, 13> resulttype;  //deprecated v1 type:
0124 
0125   //define zero
0126   inputtype fillzero = 0.0;
0127 
0128   //AD vector declaration, will fill later
0129   inputtype ADModelInput[NInputs] = {};
0130 
0131   //initializing vector by type for my sanity
0132   inputtype MuInput[MuVecSize];
0133   inputtype JetInput[JVecSize];
0134   inputtype EgammaInput[EGVecSize];
0135   inputtype EtSumInput[EtSumVecSize];
0136 
0137   //declare result vectors +score
0138   resulttype result;
0139   losstype loss;
0140   pairtype ADModelResult;  //model outputs a pair of the (result vector, loss)
0141   float score = -1.0;      //not sure what the best default is hm??
0142 
0143   //check number of input objects we actually have (muons, jets etc)
0144   int NCandMu = candMuVec->size(useBx);
0145   int NCandJet = candJetVec->size(useBx);
0146   int NCandEG = candEGVec->size(useBx);
0147   int NCandEtSum = candEtSumVec->size(useBx);
0148 
0149   //initialize arrays to zero (std::fill(first, last, value);)
0150   std::fill(EtSumInput, EtSumInput + EtSumVecSize, fillzero);
0151   std::fill(MuInput, MuInput + MuVecSize, fillzero);
0152   std::fill(JetInput, JetInput + JVecSize, fillzero);
0153   std::fill(EgammaInput, EgammaInput + EGVecSize, fillzero);
0154   std::fill(ADModelInput, ADModelInput + NInputs, fillzero);
0155 
0156   //then fill the object vectors
0157   //NOTE assume candidates are already sorted by pt
0158   //loop over EtSums first, easy because there is max 1 of them
0159   if (NCandEtSum > 0) {  //check if not empty
0160     for (int iEtSum = 0; iEtSum < NCandEtSum; iEtSum++) {
0161       if ((candEtSumVec->at(useBx, iEtSum))->getType() == l1t::EtSum::EtSumType::kMissingEt) {
0162         EtSumInput[0] =
0163             ((candEtSumVec->at(useBx, iEtSum))->hwPt()) / 2;  //have to do hwPt/2 in order to match original et inputs
0164         // EtSumInput[1] = (candEtSumVec->at(useBx, iEtSum))->hwEta(); //this one is zero, so leave it zero
0165         EtSumInput[2] = (candEtSumVec->at(useBx, iEtSum))->hwPhi();
0166       }
0167     }
0168   }
0169 
0170   //next egammas
0171   if (NCandEG > 0) {  //check if not empty
0172     for (int iEG = 0; iEG < NCandEG; iEG++) {
0173       if (iEG < NEgammas) {  //stop if fill the Nobjects we need
0174         EgammaInput[0 + (3 * iEG)] = ((candEGVec->at(useBx, iEG))->hwPt()) /
0175                                      2;  //index 0,3,6,9 //have to do hwPt/2 in order to match original et inputs
0176         EgammaInput[1 + (3 * iEG)] = (candEGVec->at(useBx, iEG))->hwEta();  //index 1,4,7,10
0177         EgammaInput[2 + (3 * iEG)] = (candEGVec->at(useBx, iEG))->hwPhi();  //index 2,5,8,11
0178       }
0179     }
0180   }
0181 
0182   //next muons
0183   if (NCandMu > 0) {  //check if not empty
0184     for (int iMu = 0; iMu < NCandMu; iMu++) {
0185       if (iMu < NMuons) {  //stop if fill the Nobjects we need
0186         MuInput[0 + (3 * iMu)] = ((candMuVec->at(useBx, iMu))->hwPt()) /
0187                                  2;  //index 0,3,6,9 //have to do hwPt/2 in order to match original et inputs
0188         MuInput[1 + (3 * iMu)] = (candMuVec->at(useBx, iMu))->hwEta();  //index 1,4,7,10
0189         MuInput[2 + (3 * iMu)] = (candMuVec->at(useBx, iMu))->hwPhi();  //index 2,5,8,11
0190       }
0191     }
0192   }
0193 
0194   //next jets
0195   if (NCandJet > 0) {  //check if not empty
0196     for (int iJet = 0; iJet < NCandJet; iJet++) {
0197       if (iJet < NJets) {  //stop if fill the Nobjects we need
0198         JetInput[0 + (3 * iJet)] = ((candJetVec->at(useBx, iJet))->hwPt()) /
0199                                    2;  //index 0,3,6,9...27 //have to do hwPt/2 in order to match original et inputs
0200         JetInput[1 + (3 * iJet)] = (candJetVec->at(useBx, iJet))->hwEta();  //index 1,4,7,10...28
0201         JetInput[2 + (3 * iJet)] = (candJetVec->at(useBx, iJet))->hwPhi();  //index 2,5,8,11...29
0202       }
0203     }
0204   }
0205 
0206   //now put it all together-> EtSum+EGamma+Muon+Jet into ADModelInput
0207   int index = 0;
0208   for (int idET = 0; idET < EtSumVecSize; idET++) {
0209     ADModelInput[index++] = EtSumInput[idET];
0210   }
0211   for (int idEG = 0; idEG < EGVecSize; idEG++) {
0212     ADModelInput[index++] = EgammaInput[idEG];
0213   }
0214   for (int idMu = 0; idMu < MuVecSize; idMu++) {
0215     ADModelInput[index++] = MuInput[idMu];
0216   }
0217   for (int idJ = 0; idJ < JVecSize; idJ++) {
0218     ADModelInput[index++] = JetInput[idJ];
0219   }
0220 
0221   //now run the inference
0222   model->prepare_input(ADModelInput);  //scaling internal here
0223   model->predict();
0224   model->read_result(&ADModelResult);  // this should be the square sum model result
0225 
0226   result = ADModelResult.first;
0227   loss = ADModelResult.second;
0228   score = ((loss).to_float()) * 16.0;  //scaling to match threshold
0229 
0230   //number of objects/thrsholds to check
0231   int iCondition = 0;  // number of conditions: there is only one
0232   int nObjInCond = m_gtAXOL1TLTemplate->nrObjects();
0233 
0234   if (iCondition >= nObjInCond || iCondition < 0) {
0235     return false;
0236   }
0237 
0238   const AXOL1TLTemplate::ObjectParameter objPar = (*(m_gtAXOL1TLTemplate->objectParameter()))[iCondition];
0239 
0240   // condGEqVal indicates the operator used for the condition (>=, =): true for >=
0241   bool condGEqVal = m_gtAXOL1TLTemplate->condGEq();
0242   bool passCondition = false;
0243 
0244   passCondition = checkCut(objPar.minAXOL1TLThreshold, score, condGEqVal);
0245 
0246   condResult |= passCondition;  //condresult true if passCondition true else it is false
0247 
0248   //return result
0249   return condResult;
0250 }
0251 
0252 //in order to set model version from config
0253 void l1t::AXOL1TLCondition::setModelVersion(const std::string modelversionname) {
0254   m_AXOL1TLmodelversion = modelversionname;
0255 }
0256 
0257 void l1t::AXOL1TLCondition::print(std::ostream& myCout) const {
0258   myCout << "Dummy Print for AXOL1TLCondition" << std::endl;
0259   m_gtAXOL1TLTemplate->print(myCout);
0260 
0261   ConditionEvaluation::print(myCout);
0262 }