Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-11 03:34:26 

0001 /* -*- C++ -*-
0002 
0003 Package: L1CaloTrigger
0004 Class: L1NNCaloTauEmulator
0005 Frinedly name: The TauMinator
0006 
0007 \class L1NNCaloTauEmulator L1NNCaloTauEmulator.cc
0008 
0009 Description: 
0010 Perform firmware-exact emulation of the l1tNNCaloTauProducer
0011 that implements the NN Calo Tau.
0012 (Perform reconstruction and identification of tau 
0013 candidates at L1 Trigger with a CNN.)
0014 
0015 Implementation:
0016 The implementation is done forseeing the integration
0017 of the algorithm in the GCT Sum card. This means that
0018 the full detector information can be accessed at the same
0019 time (ECAL, HCAL, HGCAL full eta-phi coverage). This will
0020 come in the form of arrays of towers and clusters.
0021 Given that the emulators of the upstream algortihms are
0022 not fully determined yet, this emulator takes as input
0023 the simulation-based information, manipulates with software
0024 precision to pruduce the arrays of towers and clusters as 
0025 they should be available in the GCT sum card.
0026 Only then the actual fixed point algorithm emulation arrives.
0027 
0028 ** INFO : THE NNs ARE APPLIED USING THE TENSORFLOW SOFTWARE
0029           the implementation of full emulation via hls4ml is ongoing
0030           (it has already been shown in other contexts that tensorflow 
0031           softwrae and full emulation are very close to each other)
0032 
0033 Original Author: Jona Motta
0034 Created: Tue June 7th 2023
0035 
0036 */
0037 
0038 #include <iostream>
0039 #include <vector>
0040 #include <cmath>
0041 
0042 #include "boost/property_tree/ptree.hpp"
0043 #include "boost/property_tree/json_parser.hpp"
0044 
0045 #include "ap_int.h"
0046 #include "ap_fixed.h"
0047 // #include "hls4ml/emulator.h"
0048 
0049 #include "FWCore/Framework/interface/Frameworkfwd.h"
0050 #include "FWCore/Framework/interface/stream/EDProducer.h"
0051 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0052 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0053 #include "FWCore/ServiceRegistry/interface/Service.h"
0054 #include "FWCore/Framework/interface/Event.h"
0055 #include "FWCore/Framework/interface/ESHandle.h"
0056 #include "FWCore/Framework/interface/MakerMacros.h"
0057 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0058 
0059 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h"
0060 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0061 #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
0062 #include "DataFormats/L1TParticleFlow/interface/PFCluster.h"
0063 #include "DataFormats/L1THGCal/interface/HGCalTower.h"
0064 #include "DataFormats/Math/interface/deltaPhi.h"
0065 #include "DataFormats/L1Trigger/interface/Tau.h"
0066 
0067 #include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h"
0068 #include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h"
0069 
0070 #include "L1Trigger/L1THGCal/interface/backend/HGCalTriggerClusterIdentificationBase.h"
0071 #include "L1Trigger/Phase2L1ParticleFlow/interface/HGC3DClusterEgID.h"
0072 #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
0073 
0074 #include "PhysicsTools/TensorFlow/interface/TensorFlow.h"
0075 
0076 struct NNmodels_GlobalCache {
0077   std::string CNNmodel_CB_path;
0078   std::string DNNident_CB_path;
0079   std::string DNNcalib_CB_path;
0080 
0081   std::string CNNmodel_CE_path;
0082   std::string DNNident_CE_path;
0083   std::string DNNcalib_CE_path;
0084   std::string FeatScaler_CE_path;
0085   boost::property_tree::ptree FeatScaler_CE;
0086 
0087   tensorflow::GraphDef* CNNmodel_CB;
0088   tensorflow::GraphDef* DNNident_CB;
0089   tensorflow::GraphDef* DNNcalib_CB;
0090 
0091   tensorflow::Session* CNNmodel_CBsession;
0092   tensorflow::Session* DNNident_CBsession;
0093   tensorflow::Session* DNNcalib_CBsession;
0094 
0095   tensorflow::GraphDef* CNNmodel_CE;
0096   tensorflow::GraphDef* DNNident_CE;
0097   tensorflow::GraphDef* DNNcalib_CE;
0098 
0099   tensorflow::Session* CNNmodel_CEsession;
0100   tensorflow::Session* DNNident_CEsession;
0101   tensorflow::Session* DNNcalib_CEsession;
0102 };
0103 
0104 class L1NNCaloTauEmulator : public edm::stream::EDProducer<edm::GlobalCache<NNmodels_GlobalCache>> {
0105 public:
0106   explicit L1NNCaloTauEmulator(const edm::ParameterSet&, const NNmodels_GlobalCache*);
0107 
0108   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0109   static std::unique_ptr<NNmodels_GlobalCache> initializeGlobalCache(const edm::ParameterSet&);
0110   static void globalEndJob(const NNmodels_GlobalCache*){/*do nothing*/};
0111 
0112 private:
0113   // ----fixed LSBs, Nbits, scales, and types----
0114   static constexpr int INTPHI_PI = 36;
0115   static constexpr int INTPHI_2PI = 2 * INTPHI_PI;
0116   static constexpr float IETAPHI_LSB = M_PI / INTPHI_PI;
0117 
0118   static constexpr int FINEINTPHI_PI = 720;
0119   static constexpr int FINEINTPHI_2PI = 2 * FINEINTPHI_PI;
0120   static constexpr float ETAPHI_LSB = M_PI / FINEINTPHI_PI;
0121 
0122   static constexpr float SHAPEFEAT_LSB = 0.0000153;  // pow(2, -16)
0123   static constexpr float SZZ_LSB = SHAPEFEAT_LSB * 100;
0124   static constexpr float ETAHGCAL_OFFSET = 1.321;  // inferred from hgcal schematics
0125   static constexpr float IETAHGCAL_LSBp = 0.0808;  // inferred from simulation
0126   static constexpr float IETAHGCAL_LSB = 0.0845;   // inferred from simulation
0127   static constexpr float PUID_LSB = 0.00390625;    // pow(2, -8)
0128   static constexpr float MEANZ_OFFSET = 321.05;    // inferred from hgcal schematics
0129   static constexpr int IETAHGCAL_OFFSET = 17;
0130   static constexpr float MEANZ_LSB = 0.5;
0131   static constexpr float PTET_LSB = 0.25;
0132   static constexpr float CM2MM = 10;
0133   static constexpr int R2cone = 0.25 / ETAPHI_LSB / ETAPHI_LSB;
0134 
0135   static constexpr int SHAPEFEAT_W = 16;  // maximum forseen per shape
0136   static constexpr int DETAPHI_W = 12;
0137   static constexpr int DIETAPHI_W = 8;
0138   static constexpr int IETAPHI_W = 7;
0139   static constexpr int SHOWLEN_W = 6;
0140   static constexpr int ETAPHI_W = 11;  // precision going to correlator
0141   static constexpr int MEANZ_W = 12;
0142   static constexpr int PUID_W = 9;
0143 
0144   static constexpr int PT_W = 14;
0145   static constexpr int PT_I = 12;
0146   static constexpr int ET_W = 10;
0147   static constexpr int ET_I = 8;
0148 
0149   // forseen precision of the HLS4ML emulation of the NNs
0150   static constexpr int CALIBPT_W = 10;
0151   static constexpr int CALIBPT_I = 9;
0152   static constexpr int ID_W = 8;
0153   static constexpr int ID_I = 1;
0154 
0155   typedef ap_ufixed<PT_W, PT_I, AP_TRN, AP_SAT> Pt_t;
0156   typedef ap_ufixed<ET_W, ET_I, AP_TRN, AP_SAT> Et_t;
0157 
0158   typedef ap_ufixed<CALIBPT_W, CALIBPT_I, AP_TRN, AP_SAT> CalibPt_t;
0159   typedef ap_ufixed<ID_W, ID_I, AP_TRN, AP_SAT> Id_t;
0160 
0161   typedef ap_uint<SHAPEFEAT_W> ShapeFeat_t;
0162   typedef ap_int<DIETAPHI_W> dIEtaPhi_t;
0163   typedef ap_int<DETAPHI_W> dEtaPhi_t;
0164   typedef ap_uint<SHOWLEN_W> ShowLen_t;
0165   typedef ap_int<ETAPHI_W> EtaPhi_t;
0166   typedef ap_uint<IETAPHI_W> IPhi_t;
0167   typedef ap_int<IETAPHI_W> IEta_t;
0168   typedef ap_uint<MEANZ_W> Meanz_t;
0169   typedef ap_int<PUID_W> PUid_t;
0170 
0171   // ----fixed dimensions of tower clusters----
0172   const int seedIdx = 22;
0173   const int IEta_dim = 5;
0174   const int IPhi_dim = 9;
0175   const int Eta_limit = 33;
0176 
0177   //----edm control---
0178   void produce(edm::Event&, const edm::EventSetup&) override;
0179 
0180   //----private functions----
0181   template <class outPrecision, class inPrecision>
0182   outPrecision dPhi(inPrecision iPhi_1, inPrecision iPhi_2);
0183   dIEtaPhi_t tower_dIEta(IEta_t iEta_1, IEta_t iEta_2);
0184   dEtaPhi_t tw2cl_dPhi(EtaPhi_t iPhi_1, IPhi_t iPhi_2);
0185   dEtaPhi_t tw2cl_dEta(EtaPhi_t iEta_1, IEta_t iEta_2);
0186   IEta_t makeEndcapHwIEta(float eta);
0187   IPhi_t makeEndcapHwIPhi(float phi);
0188   float apfixedQuantizer(float inputF, float LSB, int nbits);
0189   int apintQuantizer(float inputF, float LSB, int nbits);
0190   float inputScaler(float inputF, std::string feature);
0191   float correctInputEtaCl3d(float eta);
0192   float correctInputMeanzCl3d(float meanz);
0193 
0194   inline float floatPt(Pt_t pt) { return pt.to_float(); }
0195   inline float floatEt(Et_t et) { return et.to_float(); }
0196   inline float floatEta(EtaPhi_t eta) { return eta.to_float() * ETAPHI_LSB; }
0197   inline float floatPhi(EtaPhi_t phi) { return phi.to_float() * ETAPHI_LSB; }
0198   inline float floatShape(ShapeFeat_t shape) { return shape.to_float() * SHAPEFEAT_LSB; };
0199   inline float floatSzz(ShapeFeat_t szz) { return szz.to_float() * SZZ_LSB; };
0200   inline float floatMeanZ(Meanz_t meanz) { return meanz.to_float() * MEANZ_LSB + MEANZ_OFFSET; };
0201   inline float floatMeanZHgcalCoord(Meanz_t meanz) { return meanz.to_float() * MEANZ_LSB; };
0202   inline float floatPuId(PUid_t pu) { return pu.to_float() * PUID_LSB; };
0203   float floatIEta(IEta_t eta);
0204   float floatIPhi(IPhi_t phi);
0205 
0206   template <int W>
0207   ap_int<W> ap_abs(ap_int<W> x);
0208   template <int W, int I, ap_q_mode _AP_Q, ap_o_mode _AP_O>
0209   ap_ufixed<W, I> ap_abs(ap_fixed<W, I, _AP_Q, _AP_O> x);
0210 
0211   //----tokens and handles----
0212   edm::EDGetTokenT<l1tp2::CaloTowerCollection> l1TowersToken;
0213   edm::Handle<l1tp2::CaloTowerCollection> l1CaloTowerHandle;
0214 
0215   edm::EDGetToken hgcalTowersToken;
0216   edm::Handle<l1t::HGCalTowerBxCollection> hgcalTowersHandle;
0217 
0218   edm::EDGetTokenT<l1t::HGCalMulticlusterBxCollection> HGClusterToken;
0219   edm::Handle<l1t::HGCalMulticlusterBxCollection> HGClusterHandle;
0220 
0221   //----private variables----
0222   enum class UseEmInterp { No, EmOnly, AllKeepHad, AllKeepTot };
0223   UseEmInterp scenario;
0224   StringCutObjectSelector<l1t::HGCalMulticluster> preEmId;
0225   l1tpf::HGC3DClusterEgID VsPuId;
0226 
0227   double EcalEtMinForClustering;
0228   double HcalEtMinForClustering;
0229   double EtMinForSeeding;
0230   double EtaRestriction;
0231   double CB_CE_split;
0232   double PuidThr;
0233 
0234   double IdWp90_CB;
0235   double IdWp95_CB;
0236   double IdWp99_CB;
0237 
0238   double IdWp90_CE;
0239   double IdWp95_CE;
0240   double IdWp99_CE;
0241 
0242   PUid_t intPuidThr;
0243   IEta_t intEtaRestriction;
0244   IEta_t intCB_CE_split;
0245 
0246   bool DEBUG;
0247 
0248   // Class for the towers info as they should be in GCT
0249   class SimpleTowerHit {
0250   public:
0251     IEta_t towerIeta = 0;
0252     IPhi_t towerIphi = 0;
0253     Et_t towerEm = 0.;
0254     Et_t towerHad = 0.;
0255     Et_t l1egTowerEt = 0.;
0256     Et_t towerEt = 0.;
0257     ap_uint<1> isBarrel = 0x1;
0258     ap_uint<1> stale = 0x0;
0259     ap_uint<1> stale4seed = 0x0;
0260   };
0261 
0262   // Class for the clusters info as they should arrive from HGCAL
0263   class SimpleHGCluster {
0264   public:
0265     Pt_t pt;
0266     EtaPhi_t eta;
0267     EtaPhi_t phi;
0268     ShowLen_t showerlength;
0269     ShowLen_t coreshowerlength;
0270     ShapeFeat_t spptot;
0271     ShapeFeat_t szz;
0272     ShapeFeat_t srrtot;
0273     Meanz_t meanz;
0274     PUid_t PUid;
0275     ap_uint<1> stale = 0x0;
0276   };
0277 
0278   // Classes for the tower clusters
0279   class SimplifiedTower {
0280   public:
0281     Et_t towerEm = 0.;
0282     Et_t towerHad = 0.;
0283     Et_t l1egTowerEt = 0.;
0284 
0285     void fill(SimpleTowerHit Tower) {
0286       towerEm = Tower.towerEm;
0287       towerHad = Tower.towerHad;
0288       l1egTowerEt = Tower.l1egTowerEt;
0289     }
0290   };
0291 
0292   class InputTowerCluster {
0293   public:
0294     SimplifiedTower towerHits[45];
0295     ap_uint<1> barrelSeeded = 0x0;
0296     ap_uint<1> filled[45];
0297 
0298     void fill(int idx, SimpleTowerHit Tower) {
0299       towerHits[idx].fill(Tower);
0300       filled[idx] = 0x1;
0301     }
0302 
0303     void init() {
0304       SimplifiedTower emptyT;
0305       std::fill(towerHits, towerHits + 44, emptyT);
0306       std::fill(filled, filled + 44, 0x0);
0307     }
0308   };
0309 
0310   class InputTowerCluster_pstn {
0311   public:
0312     IEta_t seedIeta = 0;
0313     IPhi_t seedIphi = 0;
0314 
0315     void fill(SimpleTowerHit Tower) {
0316       seedIeta = Tower.towerIeta;
0317       seedIphi = Tower.towerIphi;
0318     }
0319   };
0320 
0321   // INFO : now variables are in GCT precision, they should be in NN precision
0322   // after scaling, i.e. something like ap_fixed<16, 6, AP_TRN, AP_SAT>, when the
0323   // full hls4ml emulation is available
0324   class InputHGCluster {
0325   public:
0326     Pt_t pt;
0327     EtaPhi_t eta;
0328     ShowLen_t showerlength;
0329     ShowLen_t coreshowerlength;
0330     ShapeFeat_t spptot;
0331     ShapeFeat_t szz;
0332     ShapeFeat_t srrtot;
0333     Meanz_t meanz;
0334 
0335     void fill(SimpleHGCluster Cluster) {
0336       pt = Cluster.pt;
0337       eta = Cluster.eta;
0338       showerlength = Cluster.showerlength;
0339       coreshowerlength = Cluster.coreshowerlength;
0340       spptot = Cluster.spptot;
0341       szz = Cluster.szz;
0342       srrtot = Cluster.srrtot;
0343       meanz = Cluster.meanz;
0344     }
0345   };
0346 
0347   l1t::Tau MakeTauCandidate(bool isBarrel,
0348                             int clNxMIdx,
0349                             std::vector<tensorflow::Tensor> outputsIdent,
0350                             std::vector<tensorflow::Tensor> outputsCalib,
0351                             std::vector<InputTowerCluster_pstn> clustersNxM_pstn);
0352 };
0353 
0354 /*
0355 ████████ ██   ██ ██████     ████████  █████  ██   ██ ███    ███ ██ ███    ██  █████  ████████  ██████  ██████  
0356    ██    ██   ██ ██            ██    ██   ██ ██   ██ ████  ████ ██ ████   ██ ██   ██    ██    ██    ██ ██   ██ 
0357    ██    ███████ █████         ██    ███████ ██   ██ ██ ████ ██ ██ ██ ██  ██ ███████    ██    ██    ██ ██████  
0358    ██    ██   ██ ██            ██    ██   ██ ██   ██ ██  ██  ██ ██ ██  ██ ██ ██   ██    ██    ██    ██ ██   ██ 
0359    ██    ██   ██ ██████        ██    ██   ██ ███████ ██      ██ ██ ██   ████ ██   ██    ██     ██████  ██    ██
0360 */
0361 
0362 std::unique_ptr<NNmodels_GlobalCache> L1NNCaloTauEmulator::initializeGlobalCache(const edm::ParameterSet& iConfig) {
0363   edm::LogInfo("Initialization") << "Init NN models Global Cache " << std::endl;
0364 
0365   std::unique_ptr<NNmodels_GlobalCache> GlobalCache(new NNmodels_GlobalCache);
0366 
0367   GlobalCache->CNNmodel_CB_path = iConfig.getParameter<std::string>("CNNmodel_CB_path");
0368   GlobalCache->DNNident_CB_path = iConfig.getParameter<std::string>("DNNident_CB_path");
0369   GlobalCache->DNNcalib_CB_path = iConfig.getParameter<std::string>("DNNcalib_CB_path");
0370   GlobalCache->CNNmodel_CE_path = iConfig.getParameter<std::string>("CNNmodel_CE_path");
0371   GlobalCache->DNNident_CE_path = iConfig.getParameter<std::string>("DNNident_CE_path");
0372   GlobalCache->DNNcalib_CE_path = iConfig.getParameter<std::string>("DNNcalib_CE_path");
0373   GlobalCache->FeatScaler_CE_path = iConfig.getParameter<std::string>("FeatScaler_CE_path");
0374 
0375   // Create sessions for Tensorflow inferece
0376   (GlobalCache->CNNmodel_CB) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->CNNmodel_CB_path)).fullPath());
0377   (GlobalCache->CNNmodel_CBsession) = tensorflow::createSession((GlobalCache->CNNmodel_CB));
0378 
0379   (GlobalCache->DNNident_CB) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNident_CB_path)).fullPath());
0380   (GlobalCache->DNNident_CBsession) = tensorflow::createSession((GlobalCache->DNNident_CB));
0381 
0382   (GlobalCache->DNNcalib_CB) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNcalib_CB_path)).fullPath());
0383   (GlobalCache->DNNcalib_CBsession) = tensorflow::createSession((GlobalCache->DNNcalib_CB));
0384 
0385   (GlobalCache->CNNmodel_CE) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->CNNmodel_CE_path)).fullPath());
0386   (GlobalCache->CNNmodel_CEsession) = tensorflow::createSession((GlobalCache->CNNmodel_CE));
0387 
0388   (GlobalCache->DNNident_CE) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNident_CE_path)).fullPath());
0389   (GlobalCache->DNNident_CEsession) = tensorflow::createSession((GlobalCache->DNNident_CE));
0390 
0391   (GlobalCache->DNNcalib_CE) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNcalib_CE_path)).fullPath());
0392   (GlobalCache->DNNcalib_CEsession) = tensorflow::createSession((GlobalCache->DNNcalib_CE));
0393 
0394   // Read features scaler
0395   boost::property_tree::read_json(edm::FileInPath((GlobalCache->FeatScaler_CE_path)).fullPath(),
0396                                   (GlobalCache->FeatScaler_CE));
0397 
0398   return GlobalCache;
0399 }
0400 
0401 // ----Constructor and Destructor -----
0402 L1NNCaloTauEmulator::L1NNCaloTauEmulator(const edm::ParameterSet& iConfig, const NNmodels_GlobalCache* globalCache)
0403     : l1TowersToken(consumes<l1tp2::CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("l1CaloTowers"))),
0404       hgcalTowersToken(consumes<l1t::HGCalTowerBxCollection>(iConfig.getParameter<edm::InputTag>("hgcalTowers"))),
0405 
0406       HGClusterToken(
0407           consumes<l1t::HGCalMulticlusterBxCollection>(iConfig.getParameter<edm::InputTag>("HgcalClusters"))),
0408       scenario(UseEmInterp::No),
0409       preEmId(iConfig.getParameter<std::string>("preEmId")),
0410       VsPuId(iConfig.getParameter<edm::ParameterSet>("VsPuId")),
0411 
0412       EcalEtMinForClustering(iConfig.getParameter<double>("EcalEtMinForClustering")),
0413       HcalEtMinForClustering(iConfig.getParameter<double>("HcalEtMinForClustering")),
0414       EtMinForSeeding(iConfig.getParameter<double>("EtMinForSeeding")),
0415       EtaRestriction(iConfig.getParameter<double>("EtaRestriction")),
0416       CB_CE_split(iConfig.getParameter<double>("CB_CE_split")),
0417       PuidThr(iConfig.getParameter<double>("PuidThr")),
0418 
0419       IdWp90_CB(iConfig.getParameter<double>("IdWp90_CB")),
0420       IdWp95_CB(iConfig.getParameter<double>("IdWp95_CB")),
0421       IdWp99_CB(iConfig.getParameter<double>("IdWp99_CB")),
0422 
0423       IdWp90_CE(iConfig.getParameter<double>("IdWp90_CE")),
0424       IdWp95_CE(iConfig.getParameter<double>("IdWp95_CE")),
0425       IdWp99_CE(iConfig.getParameter<double>("IdWp99_CE")),
0426 
0427       DEBUG(iConfig.getParameter<bool>("DEBUG")) {
0428   // Initialize HGCAL BDTs
0429   if (!VsPuId.method().empty()) {
0430     VsPuId.prepareTMVA();
0431   }
0432 
0433   intPuidThr = apintQuantizer(PuidThr, PUID_LSB, PUID_W);
0434   intEtaRestriction = apintQuantizer(EtaRestriction, IETAPHI_LSB, IETAPHI_W);
0435   intCB_CE_split = apintQuantizer(CB_CE_split, IETAPHI_LSB, IETAPHI_W) + 1;
0436 
0437   // Create produced outputs
0438   produces<BXVector<l1t::Tau>>("L1NNCaloTauCollectionBXV");
0439 
0440   // Settings output
0441   edm::LogInfo("Settings") << "EtaRestriction = " << EtaRestriction << " (" << intEtaRestriction << ")"
0442                            << " , CB_CE_split = " << CB_CE_split << "(" << intCB_CE_split
0443                            << ") , EtMinForSeeding = " << EtMinForSeeding
0444                            << " , HcalTpEtMin = " << HcalEtMinForClustering
0445                            << " , EcalTpEtMin = " << EcalEtMinForClustering << " , PuidThr = " << PuidThr << "("
0446                            << intPuidThr << ")" << std::endl;
0447 }
0448 
0449 void L1NNCaloTauEmulator::produce(edm::Event& iEvent, const edm::EventSetup& eSetup) {
0450   // Output collection
0451   std::unique_ptr<BXVector<l1t::Tau>> L1NNCaloTauCollectionBXV(new l1t::TauBxCollection);
0452 
0453   // Create and Fill collection of all calotowers and their attributes
0454   std::vector<SimpleTowerHit> l1CaloTowers;
0455 
0456   iEvent.getByToken(l1TowersToken, l1CaloTowerHandle);
0457   int warnings = 0;
0458   for (auto& hit : *l1CaloTowerHandle.product()) {
0459     // Skip this weird towers and store warning
0460     if (hit.towerIEta() == -1016 && hit.towerIPhi() == -962) {
0461       warnings += 1;
0462       continue;
0463     }
0464 
0465     SimpleTowerHit l1Hit;
0466     l1Hit.isBarrel = 0x1;
0467     l1Hit.l1egTowerEt = apfixedQuantizer(hit.l1egTowerEt(), PTET_LSB, ET_W);
0468     l1Hit.towerEm = apfixedQuantizer(hit.ecalTowerEt(), PTET_LSB, ET_W);
0469     l1Hit.towerHad = apfixedQuantizer(hit.hcalTowerEt(), PTET_LSB, ET_W);
0470     l1Hit.towerEt = apfixedQuantizer(hit.ecalTowerEt() + hit.hcalTowerEt() + hit.l1egTowerEt(), PTET_LSB, ET_W);
0471     l1Hit.towerIeta = hit.towerIEta();
0472     l1Hit.towerIphi = hit.towerIPhi();
0473 
0474     l1CaloTowers.push_back(l1Hit);
0475   }
0476   if (warnings != 0 && DEBUG) {
0477     edm::LogWarning("BrokenTowers") << " ** WARNING : FOUND " << warnings
0478                                     << " TOWERS WITH towerIeta=-1016 AND towerIphi=-962" << std::endl;
0479   }
0480 
0481   iEvent.getByToken(hgcalTowersToken, hgcalTowersHandle);
0482   for (auto& hit : *hgcalTowersHandle.product()) {
0483     SimpleTowerHit l1Hit;
0484     l1Hit.isBarrel = 0x0;
0485     l1Hit.l1egTowerEt = 0.0;
0486     l1Hit.towerEm = apfixedQuantizer(hit.etEm(), PTET_LSB, ET_W);
0487     l1Hit.towerHad = apfixedQuantizer(hit.etHad(), PTET_LSB, ET_W);
0488     l1Hit.towerEt = apfixedQuantizer(hit.etEm() + hit.etHad(), PTET_LSB, ET_W);
0489     l1Hit.towerIeta = makeEndcapHwIEta(hit.eta());
0490     l1Hit.towerIphi = makeEndcapHwIPhi(hit.phi());
0491 
0492     l1CaloTowers.push_back(l1Hit);
0493   }
0494 
0495   // Sort the ECAL+HCAL+L1EGs tower sums based on total ET
0496   std::sort(begin(l1CaloTowers), end(l1CaloTowers), [](const SimpleTowerHit& a, SimpleTowerHit& b) {
0497     return a.towerEt > b.towerEt;
0498   });
0499 
0500   // Create and Fill the collection of 3D clusters and their attributes
0501   std::vector<SimpleHGCluster> AllHGClusters;
0502   iEvent.getByToken(HGClusterToken, HGClusterHandle);
0503 
0504   for (auto cl3dIt = HGClusterHandle->begin(0); cl3dIt != HGClusterHandle->end(0); ++cl3dIt) {
0505     auto& cl3d = *cl3dIt;
0506 
0507     // Implement cl3d PU ID as done in
0508     // https://github.com/cms-sw/cmssw/blob/master/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc#L120
0509     bool isEM = preEmId(*cl3dIt);
0510     l1t::PFCluster cluster(cl3d.pt(), cl3d.eta(), cl3d.phi(), cl3d.hOverE());
0511     if (scenario == UseEmInterp::EmOnly)  // for emID objs, use EM interp as pT and set H = 0
0512     {
0513       if (isEM) {
0514         float pt_new = cl3d.iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
0515         float hoe_new = 0.;
0516         cluster = l1t::PFCluster(pt_new, cl3d.eta(), cl3d.phi(), hoe_new, isEM);
0517       }
0518     } else if (scenario == UseEmInterp::AllKeepHad)  // for all objs, replace EM part with EM interp, preserve H
0519     {
0520       float had_old = cl3d.pt() - cluster.emEt();
0521       float em_new = cl3d.iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
0522       float pt_new = had_old + em_new;
0523       float hoe_new = em_new > 0 ? (had_old / em_new) : -1;
0524       cluster = l1t::PFCluster(pt_new, cl3d.eta(), cl3d.phi(), hoe_new, isEM);
0525     } else if (scenario == UseEmInterp::AllKeepTot)  // for all objs, replace EM part with EM interp, preserve pT
0526     {
0527       float em_new = cl3d.iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
0528       float hoe_new = em_new > 0 ? (cl3d.pt() / em_new - 1) : -1;
0529       cluster = l1t::PFCluster(cl3d.pt(), cl3d.eta(), cl3d.phi(), hoe_new, isEM);
0530     }
0531 
0532     float idScore = -1.;
0533     if (!VsPuId.method().empty()) {
0534       idScore = VsPuId.passID(*cl3dIt, cluster);
0535       idScore = cluster.egVsPUMVAOut();
0536     }
0537 
0538     float eta_hgcalCoord = correctInputEtaCl3d(cl3d.eta());
0539     float meanz_hgcalCoord = correctInputMeanzCl3d(cl3d.zBarycenter());
0540 
0541     SimpleHGCluster HGCluster;
0542     HGCluster.pt = apfixedQuantizer(cl3d.pt(), PTET_LSB, PT_W);
0543     HGCluster.eta = apintQuantizer(eta_hgcalCoord, ETAPHI_LSB, ETAPHI_W);
0544     HGCluster.phi = apintQuantizer(cl3d.phi(), ETAPHI_LSB, ETAPHI_W);
0545     HGCluster.showerlength = cl3d.showerLength();
0546     HGCluster.coreshowerlength = cl3d.coreShowerLength();
0547     HGCluster.spptot = apintQuantizer(cl3d.sigmaPhiPhiTot(), SHAPEFEAT_LSB, SHAPEFEAT_W);
0548     HGCluster.szz = apintQuantizer(cl3d.sigmaZZ(), SZZ_LSB, SHAPEFEAT_W);
0549     HGCluster.srrtot = apintQuantizer(cl3d.sigmaRRTot(), SHAPEFEAT_LSB, SHAPEFEAT_W);
0550     HGCluster.meanz = apintQuantizer(meanz_hgcalCoord, MEANZ_LSB, MEANZ_W);
0551     HGCluster.PUid = apintQuantizer(idScore, PUID_LSB, PUID_W);
0552 
0553     AllHGClusters.push_back(HGCluster);
0554   }
0555 
0556   // Order the collection in pt (the input to the GCT will be pt ordered)
0557   std::sort(begin(AllHGClusters), end(AllHGClusters), [](const SimpleHGCluster& a, SimpleHGCluster& b) {
0558     return a.pt > b.pt;
0559   });
0560 
0561   /*
0562   // END OF SOFTWARE PRECISION SECTION
0563   // up to here treated inputs from simulation with SW precision
0564   // to massage them into the HW precision varibales as they are
0565   // forseen (roughly) to be available at the GCT Sum card level
0566   // ------------------------------------------------------------- */
0567 
0568   // Make NxM TowerClusters and HGClusters collections for TauMinator
0569   std::vector<InputTowerCluster> l1TowerClustersNxM_CB;
0570   std::vector<InputTowerCluster_pstn> l1TowerClustersNxM_CB_pstn;
0571   std::vector<InputTowerCluster> l1TowerClustersNxM_CE;
0572   std::vector<InputTowerCluster_pstn> l1TowerClustersNxM_CE_pstn;
0573   std::vector<InputHGCluster> HGClusters;
0574 
0575   // Supporting collection of endcap clusters before cl3d matching
0576   std::vector<InputTowerCluster> AllL1TowerClustersNxM_CE;
0577   std::vector<InputTowerCluster_pstn> AllL1TowerClustersNxM_CE_pstn;
0578 
0579   int Nclusters_CB = 0;
0580   int AllNclusters_CE = 0;
0581   bool caloTauSeedingFinished = false;
0582   // Loop for seeding of clNxM objects
0583   while (!caloTauSeedingFinished) {
0584     InputTowerCluster clNxM;
0585     clNxM.init();
0586     InputTowerCluster_pstn clNxM_pstn;
0587     bool seeded = false;
0588 
0589     for (auto& l1CaloTower : l1CaloTowers) {
0590       // Skip seeding in towers that would make the cluster extend in HF
0591       // Skip l1CaloTowers which are already used by this clusters' mask
0592       if (ap_abs(l1CaloTower.towerIeta) > Eta_limit || ap_abs(l1CaloTower.towerIeta) > intEtaRestriction ||
0593           l1CaloTower.stale4seed) {
0594         continue;
0595       }
0596 
0597       // If not seded do the seeding
0598       if (!seeded) {
0599         // The leading unused tower has ET < min, stop jet clustering
0600         if (l1CaloTower.towerEt < EtMinForSeeding) {
0601           caloTauSeedingFinished = true;
0602           continue;
0603         }
0604 
0605         clNxM.fill(seedIdx, l1CaloTower);
0606         clNxM_pstn.fill(l1CaloTower);
0607         if (l1CaloTower.isBarrel) {
0608           clNxM.barrelSeeded = 0x1;
0609         }
0610 
0611         l1CaloTower.stale4seed = 0x1;
0612         l1CaloTower.stale = 0x1;
0613         seeded = true;
0614 
0615         continue;
0616       }
0617 
0618       dIEtaPhi_t d_iEta = tower_dIEta(l1CaloTower.towerIeta, clNxM_pstn.seedIeta);
0619       dIEtaPhi_t d_iPhi = dPhi<dIEtaPhi_t, IPhi_t>(l1CaloTower.towerIphi, clNxM_pstn.seedIphi);
0620 
0621       // Stale tower for seeding if it would lead to overalp between clusters
0622       if ((ap_abs(d_iEta) <= IEta_dim - 1 && ap_abs(d_iPhi) <= IPhi_dim - 1)) {
0623         l1CaloTower.stale4seed = 0x1;
0624       }
0625 
0626     }  // End for loop over TPs
0627 
0628     // Pushback seeds split in barrel and endcap
0629     if (seeded) {
0630       if (ap_abs(clNxM_pstn.seedIeta) <= intCB_CE_split) {
0631         l1TowerClustersNxM_CB.push_back(clNxM);
0632         l1TowerClustersNxM_CB_pstn.push_back(clNxM_pstn);
0633         Nclusters_CB++;
0634       } else {
0635         AllL1TowerClustersNxM_CE.push_back(clNxM);
0636         AllL1TowerClustersNxM_CE_pstn.push_back(clNxM_pstn);
0637         AllNclusters_CE++;
0638       }
0639     }
0640 
0641   }  // End while loop of TowerClusters seeding
0642 
0643   // Loop for barrel NxM TowerClusters clustering starting from the seeds
0644   for (int clNxMIdx = 0; clNxMIdx < Nclusters_CB; clNxMIdx++) {
0645     for (auto& l1CaloTower : l1CaloTowers) {
0646       // Skip l1CaloTowers which are already used
0647       if (l1CaloTower.stale) {
0648         continue;
0649       }
0650 
0651       dIEtaPhi_t d_iEta = tower_dIEta(l1CaloTower.towerIeta, l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIeta);
0652       dIEtaPhi_t d_iPhi =
0653           dPhi<dIEtaPhi_t, IPhi_t>(l1CaloTower.towerIphi, l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIphi);
0654       int hitIdx = d_iEta * 9 + d_iPhi + seedIdx;
0655 
0656       // Cluster all towers in a NxM towers mask
0657       if ((ap_abs(d_iEta) <= (IEta_dim - 1) / 2 && ap_abs(d_iPhi) <= (IPhi_dim - 1) / 2)) {
0658         l1TowerClustersNxM_CB[clNxMIdx].fill(hitIdx, l1CaloTower);
0659         l1CaloTower.stale = 0x1;
0660       }
0661 
0662     }  // End for loop over TPs
0663 
0664   }  // End while loop of barrel TowerClusters creation
0665 
0666   // In the endcap cross-loop over clNxM and cl3d to match them
0667   // (we can do it before full clustering just using the seed info)
0668   int Nclusters_CE = 0;
0669   for (int clNxMIdx = 0; clNxMIdx < AllNclusters_CE; clNxMIdx++) {
0670     bool matched = false;
0671     for (auto& HGCluster : AllHGClusters) {
0672       // In case the clNxM or HGCluster have already been matched just continue through the list to the end
0673       // only use cl3ds above 4GeV and above -0.10 pu id
0674       if (matched || HGCluster.stale || HGCluster.pt < Pt_t(4.) || HGCluster.PUid < intPuidThr) {
0675         continue;
0676       }
0677 
0678       dEtaPhi_t d_iEta = tw2cl_dEta(HGCluster.eta, AllL1TowerClustersNxM_CE_pstn[clNxMIdx].seedIeta);
0679       dEtaPhi_t d_iPhi = tw2cl_dPhi(HGCluster.phi, AllL1TowerClustersNxM_CE_pstn[clNxMIdx].seedIphi);
0680       matched = d_iEta * d_iEta + d_iPhi * d_iPhi < R2cone;
0681 
0682       if (matched) {
0683         HGCluster.stale = 0x1;
0684         InputHGCluster cl3d;
0685         cl3d.fill(HGCluster);
0686         HGClusters.push_back(cl3d);
0687         l1TowerClustersNxM_CE.push_back(AllL1TowerClustersNxM_CE[clNxMIdx]);
0688         l1TowerClustersNxM_CE_pstn.push_back(AllL1TowerClustersNxM_CE_pstn[clNxMIdx]);
0689         Nclusters_CE++;
0690       }
0691 
0692     }  // End for loop over cl3ds
0693 
0694   }  // End for loop over clNxM
0695 
0696   // Loop for endcap matched NxM TowerClusters clustering starting from the seeds just found
0697   for (int clNxMIdx = 0; clNxMIdx < Nclusters_CE; clNxMIdx++) {
0698     for (auto& l1CaloTower : l1CaloTowers) {
0699       // Skip l1CaloTowers which are already used
0700       if (l1CaloTower.stale) {
0701         continue;
0702       }
0703 
0704       dIEtaPhi_t d_iEta = tower_dIEta(l1CaloTower.towerIeta, l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIeta);
0705       dIEtaPhi_t d_iPhi =
0706           dPhi<dIEtaPhi_t, IPhi_t>(l1CaloTower.towerIphi, l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIphi);
0707       int hitIdx = d_iEta * 9 + d_iPhi + seedIdx;
0708 
0709       // Cluster all towers in a NxM towers mask
0710       if ((ap_abs(d_iEta) <= (IEta_dim - 1) / 2 && ap_abs(d_iPhi) <= (IPhi_dim - 1) / 2)) {
0711         l1TowerClustersNxM_CE[clNxMIdx].fill(hitIdx, l1CaloTower);
0712         l1CaloTower.stale = 0x1;
0713       }
0714 
0715     }  // End for loop over TPs
0716 
0717   }  // End while loop of barrel TowerClusters creation
0718 
0719   // Barrel TauMinator application
0720   tensorflow::setLogging("2");
0721   int batchSize_CB = (int)(Nclusters_CB);
0722   tensorflow::TensorShape imageShape_CB({batchSize_CB, IEta_dim, IPhi_dim, 2});
0723   tensorflow::TensorShape positionShape_CB({batchSize_CB, 2});
0724   tensorflow::Tensor TowerClusterImage_CB(tensorflow::DT_FLOAT, imageShape_CB);
0725   tensorflow::Tensor TowerClusterPosition_CB(tensorflow::DT_FLOAT, positionShape_CB);
0726 
0727   for (int clNxMIdx = 0; clNxMIdx < Nclusters_CB; clNxMIdx++) {
0728     // Fill inputs for Tensorflow inference
0729     for (int eta = 0; eta < IEta_dim; ++eta) {
0730       for (int phi = 0; phi < IPhi_dim; ++phi) {
0731         int towerIdx = eta * IPhi_dim + phi;
0732         TowerClusterImage_CB.tensor<float, 4>()(clNxMIdx, eta, phi, 0) =
0733             (l1TowerClustersNxM_CB[clNxMIdx].towerHits[towerIdx].l1egTowerEt.to_float() +
0734              l1TowerClustersNxM_CB[clNxMIdx].towerHits[towerIdx].towerEm.to_float());
0735         TowerClusterImage_CB.tensor<float, 4>()(clNxMIdx, eta, phi, 1) =
0736             (l1TowerClustersNxM_CB[clNxMIdx].towerHits[towerIdx].towerHad.to_float());
0737       }
0738     }
0739 
0740     TowerClusterPosition_CB.tensor<float, 2>()(clNxMIdx, 0) = floatIEta(l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIeta);
0741     TowerClusterPosition_CB.tensor<float, 2>()(clNxMIdx, 1) = floatIPhi(l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIphi);
0742   }
0743 
0744   if (batchSize_CB >
0745       0)  // from CMSSW_14_0_X tensorflow does not seem to be able to deal with a tensor of dimension 0 anymore
0746   {
0747     // Apply CNN model
0748     tensorflow::NamedTensorList CNNmodel_CBinputList = {{"TowerClusterImage", TowerClusterImage_CB},
0749                                                         {"TowerClusterPosition", TowerClusterPosition_CB}};
0750     std::vector<tensorflow::Tensor> CNNmodel_CBoutputs;
0751     tensorflow::run((globalCache()->CNNmodel_CBsession),
0752                     CNNmodel_CBinputList,
0753                     {"TauMinator_CB_conv/middleMan/concat"},
0754                     &CNNmodel_CBoutputs);
0755     tensorflow::NamedTensorList DNN_CBinputsList = {{"middleMan", CNNmodel_CBoutputs[0]}};
0756 
0757     // Apply DNN for identification
0758     std::vector<tensorflow::Tensor> DNN_CBoutputsIdent;
0759     tensorflow::run((globalCache()->DNNident_CBsession),
0760                     DNN_CBinputsList,
0761                     {"TauMinator_CB_ident/sigmoid_IDout/Sigmoid"},
0762                     &DNN_CBoutputsIdent);
0763 
0764     // Apply DNN for calibration
0765     std::vector<tensorflow::Tensor> DNN_CBoutputsCalib;
0766     tensorflow::run((globalCache()->DNNcalib_CBsession),
0767                     DNN_CBinputsList,
0768                     {"TauMinator_CB_calib/DNNout/MatMul"},
0769                     &DNN_CBoutputsCalib);
0770 
0771     // Fill the output collection of L1 taus with the barrel candidates
0772     for (int clNxMIdx = 0; clNxMIdx < Nclusters_CB; clNxMIdx++) {
0773       l1t::Tau l1Tau =
0774           MakeTauCandidate(true, clNxMIdx, DNN_CBoutputsIdent, DNN_CBoutputsCalib, l1TowerClustersNxM_CB_pstn);
0775       if (l1Tau.pt() < 0) {
0776         continue;
0777       }
0778       L1NNCaloTauCollectionBXV->push_back(0, l1Tau);
0779     }
0780   }
0781 
0782   // Endcap TauMinator application
0783   int batchSize_CE = (int)(Nclusters_CE);
0784   tensorflow::TensorShape imageShape_CE({batchSize_CE, IEta_dim, IPhi_dim, 2});
0785   tensorflow::TensorShape positionShape_CE({batchSize_CE, 2});
0786   tensorflow::TensorShape cl3dfeatShape_CE({batchSize_CE, 8});
0787   tensorflow::Tensor TowerClusterImage_CE(tensorflow::DT_FLOAT, imageShape_CE);
0788   tensorflow::Tensor TowerClusterPosition_CE(tensorflow::DT_FLOAT, positionShape_CE);
0789   tensorflow::Tensor Cl3dShapeFeatures_CE(tensorflow::DT_FLOAT, cl3dfeatShape_CE);
0790 
0791   for (int clNxMIdx = 0; clNxMIdx < Nclusters_CE; clNxMIdx++) {
0792     // Indexing of cl3ds is the same as the one of clNxMs
0793     InputHGCluster HGClu = HGClusters[clNxMIdx];
0794 
0795     // Fill inputs for Tensorflow inference
0796     for (int eta = 0; eta < IEta_dim; ++eta) {
0797       for (int phi = 0; phi < IPhi_dim; ++phi) {
0798         int towerIdx = eta * IPhi_dim + phi;
0799         TowerClusterImage_CE.tensor<float, 4>()(clNxMIdx, eta, phi, 0) =
0800             (l1TowerClustersNxM_CE[clNxMIdx].towerHits[towerIdx].l1egTowerEt.to_float() +
0801              l1TowerClustersNxM_CE[clNxMIdx].towerHits[towerIdx].towerEm.to_float());
0802         TowerClusterImage_CE.tensor<float, 4>()(clNxMIdx, eta, phi, 1) =
0803             (l1TowerClustersNxM_CE[clNxMIdx].towerHits[towerIdx].towerHad.to_float());
0804       }
0805     }
0806 
0807     TowerClusterPosition_CE.tensor<float, 2>()(clNxMIdx, 0) = floatIEta(l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIeta);
0808     TowerClusterPosition_CE.tensor<float, 2>()(clNxMIdx, 1) = floatIPhi(l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIphi);
0809 
0810     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 0) = inputScaler(HGClu.pt.to_float(), "pt");
0811     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 1) = inputScaler(abs(floatEta(HGClu.eta)), "eta");
0812     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 2) = inputScaler(HGClu.showerlength.to_float(), "showerlength");
0813     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 3) =
0814         inputScaler(HGClu.coreshowerlength.to_float(), "coreshowerlength");
0815     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 4) = inputScaler(floatShape(HGClu.spptot), "spptot");
0816     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 5) = inputScaler(floatSzz(HGClu.szz), "szz");
0817     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 6) = inputScaler(floatShape(HGClu.srrtot), "srrtot");
0818     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 7) = inputScaler(floatMeanZHgcalCoord(HGClu.meanz), "meanz");
0819   }
0820 
0821   if (batchSize_CE >
0822       0)  // from CMSSW_14_0_X tensorflow does not seem to be able to deal with a tensor of dimension 0 anymore
0823   {
0824     // Apply CNN model
0825     tensorflow::NamedTensorList CNNmodel_CEinputList = {{"TowerClusterImage", TowerClusterImage_CE},
0826                                                         {"TowerClusterPosition", TowerClusterPosition_CE},
0827                                                         {"AssociatedCl3dFeatures", Cl3dShapeFeatures_CE}};
0828     std::vector<tensorflow::Tensor> CNNmodel_CEoutputs;
0829     tensorflow::run((globalCache()->CNNmodel_CEsession),
0830                     CNNmodel_CEinputList,
0831                     {"TauMinator_CE_conv/middleMan/concat"},
0832                     &CNNmodel_CEoutputs);
0833     tensorflow::NamedTensorList DNN_CEinputsList = {{"middleMan", CNNmodel_CEoutputs[0]}};
0834 
0835     // Apply DNN for identification
0836     std::vector<tensorflow::Tensor> DNN_CEoutputsIdent;
0837     tensorflow::run((globalCache()->DNNident_CEsession),
0838                     DNN_CEinputsList,
0839                     {"TauMinator_CE_ident/sigmoid_IDout/Sigmoid"},
0840                     &DNN_CEoutputsIdent);
0841 
0842     // Apply DNN for calibration
0843     std::vector<tensorflow::Tensor> DNN_CEoutputsCalib;
0844     tensorflow::run((globalCache()->DNNcalib_CEsession),
0845                     DNN_CEinputsList,
0846                     {"TauMinator_CE_calib/LIN_DNNout/Relu"},
0847                     &DNN_CEoutputsCalib);
0848 
0849     // Fill the output collection of L1 taus with the endcap candidates
0850     for (int clNxMIdx = 0; clNxMIdx < Nclusters_CE; clNxMIdx++) {
0851       l1t::Tau l1Tau =
0852           MakeTauCandidate(false, clNxMIdx, DNN_CEoutputsIdent, DNN_CEoutputsCalib, l1TowerClustersNxM_CE_pstn);
0853       if (l1Tau.pt() < 0) {
0854         continue;
0855       }
0856       L1NNCaloTauCollectionBXV->push_back(0, l1Tau);
0857     }
0858   }
0859 
0860   // Fill output
0861   iEvent.put(std::move(L1NNCaloTauCollectionBXV), "L1NNCaloTauCollectionBXV");
0862 
0863 }  // End of produce function
0864 
0865 template <class outPrecision, class inPrecision>
0866 outPrecision L1NNCaloTauEmulator::dPhi(inPrecision iPhi_1, inPrecision iPhi_2) {
0867   outPrecision dphi = iPhi_1 - iPhi_2;
0868 
0869   outPrecision dphi0 = dphi > outPrecision(INTPHI_PI) ? outPrecision(dphi - INTPHI_2PI) : dphi;
0870   outPrecision dphi1 = dphi <= outPrecision(-INTPHI_PI) ? outPrecision(dphi + INTPHI_2PI) : dphi;
0871 
0872   outPrecision result = dphi > outPrecision(0) ? dphi0 : dphi1;
0873 
0874   return result;
0875 }
0876 
0877 L1NNCaloTauEmulator::dIEtaPhi_t L1NNCaloTauEmulator::tower_dIEta(IEta_t iEta_1, IEta_t iEta_2) {
0878   ap_int<12> mult = iEta_1 * iEta_2;
0879   dIEtaPhi_t result = iEta_1 - iEta_2;
0880   if (mult < 0) {
0881     result = iEta_1 > 0 ? result - 1 : result + 1;
0882   }
0883 
0884   return result;
0885 }
0886 
0887 L1NNCaloTauEmulator::dEtaPhi_t L1NNCaloTauEmulator::tw2cl_dPhi(EtaPhi_t iPhi_1, IPhi_t iPhi_2) {
0888   EtaPhi_t shiftediPhi_2 = iPhi_2 <= IPhi_t(36) ? EtaPhi_t(iPhi_2) : EtaPhi_t(iPhi_2 - INTPHI_2PI + 1);
0889 
0890   EtaPhi_t fineiPhi_2 = shiftediPhi_2 * (IETAPHI_LSB / ETAPHI_LSB);
0891   // subrtaction of half rescaling corrects from edge to center of tower
0892   fineiPhi_2 = fineiPhi_2 > EtaPhi_t(0) ? EtaPhi_t(fineiPhi_2 - (IETAPHI_LSB / ETAPHI_LSB) / 2)
0893                                         : EtaPhi_t(fineiPhi_2 + (IETAPHI_LSB / ETAPHI_LSB) / 2);
0894 
0895   return dPhi<dEtaPhi_t, EtaPhi_t>(iPhi_1, fineiPhi_2);
0896 }
0897 
0898 L1NNCaloTauEmulator::dEtaPhi_t L1NNCaloTauEmulator::tw2cl_dEta(EtaPhi_t iEta_1, IEta_t iEta_2) {
0899   // change from hgcal frame to barrel-centered frame
0900   EtaPhi_t framechangeCl3d = 303;  // 303*pi/720 = 1.322
0901   iEta_1 = iEta_1 > EtaPhi_t(0) ? EtaPhi_t(iEta_1 + framechangeCl3d) : EtaPhi_t(iEta_1 - framechangeCl3d);
0902 
0903   // the actual depth is 330 but 329 corrects for 0.0808 tower
0904   EtaPhi_t barrelEtaDepth = 329;
0905   EtaPhi_t fineiEta_2 = barrelEtaDepth + (iEta_2 - IETAHGCAL_OFFSET) * (IETAHGCAL_LSB / ETAPHI_LSB);
0906 
0907   return iEta_1 - fineiEta_2;
0908 }
0909 
0910 L1NNCaloTauEmulator::IEta_t L1NNCaloTauEmulator::makeEndcapHwIEta(float eta) {
0911   IEta_t ieta = floor(eta / IETAHGCAL_LSB);
0912   // +1 because flooring gets it 1 unit lower when negative
0913   ieta = ieta < IEta_t(0) ? IEta_t(ieta + 1) : ieta;
0914 
0915   return ieta;
0916 }
0917 
0918 L1NNCaloTauEmulator::IPhi_t L1NNCaloTauEmulator::makeEndcapHwIPhi(float phi) {
0919   phi = phi < 0 ? phi + 2 * M_PI : phi;
0920 
0921   // +1 because tower 0 does not exist
0922   return floor(phi / IETAPHI_LSB) + 1;
0923 }
0924 
0925 template <int W>
0926 ap_int<W> L1NNCaloTauEmulator::ap_abs(ap_int<W> x) {
0927   ap_int<W> result;
0928   if (x < 0) {
0929     result = -x;
0930   } else {
0931     result = x;
0932   }
0933 
0934   return result;
0935 }
0936 
0937 template <int W, int I, ap_q_mode _AP_Q, ap_o_mode _AP_O>
0938 ap_ufixed<W, I> L1NNCaloTauEmulator::ap_abs(ap_fixed<W, I, _AP_Q, _AP_O> x) {
0939   ap_ufixed<W, I> result;
0940   if (x < 0) {
0941     result = -x;
0942   } else {
0943     result = x;
0944   }
0945 
0946   return result;
0947 }
0948 
0949 float L1NNCaloTauEmulator::apfixedQuantizer(float inputF, float LSB, int nbits) {
0950   return min(floor(inputF / LSB), float(pow(2, nbits) - 1)) * LSB;
0951 }
0952 
0953 int L1NNCaloTauEmulator::apintQuantizer(float inputF, float LSB, int nbits) {
0954   return min(floor(inputF / LSB), float(pow(2, nbits) - 1));
0955 }
0956 
0957 float L1NNCaloTauEmulator::inputScaler(float inputF, std::string feature) {
0958   float mean = (globalCache()->FeatScaler_CE).get_child(feature).get<float>("mean");
0959   float std = (globalCache()->FeatScaler_CE).get_child(feature).get<float>("std");
0960 
0961   return (inputF - mean) / std;
0962 }
0963 
0964 float L1NNCaloTauEmulator::correctInputEtaCl3d(float eta) {
0965   return eta > 0 ? eta - ETAHGCAL_OFFSET : eta + ETAHGCAL_OFFSET;
0966 }
0967 
0968 float L1NNCaloTauEmulator::correctInputMeanzCl3d(float meanz) { return CM2MM * (abs(meanz) - MEANZ_OFFSET); }
0969 
0970 float L1NNCaloTauEmulator::floatIEta(IEta_t eta) {
0971   // transform eta of towers from integer to float, correcting for different barrel/endcap LSB
0972   float feta;
0973   if (abs(eta) > IETAHGCAL_OFFSET) {
0974     if (eta > 0) {
0975       feta = IETAHGCAL_OFFSET * IETAPHI_LSB - (IETAHGCAL_LSB - IETAHGCAL_LSBp) +
0976              (eta.to_float() - IETAHGCAL_OFFSET) * IETAHGCAL_LSB;
0977     } else {
0978       feta = -IETAHGCAL_OFFSET * IETAPHI_LSB + (IETAHGCAL_LSB - IETAHGCAL_LSBp) +
0979              (eta.to_float() + IETAHGCAL_OFFSET) * IETAHGCAL_LSB;
0980     }
0981   } else {
0982     feta = eta.to_float() * IETAPHI_LSB;
0983   }
0984 
0985   // shift by half a tower to consider the tower center instead of the edge
0986   return feta > 0 ? feta - IETAPHI_LSB / 2 : feta + IETAPHI_LSB / 2;
0987 }
0988 
0989 float L1NNCaloTauEmulator::floatIPhi(IPhi_t phi) {
0990   float fphi = phi.to_float();
0991   // add 2pi + 1 because tower 0 does not exist
0992   fphi = fphi > INTPHI_PI ? fphi - INTPHI_2PI + 1 : fphi;
0993   fphi *= IETAPHI_LSB;
0994 
0995   // shift by half a tower to consider the tower center instead of the edge
0996   return fphi > 0 ? fphi - IETAPHI_LSB / 2 : fphi + IETAPHI_LSB / 2;
0997 }
0998 
0999 l1t::Tau L1NNCaloTauEmulator::MakeTauCandidate(
1000     bool isBarrel,
1001     int clNxMIdx,
1002     std::vector<tensorflow::Tensor> outputsIdent,
1003     std::vector<tensorflow::Tensor> outputsCalib,
1004     std::vector<L1NNCaloTauEmulator::InputTowerCluster_pstn> clustersNxM_pstn) {
1005   int seedIeta = clustersNxM_pstn[clNxMIdx].seedIeta;
1006   int seedIphi = clustersNxM_pstn[clNxMIdx].seedIphi;
1007 
1008   if (seedIeta > intEtaRestriction) {
1009     return l1t::Tau(reco::Candidate::PolarLorentzVector(-1, 0, 0, 0), -1, 0, 0, 0, 0);
1010     ;
1011   }
1012 
1013   float tau_IDscore = outputsIdent[0].matrix<float>()(0, clNxMIdx);
1014   float tau_calibPt = outputsCalib[0].matrix<float>()(0, clNxMIdx);
1015   float tau_eta = floatIEta(seedIeta);
1016   float tau_phi = floatIPhi(seedIphi);
1017 
1018   // Assign increasing quality to higher scoring candidates
1019   int quality = 0;
1020   if (isBarrel) {
1021     // 99% WP
1022     if (tau_IDscore > IdWp99_CB) {
1023       quality = 1;
1024     }
1025     // 95% WP
1026     if (tau_IDscore > IdWp95_CB) {
1027       quality = 2;
1028     }
1029     // 90% WP
1030     if (tau_IDscore > IdWp90_CB) {
1031       quality = 3;
1032     }
1033   } else {
1034     // 99% WP
1035     if (tau_IDscore > IdWp99_CE) {
1036       quality = 1;
1037     }
1038     // 95% WP
1039     if (tau_IDscore > IdWp95_CE) {
1040       quality = 2;
1041     }
1042     // 90% WP
1043     if (tau_IDscore > IdWp90_CE) {
1044       quality = 3;
1045     }
1046   }
1047 
1048   reco::Candidate::PolarLorentzVector tauP4 = reco::Candidate::PolarLorentzVector(tau_calibPt, tau_eta, tau_phi, 0);
1049 
1050   // store ID score multiplied by 10E4 to have good precision even using the Phase1 tau int iso format
1051   // (this is stored just in case for possible additional offline studies)
1052   // tau initialisation =  (p4,    pt,          eta,     phi,     qual,    iso)
1053   l1t::Tau l1Tau = l1t::Tau(tauP4, tau_calibPt, tau_eta, tau_phi, quality, tau_IDscore * 10E4);
1054   l1Tau.setTowerIEta(seedIeta);
1055   l1Tau.setTowerIPhi(seedIphi);
1056 
1057   return l1Tau;
1058 }
1059 
1060 void L1NNCaloTauEmulator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1061   edm::ParameterSetDescription desc;
1062 
1063   desc.add<edm::InputTag>("l1CaloTowers", edm::InputTag("l1tEGammaClusterEmuProducer", "L1CaloTowerCollection"));
1064   desc.add<edm::InputTag>("hgcalTowers", edm::InputTag("l1tHGCalTowerProducer", "HGCalTowerProcessor"));
1065   desc.add<edm::InputTag>("HgcalClusters",
1066                           edm::InputTag("l1tHGCalBackEndLayer2Producer", "HGCalBackendLayer2Processor3DClustering"));
1067 
1068   desc.add<std::string>("preEmId", "hOverE < 0.3 && hOverE >= 0");
1069   {
1070     edm::ParameterSetDescription psd0;
1071     psd0.add<bool>("isPUFilter", true);
1072     psd0.add<std::string>("preselection", "");
1073     psd0.add<std::string>("method", "BDT");
1074     {
1075       edm::ParameterSetDescription vpsd2;
1076       vpsd2.add<std::string>("name");
1077       vpsd2.add<std::string>("value");
1078       std::vector<edm::ParameterSet> temp2;
1079       temp2.reserve(5);
1080       {
1081         edm::ParameterSet temp3;
1082         temp3.addParameter<std::string>("name", "eMax");
1083         temp3.addParameter<std::string>("value", "eMax()");
1084         temp2.push_back(temp3);
1085       }
1086       {
1087         edm::ParameterSet temp3;
1088         temp3.addParameter<std::string>("name", "eMaxOverE");
1089         temp3.addParameter<std::string>("value", "eMax()/energy()");
1090         temp2.push_back(temp3);
1091       }
1092       {
1093         edm::ParameterSet temp3;
1094         temp3.addParameter<std::string>("name", "sigmaPhiPhiTot");
1095         temp3.addParameter<std::string>("value", "sigmaPhiPhiTot()");
1096         temp2.push_back(temp3);
1097       }
1098       {
1099         edm::ParameterSet temp3;
1100         temp3.addParameter<std::string>("name", "sigmaRRTot");
1101         temp3.addParameter<std::string>("value", "sigmaRRTot()");
1102         temp2.push_back(temp3);
1103       }
1104       {
1105         edm::ParameterSet temp3;
1106         temp3.addParameter<std::string>("name", "triggerCells90percent");
1107         temp3.addParameter<std::string>("value", "triggerCells90percent()");
1108         temp2.push_back(temp3);
1109       }
1110       psd0.addVPSet("variables", vpsd2, temp2);
1111     }
1112     psd0.add<std::string>(
1113         "weightsFile", "L1Trigger/Phase2L1ParticleFlow/data/hgcal_egID/Photon_Pion_vs_Neutrino_BDTweights_1116.xml.gz");
1114     psd0.add<std::string>("wp", "-0.10");
1115     desc.add<edm::ParameterSetDescription>("VsPuId", psd0);
1116   }
1117 
1118   desc.add<double>("EcalEtMinForClustering", 0.0);
1119   desc.add<double>("HcalEtMinForClustering", 0.0);
1120   desc.add<double>("EtMinForSeeding", 2.5);
1121   desc.add<double>("EtaRestriction", 2.4);
1122   desc.add<double>("CB_CE_split", 1.55);
1123   desc.add<double>("PuidThr", -0.1);
1124 
1125   desc.add<std::string>("CNNmodel_CB_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/CNNmodel_CB.pb");
1126   desc.add<std::string>("DNNident_CB_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNident_CB.pb");
1127   desc.add<std::string>("DNNcalib_CB_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNcalib_CB.pb");
1128   desc.add<std::string>("CNNmodel_CE_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/CNNmodel_CE.pb");
1129   desc.add<std::string>("DNNident_CE_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNident_CE.pb");
1130   desc.add<std::string>("DNNcalib_CE_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNcalib_CE.pb");
1131   desc.add<std::string>("FeatScaler_CE_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/Cl3dFeatScaler_CE.json");
1132 
1133   desc.add<double>("IdWp90_CB", 0.706);
1134   desc.add<double>("IdWp95_CB", 0.3432);
1135   desc.add<double>("IdWp99_CB", 0.0337);
1136   desc.add<double>("IdWp90_CE", 0.5711);
1137   desc.add<double>("IdWp95_CE", 0.2742);
1138   desc.add<double>("IdWp99_CE", 0.0394);
1139 
1140   desc.add<bool>("DEBUG", false);
1141 
1142   descriptions.add("l1tNNCaloTauEmulator", desc);
1143 }
1144 
1145 DEFINE_FWK_MODULE(L1NNCaloTauEmulator);
[ 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!