Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:20:09

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   // Class for the towers info as they should be in GCT
0247   class SimpleTowerHit {
0248   public:
0249     IEta_t towerIeta = 0;
0250     IPhi_t towerIphi = 0;
0251     Et_t towerEm = 0.;
0252     Et_t towerHad = 0.;
0253     Et_t l1egTowerEt = 0.;
0254     Et_t towerEt = 0.;
0255     ap_uint<1> isBarrel = 0x1;
0256     ap_uint<1> stale = 0x0;
0257     ap_uint<1> stale4seed = 0x0;
0258   };
0259 
0260   // Class for the clusters info as they should arrive from HGCAL
0261   class SimpleHGCluster {
0262   public:
0263     Pt_t pt;
0264     EtaPhi_t eta;
0265     EtaPhi_t phi;
0266     ShowLen_t showerlength;
0267     ShowLen_t coreshowerlength;
0268     ShapeFeat_t spptot;
0269     ShapeFeat_t szz;
0270     ShapeFeat_t srrtot;
0271     Meanz_t meanz;
0272     PUid_t PUid;
0273     ap_uint<1> stale = 0x0;
0274   };
0275 
0276   // Classes for the tower clusters
0277   class SimplifiedTower {
0278   public:
0279     Et_t towerEm = 0.;
0280     Et_t towerHad = 0.;
0281     Et_t l1egTowerEt = 0.;
0282 
0283     void fill(SimpleTowerHit Tower) {
0284       towerEm = Tower.towerEm;
0285       towerHad = Tower.towerHad;
0286       l1egTowerEt = Tower.l1egTowerEt;
0287     }
0288   };
0289 
0290   class InputTowerCluster {
0291   public:
0292     SimplifiedTower towerHits[45];
0293     ap_uint<1> barrelSeeded = 0x0;
0294     ap_uint<1> filled[45];
0295 
0296     void fill(int idx, SimpleTowerHit Tower) {
0297       towerHits[idx].fill(Tower);
0298       filled[idx] = 0x1;
0299     }
0300 
0301     void init() {
0302       SimplifiedTower emptyT;
0303       std::fill(towerHits, towerHits + 44, emptyT);
0304       std::fill(filled, filled + 44, 0x0);
0305     }
0306   };
0307 
0308   class InputTowerCluster_pstn {
0309   public:
0310     IEta_t seedIeta = 0;
0311     IPhi_t seedIphi = 0;
0312 
0313     void fill(SimpleTowerHit Tower) {
0314       seedIeta = Tower.towerIeta;
0315       seedIphi = Tower.towerIphi;
0316     }
0317   };
0318 
0319   // INFO : now variables are in GCT precision, they should be in NN precision
0320   // after scaling, i.e. something like ap_fixed<16, 6, AP_TRN, AP_SAT>, when the
0321   // full hls4ml emulation is available
0322   class InputHGCluster {
0323   public:
0324     Pt_t pt;
0325     EtaPhi_t eta;
0326     ShowLen_t showerlength;
0327     ShowLen_t coreshowerlength;
0328     ShapeFeat_t spptot;
0329     ShapeFeat_t szz;
0330     ShapeFeat_t srrtot;
0331     Meanz_t meanz;
0332 
0333     void fill(SimpleHGCluster Cluster) {
0334       pt = Cluster.pt;
0335       eta = Cluster.eta;
0336       showerlength = Cluster.showerlength;
0337       coreshowerlength = Cluster.coreshowerlength;
0338       spptot = Cluster.spptot;
0339       szz = Cluster.szz;
0340       srrtot = Cluster.srrtot;
0341       meanz = Cluster.meanz;
0342     }
0343   };
0344 
0345   l1t::Tau MakeTauCandidate(bool isBarrel,
0346                             int clNxMIdx,
0347                             std::vector<tensorflow::Tensor> outputsIdent,
0348                             std::vector<tensorflow::Tensor> outputsCalib,
0349                             std::vector<InputTowerCluster_pstn> clustersNxM_pstn);
0350 };
0351 
0352 /*
0353 ████████ ██   ██ ██████     ████████  █████  ██   ██ ███    ███ ██ ███    ██  █████  ████████  ██████  ██████  
0354    ██    ██   ██ ██            ██    ██   ██ ██   ██ ████  ████ ██ ████   ██ ██   ██    ██    ██    ██ ██   ██ 
0355    ██    ███████ █████         ██    ███████ ██   ██ ██ ████ ██ ██ ██ ██  ██ ███████    ██    ██    ██ ██████  
0356    ██    ██   ██ ██            ██    ██   ██ ██   ██ ██  ██  ██ ██ ██  ██ ██ ██   ██    ██    ██    ██ ██   ██ 
0357    ██    ██   ██ ██████        ██    ██   ██ ███████ ██      ██ ██ ██   ████ ██   ██    ██     ██████  ██    ██
0358 */
0359 
0360 std::unique_ptr<NNmodels_GlobalCache> L1NNCaloTauEmulator::initializeGlobalCache(const edm::ParameterSet& iConfig) {
0361   edm::LogInfo("Initialization") << "Init NN models Global Cache " << std::endl;
0362 
0363   std::unique_ptr<NNmodels_GlobalCache> GlobalCache(new NNmodels_GlobalCache);
0364 
0365   GlobalCache->CNNmodel_CB_path = iConfig.getParameter<std::string>("CNNmodel_CB_path");
0366   GlobalCache->DNNident_CB_path = iConfig.getParameter<std::string>("DNNident_CB_path");
0367   GlobalCache->DNNcalib_CB_path = iConfig.getParameter<std::string>("DNNcalib_CB_path");
0368   GlobalCache->CNNmodel_CE_path = iConfig.getParameter<std::string>("CNNmodel_CE_path");
0369   GlobalCache->DNNident_CE_path = iConfig.getParameter<std::string>("DNNident_CE_path");
0370   GlobalCache->DNNcalib_CE_path = iConfig.getParameter<std::string>("DNNcalib_CE_path");
0371   GlobalCache->FeatScaler_CE_path = iConfig.getParameter<std::string>("FeatScaler_CE_path");
0372 
0373   // Create sessions for Tensorflow inferece
0374   (GlobalCache->CNNmodel_CB) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->CNNmodel_CB_path)).fullPath());
0375   (GlobalCache->CNNmodel_CBsession) = tensorflow::createSession((GlobalCache->CNNmodel_CB));
0376 
0377   (GlobalCache->DNNident_CB) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNident_CB_path)).fullPath());
0378   (GlobalCache->DNNident_CBsession) = tensorflow::createSession((GlobalCache->DNNident_CB));
0379 
0380   (GlobalCache->DNNcalib_CB) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNcalib_CB_path)).fullPath());
0381   (GlobalCache->DNNcalib_CBsession) = tensorflow::createSession((GlobalCache->DNNcalib_CB));
0382 
0383   (GlobalCache->CNNmodel_CE) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->CNNmodel_CE_path)).fullPath());
0384   (GlobalCache->CNNmodel_CEsession) = tensorflow::createSession((GlobalCache->CNNmodel_CE));
0385 
0386   (GlobalCache->DNNident_CE) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNident_CE_path)).fullPath());
0387   (GlobalCache->DNNident_CEsession) = tensorflow::createSession((GlobalCache->DNNident_CE));
0388 
0389   (GlobalCache->DNNcalib_CE) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNcalib_CE_path)).fullPath());
0390   (GlobalCache->DNNcalib_CEsession) = tensorflow::createSession((GlobalCache->DNNcalib_CE));
0391 
0392   // Read features scaler
0393   boost::property_tree::read_json(edm::FileInPath((GlobalCache->FeatScaler_CE_path)).fullPath(),
0394                                   (GlobalCache->FeatScaler_CE));
0395 
0396   return GlobalCache;
0397 }
0398 
0399 // ----Constructor and Destructor -----
0400 L1NNCaloTauEmulator::L1NNCaloTauEmulator(const edm::ParameterSet& iConfig, const NNmodels_GlobalCache* globalCache)
0401     : l1TowersToken(consumes<l1tp2::CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("l1CaloTowers"))),
0402       hgcalTowersToken(consumes<l1t::HGCalTowerBxCollection>(iConfig.getParameter<edm::InputTag>("hgcalTowers"))),
0403 
0404       HGClusterToken(
0405           consumes<l1t::HGCalMulticlusterBxCollection>(iConfig.getParameter<edm::InputTag>("HgcalClusters"))),
0406       scenario(UseEmInterp::No),
0407       preEmId(iConfig.getParameter<std::string>("preEmId")),
0408       VsPuId(iConfig.getParameter<edm::ParameterSet>("VsPuId")),
0409 
0410       EcalEtMinForClustering(iConfig.getParameter<double>("EcalEtMinForClustering")),
0411       HcalEtMinForClustering(iConfig.getParameter<double>("HcalEtMinForClustering")),
0412       EtMinForSeeding(iConfig.getParameter<double>("EtMinForSeeding")),
0413       EtaRestriction(iConfig.getParameter<double>("EtaRestriction")),
0414       CB_CE_split(iConfig.getParameter<double>("CB_CE_split")),
0415       PuidThr(iConfig.getParameter<double>("PuidThr")),
0416 
0417       IdWp90_CB(iConfig.getParameter<double>("IdWp90_CB")),
0418       IdWp95_CB(iConfig.getParameter<double>("IdWp95_CB")),
0419       IdWp99_CB(iConfig.getParameter<double>("IdWp99_CB")),
0420 
0421       IdWp90_CE(iConfig.getParameter<double>("IdWp90_CE")),
0422       IdWp95_CE(iConfig.getParameter<double>("IdWp95_CE")),
0423       IdWp99_CE(iConfig.getParameter<double>("IdWp99_CE")) {
0424   // Initialize HGCAL BDTs
0425   if (!VsPuId.method().empty()) {
0426     VsPuId.prepareTMVA();
0427   }
0428 
0429   intPuidThr = apintQuantizer(PuidThr, PUID_LSB, PUID_W);
0430   intEtaRestriction = apintQuantizer(EtaRestriction, IETAPHI_LSB, IETAPHI_W);
0431   intCB_CE_split = apintQuantizer(CB_CE_split, IETAPHI_LSB, IETAPHI_W) + 1;
0432 
0433   // Create produced outputs
0434   produces<BXVector<l1t::Tau>>("L1NNCaloTauCollectionBXV");
0435 
0436   // Settings output
0437   edm::LogInfo("Settings") << "EtaRestriction = " << EtaRestriction << " (" << intEtaRestriction << ")"
0438                            << " , CB_CE_split = " << CB_CE_split << "(" << intCB_CE_split
0439                            << ") , EtMinForSeeding = " << EtMinForSeeding
0440                            << " , HcalTpEtMin = " << HcalEtMinForClustering
0441                            << " , EcalTpEtMin = " << EcalEtMinForClustering << " , PuidThr = " << PuidThr << "("
0442                            << intPuidThr << ")" << std::endl;
0443 }
0444 
0445 void L1NNCaloTauEmulator::produce(edm::Event& iEvent, const edm::EventSetup& eSetup) {
0446   // Output collection
0447   std::unique_ptr<BXVector<l1t::Tau>> L1NNCaloTauCollectionBXV(new l1t::TauBxCollection);
0448 
0449   // Create and Fill collection of all calotowers and their attributes
0450   std::vector<SimpleTowerHit> l1CaloTowers;
0451 
0452   iEvent.getByToken(l1TowersToken, l1CaloTowerHandle);
0453   int warnings = 0;
0454   for (auto& hit : *l1CaloTowerHandle.product()) {
0455     // Skip this weird towers and store warning
0456     if (hit.towerIEta() == -1016 && hit.towerIPhi() == -962) {
0457       warnings += 1;
0458       continue;
0459     }
0460 
0461     SimpleTowerHit l1Hit;
0462     l1Hit.isBarrel = 0x1;
0463     l1Hit.l1egTowerEt = apfixedQuantizer(hit.l1egTowerEt(), PTET_LSB, ET_W);
0464     l1Hit.towerEm = apfixedQuantizer(hit.ecalTowerEt(), PTET_LSB, ET_W);
0465     l1Hit.towerHad = apfixedQuantizer(hit.hcalTowerEt(), PTET_LSB, ET_W);
0466     l1Hit.towerEt = apfixedQuantizer(hit.ecalTowerEt() + hit.hcalTowerEt() + hit.l1egTowerEt(), PTET_LSB, ET_W);
0467     l1Hit.towerIeta = hit.towerIEta();
0468     l1Hit.towerIphi = hit.towerIPhi();
0469 
0470     l1CaloTowers.push_back(l1Hit);
0471   }
0472   if (warnings != 0) {
0473     edm::LogWarning("BrokenTowers") << " ** WARNING : FOUND " << warnings
0474                                     << " TOWERS WITH towerIeta=-1016 AND towerIphi=-962" << std::endl;
0475   }
0476 
0477   iEvent.getByToken(hgcalTowersToken, hgcalTowersHandle);
0478   for (auto& hit : *hgcalTowersHandle.product()) {
0479     SimpleTowerHit l1Hit;
0480     l1Hit.isBarrel = 0x0;
0481     l1Hit.l1egTowerEt = 0.0;
0482     l1Hit.towerEm = apfixedQuantizer(hit.etEm(), PTET_LSB, ET_W);
0483     l1Hit.towerHad = apfixedQuantizer(hit.etHad(), PTET_LSB, ET_W);
0484     l1Hit.towerEt = apfixedQuantizer(hit.etEm() + hit.etHad(), PTET_LSB, ET_W);
0485     l1Hit.towerIeta = makeEndcapHwIEta(hit.eta());
0486     l1Hit.towerIphi = makeEndcapHwIPhi(hit.phi());
0487 
0488     l1CaloTowers.push_back(l1Hit);
0489   }
0490 
0491   // Sort the ECAL+HCAL+L1EGs tower sums based on total ET
0492   std::sort(begin(l1CaloTowers), end(l1CaloTowers), [](const SimpleTowerHit& a, SimpleTowerHit& b) {
0493     return a.towerEt > b.towerEt;
0494   });
0495 
0496   // Create and Fill the collection of 3D clusters and their attributes
0497   std::vector<SimpleHGCluster> AllHGClusters;
0498   iEvent.getByToken(HGClusterToken, HGClusterHandle);
0499 
0500   for (auto cl3dIt = HGClusterHandle->begin(0); cl3dIt != HGClusterHandle->end(0); ++cl3dIt) {
0501     auto& cl3d = *cl3dIt;
0502 
0503     // Implement cl3d PU ID as done in
0504     // https://github.com/cms-sw/cmssw/blob/master/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc#L120
0505     bool isEM = preEmId(*cl3dIt);
0506     l1t::PFCluster cluster(cl3d.pt(), cl3d.eta(), cl3d.phi(), cl3d.hOverE());
0507     if (scenario == UseEmInterp::EmOnly)  // for emID objs, use EM interp as pT and set H = 0
0508     {
0509       if (isEM) {
0510         float pt_new = cl3d.iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
0511         float hoe_new = 0.;
0512         cluster = l1t::PFCluster(pt_new, cl3d.eta(), cl3d.phi(), hoe_new, isEM);
0513       }
0514     } else if (scenario == UseEmInterp::AllKeepHad)  // for all objs, replace EM part with EM interp, preserve H
0515     {
0516       float had_old = cl3d.pt() - cluster.emEt();
0517       float em_new = cl3d.iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
0518       float pt_new = had_old + em_new;
0519       float hoe_new = em_new > 0 ? (had_old / em_new) : -1;
0520       cluster = l1t::PFCluster(pt_new, cl3d.eta(), cl3d.phi(), hoe_new, isEM);
0521     } else if (scenario == UseEmInterp::AllKeepTot)  // for all objs, replace EM part with EM interp, preserve pT
0522     {
0523       float em_new = cl3d.iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
0524       float hoe_new = em_new > 0 ? (cl3d.pt() / em_new - 1) : -1;
0525       cluster = l1t::PFCluster(cl3d.pt(), cl3d.eta(), cl3d.phi(), hoe_new, isEM);
0526     }
0527 
0528     float idScore = -1.;
0529     if (!VsPuId.method().empty()) {
0530       idScore = VsPuId.passID(*cl3dIt, cluster);
0531       idScore = cluster.egVsPUMVAOut();
0532     }
0533 
0534     float eta_hgcalCoord = correctInputEtaCl3d(cl3d.eta());
0535     float meanz_hgcalCoord = correctInputMeanzCl3d(cl3d.zBarycenter());
0536 
0537     SimpleHGCluster HGCluster;
0538     HGCluster.pt = apfixedQuantizer(cl3d.pt(), PTET_LSB, PT_W);
0539     HGCluster.eta = apintQuantizer(eta_hgcalCoord, ETAPHI_LSB, ETAPHI_W);
0540     HGCluster.phi = apintQuantizer(cl3d.phi(), ETAPHI_LSB, ETAPHI_W);
0541     HGCluster.showerlength = cl3d.showerLength();
0542     HGCluster.coreshowerlength = cl3d.coreShowerLength();
0543     HGCluster.spptot = apintQuantizer(cl3d.sigmaPhiPhiTot(), SHAPEFEAT_LSB, SHAPEFEAT_W);
0544     HGCluster.szz = apintQuantizer(cl3d.sigmaZZ(), SZZ_LSB, SHAPEFEAT_W);
0545     HGCluster.srrtot = apintQuantizer(cl3d.sigmaRRTot(), SHAPEFEAT_LSB, SHAPEFEAT_W);
0546     HGCluster.meanz = apintQuantizer(meanz_hgcalCoord, MEANZ_LSB, MEANZ_W);
0547     HGCluster.PUid = apintQuantizer(idScore, PUID_LSB, PUID_W);
0548 
0549     AllHGClusters.push_back(HGCluster);
0550   }
0551 
0552   // Order the collection in pt (the input to the GCT will be pt ordered)
0553   std::sort(begin(AllHGClusters), end(AllHGClusters), [](const SimpleHGCluster& a, SimpleHGCluster& b) {
0554     return a.pt > b.pt;
0555   });
0556 
0557   /*
0558   // END OF SOFTWARE PRECISION SECTION
0559   // up to here treated inputs from simulation with SW precision
0560   // to massage them into the HW precision varibales as they are
0561   // forseen (roughly) to be available at the GCT Sum card level
0562   // ------------------------------------------------------------- */
0563 
0564   // Make NxM TowerClusters and HGClusters collections for TauMinator
0565   std::vector<InputTowerCluster> l1TowerClustersNxM_CB;
0566   std::vector<InputTowerCluster_pstn> l1TowerClustersNxM_CB_pstn;
0567   std::vector<InputTowerCluster> l1TowerClustersNxM_CE;
0568   std::vector<InputTowerCluster_pstn> l1TowerClustersNxM_CE_pstn;
0569   std::vector<InputHGCluster> HGClusters;
0570 
0571   // Supporting collection of endcap clusters before cl3d matching
0572   std::vector<InputTowerCluster> AllL1TowerClustersNxM_CE;
0573   std::vector<InputTowerCluster_pstn> AllL1TowerClustersNxM_CE_pstn;
0574 
0575   int Nclusters_CB = 0;
0576   int AllNclusters_CE = 0;
0577   bool caloTauSeedingFinished = false;
0578   // Loop for seeding of clNxM objects
0579   while (!caloTauSeedingFinished) {
0580     InputTowerCluster clNxM;
0581     clNxM.init();
0582     InputTowerCluster_pstn clNxM_pstn;
0583     bool seeded = false;
0584 
0585     for (auto& l1CaloTower : l1CaloTowers) {
0586       // Skip seeding in towers that would make the cluster extend in HF
0587       // Skip l1CaloTowers which are already used by this clusters' mask
0588       if (ap_abs(l1CaloTower.towerIeta) > Eta_limit || ap_abs(l1CaloTower.towerIeta) > intEtaRestriction ||
0589           l1CaloTower.stale4seed) {
0590         continue;
0591       }
0592 
0593       // If not seded do the seeding
0594       if (!seeded) {
0595         // The leading unused tower has ET < min, stop jet clustering
0596         if (l1CaloTower.towerEt < EtMinForSeeding) {
0597           caloTauSeedingFinished = true;
0598           continue;
0599         }
0600 
0601         clNxM.fill(seedIdx, l1CaloTower);
0602         clNxM_pstn.fill(l1CaloTower);
0603         if (l1CaloTower.isBarrel) {
0604           clNxM.barrelSeeded = 0x1;
0605         }
0606 
0607         l1CaloTower.stale4seed = 0x1;
0608         l1CaloTower.stale = 0x1;
0609         seeded = true;
0610 
0611         continue;
0612       }
0613 
0614       dIEtaPhi_t d_iEta = tower_dIEta(l1CaloTower.towerIeta, clNxM_pstn.seedIeta);
0615       dIEtaPhi_t d_iPhi = dPhi<dIEtaPhi_t, IPhi_t>(l1CaloTower.towerIphi, clNxM_pstn.seedIphi);
0616 
0617       // Stale tower for seeding if it would lead to overalp between clusters
0618       if ((ap_abs(d_iEta) <= IEta_dim - 1 && ap_abs(d_iPhi) <= IPhi_dim - 1)) {
0619         l1CaloTower.stale4seed = 0x1;
0620       }
0621 
0622     }  // End for loop over TPs
0623 
0624     // Pushback seeds split in barrel and endcap
0625     if (seeded) {
0626       if (ap_abs(clNxM_pstn.seedIeta) <= intCB_CE_split) {
0627         l1TowerClustersNxM_CB.push_back(clNxM);
0628         l1TowerClustersNxM_CB_pstn.push_back(clNxM_pstn);
0629         Nclusters_CB++;
0630       } else {
0631         AllL1TowerClustersNxM_CE.push_back(clNxM);
0632         AllL1TowerClustersNxM_CE_pstn.push_back(clNxM_pstn);
0633         AllNclusters_CE++;
0634       }
0635     }
0636 
0637   }  // End while loop of TowerClusters seeding
0638 
0639   // Loop for barrel NxM TowerClusters clustering starting from the seeds
0640   for (int clNxMIdx = 0; clNxMIdx < Nclusters_CB; clNxMIdx++) {
0641     for (auto& l1CaloTower : l1CaloTowers) {
0642       // Skip l1CaloTowers which are already used
0643       if (l1CaloTower.stale) {
0644         continue;
0645       }
0646 
0647       dIEtaPhi_t d_iEta = tower_dIEta(l1CaloTower.towerIeta, l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIeta);
0648       dIEtaPhi_t d_iPhi =
0649           dPhi<dIEtaPhi_t, IPhi_t>(l1CaloTower.towerIphi, l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIphi);
0650       int hitIdx = d_iEta * 9 + d_iPhi + seedIdx;
0651 
0652       // Cluster all towers in a NxM towers mask
0653       if ((ap_abs(d_iEta) <= (IEta_dim - 1) / 2 && ap_abs(d_iPhi) <= (IPhi_dim - 1) / 2)) {
0654         l1TowerClustersNxM_CB[clNxMIdx].fill(hitIdx, l1CaloTower);
0655         l1CaloTower.stale = 0x1;
0656       }
0657 
0658     }  // End for loop over TPs
0659 
0660   }  // End while loop of barrel TowerClusters creation
0661 
0662   // In the endcap cross-loop over clNxM and cl3d to match them
0663   // (we can do it before full clustering just using the seed info)
0664   int Nclusters_CE = 0;
0665   for (int clNxMIdx = 0; clNxMIdx < AllNclusters_CE; clNxMIdx++) {
0666     bool matched = false;
0667     for (auto& HGCluster : AllHGClusters) {
0668       // In case the clNxM or HGCluster have already been matched just continue through the list to the end
0669       // only use cl3ds above 4GeV and above -0.10 pu id
0670       if (matched || HGCluster.stale || HGCluster.pt < Pt_t(4.) || HGCluster.PUid < intPuidThr) {
0671         continue;
0672       }
0673 
0674       dEtaPhi_t d_iEta = tw2cl_dEta(HGCluster.eta, AllL1TowerClustersNxM_CE_pstn[clNxMIdx].seedIeta);
0675       dEtaPhi_t d_iPhi = tw2cl_dPhi(HGCluster.phi, AllL1TowerClustersNxM_CE_pstn[clNxMIdx].seedIphi);
0676       matched = d_iEta * d_iEta + d_iPhi * d_iPhi < R2cone;
0677 
0678       if (matched) {
0679         HGCluster.stale = 0x1;
0680         InputHGCluster cl3d;
0681         cl3d.fill(HGCluster);
0682         HGClusters.push_back(cl3d);
0683         l1TowerClustersNxM_CE.push_back(AllL1TowerClustersNxM_CE[clNxMIdx]);
0684         l1TowerClustersNxM_CE_pstn.push_back(AllL1TowerClustersNxM_CE_pstn[clNxMIdx]);
0685         Nclusters_CE++;
0686       }
0687 
0688     }  // End for loop over cl3ds
0689 
0690   }  // End for loop over clNxM
0691 
0692   // Loop for endcap matched NxM TowerClusters clustering starting from the seeds just found
0693   for (int clNxMIdx = 0; clNxMIdx < Nclusters_CE; clNxMIdx++) {
0694     for (auto& l1CaloTower : l1CaloTowers) {
0695       // Skip l1CaloTowers which are already used
0696       if (l1CaloTower.stale) {
0697         continue;
0698       }
0699 
0700       dIEtaPhi_t d_iEta = tower_dIEta(l1CaloTower.towerIeta, l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIeta);
0701       dIEtaPhi_t d_iPhi =
0702           dPhi<dIEtaPhi_t, IPhi_t>(l1CaloTower.towerIphi, l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIphi);
0703       int hitIdx = d_iEta * 9 + d_iPhi + seedIdx;
0704 
0705       // Cluster all towers in a NxM towers mask
0706       if ((ap_abs(d_iEta) <= (IEta_dim - 1) / 2 && ap_abs(d_iPhi) <= (IPhi_dim - 1) / 2)) {
0707         l1TowerClustersNxM_CE[clNxMIdx].fill(hitIdx, l1CaloTower);
0708         l1CaloTower.stale = 0x1;
0709       }
0710 
0711     }  // End for loop over TPs
0712 
0713   }  // End while loop of barrel TowerClusters creation
0714 
0715   // Barrel TauMinator application
0716   tensorflow::setLogging("2");
0717   int batchSize_CB = (int)(Nclusters_CB);
0718   tensorflow::TensorShape imageShape_CB({batchSize_CB, IEta_dim, IPhi_dim, 2});
0719   tensorflow::TensorShape positionShape_CB({batchSize_CB, 2});
0720   tensorflow::Tensor TowerClusterImage_CB(tensorflow::DT_FLOAT, imageShape_CB);
0721   tensorflow::Tensor TowerClusterPosition_CB(tensorflow::DT_FLOAT, positionShape_CB);
0722 
0723   for (int clNxMIdx = 0; clNxMIdx < Nclusters_CB; clNxMIdx++) {
0724     // Fill inputs for Tensorflow inference
0725     for (int eta = 0; eta < IEta_dim; ++eta) {
0726       for (int phi = 0; phi < IPhi_dim; ++phi) {
0727         int towerIdx = eta * IPhi_dim + phi;
0728         TowerClusterImage_CB.tensor<float, 4>()(clNxMIdx, eta, phi, 0) =
0729             (l1TowerClustersNxM_CB[clNxMIdx].towerHits[towerIdx].l1egTowerEt.to_float() +
0730              l1TowerClustersNxM_CB[clNxMIdx].towerHits[towerIdx].towerEm.to_float());
0731         TowerClusterImage_CB.tensor<float, 4>()(clNxMIdx, eta, phi, 1) =
0732             (l1TowerClustersNxM_CB[clNxMIdx].towerHits[towerIdx].towerHad.to_float());
0733       }
0734     }
0735 
0736     TowerClusterPosition_CB.tensor<float, 2>()(clNxMIdx, 0) = floatIEta(l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIeta);
0737     TowerClusterPosition_CB.tensor<float, 2>()(clNxMIdx, 1) = floatIPhi(l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIphi);
0738   }
0739 
0740   if (batchSize_CB >
0741       0)  // from CMSSW_14_0_X tensorflow does not seem to be able to deal with a tensor of dimension 0 anymore
0742   {
0743     // Apply CNN model
0744     tensorflow::NamedTensorList CNNmodel_CBinputList = {{"TowerClusterImage", TowerClusterImage_CB},
0745                                                         {"TowerClusterPosition", TowerClusterPosition_CB}};
0746     std::vector<tensorflow::Tensor> CNNmodel_CBoutputs;
0747     tensorflow::run((globalCache()->CNNmodel_CBsession),
0748                     CNNmodel_CBinputList,
0749                     {"TauMinator_CB_conv/middleMan/concat"},
0750                     &CNNmodel_CBoutputs);
0751     tensorflow::NamedTensorList DNN_CBinputsList = {{"middleMan", CNNmodel_CBoutputs[0]}};
0752 
0753     // Apply DNN for identification
0754     std::vector<tensorflow::Tensor> DNN_CBoutputsIdent;
0755     tensorflow::run((globalCache()->DNNident_CBsession),
0756                     DNN_CBinputsList,
0757                     {"TauMinator_CB_ident/sigmoid_IDout/Sigmoid"},
0758                     &DNN_CBoutputsIdent);
0759 
0760     // Apply DNN for calibration
0761     std::vector<tensorflow::Tensor> DNN_CBoutputsCalib;
0762     tensorflow::run((globalCache()->DNNcalib_CBsession),
0763                     DNN_CBinputsList,
0764                     {"TauMinator_CB_calib/DNNout/MatMul"},
0765                     &DNN_CBoutputsCalib);
0766 
0767     // Fill the output collection of L1 taus with the barrel candidates
0768     for (int clNxMIdx = 0; clNxMIdx < Nclusters_CB; clNxMIdx++) {
0769       l1t::Tau l1Tau =
0770           MakeTauCandidate(true, clNxMIdx, DNN_CBoutputsIdent, DNN_CBoutputsCalib, l1TowerClustersNxM_CB_pstn);
0771       if (l1Tau.pt() < 0) {
0772         continue;
0773       }
0774       L1NNCaloTauCollectionBXV->push_back(0, l1Tau);
0775     }
0776   }
0777 
0778   // Endcap TauMinator application
0779   int batchSize_CE = (int)(Nclusters_CE);
0780   tensorflow::TensorShape imageShape_CE({batchSize_CE, IEta_dim, IPhi_dim, 2});
0781   tensorflow::TensorShape positionShape_CE({batchSize_CE, 2});
0782   tensorflow::TensorShape cl3dfeatShape_CE({batchSize_CE, 8});
0783   tensorflow::Tensor TowerClusterImage_CE(tensorflow::DT_FLOAT, imageShape_CE);
0784   tensorflow::Tensor TowerClusterPosition_CE(tensorflow::DT_FLOAT, positionShape_CE);
0785   tensorflow::Tensor Cl3dShapeFeatures_CE(tensorflow::DT_FLOAT, cl3dfeatShape_CE);
0786 
0787   for (int clNxMIdx = 0; clNxMIdx < Nclusters_CE; clNxMIdx++) {
0788     // Indexing of cl3ds is the same as the one of clNxMs
0789     InputHGCluster HGClu = HGClusters[clNxMIdx];
0790 
0791     // Fill inputs for Tensorflow inference
0792     for (int eta = 0; eta < IEta_dim; ++eta) {
0793       for (int phi = 0; phi < IPhi_dim; ++phi) {
0794         int towerIdx = eta * IPhi_dim + phi;
0795         TowerClusterImage_CE.tensor<float, 4>()(clNxMIdx, eta, phi, 0) =
0796             (l1TowerClustersNxM_CE[clNxMIdx].towerHits[towerIdx].l1egTowerEt.to_float() +
0797              l1TowerClustersNxM_CE[clNxMIdx].towerHits[towerIdx].towerEm.to_float());
0798         TowerClusterImage_CE.tensor<float, 4>()(clNxMIdx, eta, phi, 1) =
0799             (l1TowerClustersNxM_CE[clNxMIdx].towerHits[towerIdx].towerHad.to_float());
0800       }
0801     }
0802 
0803     TowerClusterPosition_CE.tensor<float, 2>()(clNxMIdx, 0) = floatIEta(l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIeta);
0804     TowerClusterPosition_CE.tensor<float, 2>()(clNxMIdx, 1) = floatIPhi(l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIphi);
0805 
0806     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 0) = inputScaler(HGClu.pt.to_float(), "pt");
0807     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 1) = inputScaler(abs(floatEta(HGClu.eta)), "eta");
0808     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 2) = inputScaler(HGClu.showerlength.to_float(), "showerlength");
0809     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 3) =
0810         inputScaler(HGClu.coreshowerlength.to_float(), "coreshowerlength");
0811     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 4) = inputScaler(floatShape(HGClu.spptot), "spptot");
0812     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 5) = inputScaler(floatSzz(HGClu.szz), "szz");
0813     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 6) = inputScaler(floatShape(HGClu.srrtot), "srrtot");
0814     Cl3dShapeFeatures_CE.tensor<float, 2>()(clNxMIdx, 7) = inputScaler(floatMeanZHgcalCoord(HGClu.meanz), "meanz");
0815   }
0816 
0817   if (batchSize_CE >
0818       0)  // from CMSSW_14_0_X tensorflow does not seem to be able to deal with a tensor of dimension 0 anymore
0819   {
0820     // Apply CNN model
0821     tensorflow::NamedTensorList CNNmodel_CEinputList = {{"TowerClusterImage", TowerClusterImage_CE},
0822                                                         {"TowerClusterPosition", TowerClusterPosition_CE},
0823                                                         {"AssociatedCl3dFeatures", Cl3dShapeFeatures_CE}};
0824     std::vector<tensorflow::Tensor> CNNmodel_CEoutputs;
0825     tensorflow::run((globalCache()->CNNmodel_CEsession),
0826                     CNNmodel_CEinputList,
0827                     {"TauMinator_CE_conv/middleMan/concat"},
0828                     &CNNmodel_CEoutputs);
0829     tensorflow::NamedTensorList DNN_CEinputsList = {{"middleMan", CNNmodel_CEoutputs[0]}};
0830 
0831     // Apply DNN for identification
0832     std::vector<tensorflow::Tensor> DNN_CEoutputsIdent;
0833     tensorflow::run((globalCache()->DNNident_CEsession),
0834                     DNN_CEinputsList,
0835                     {"TauMinator_CE_ident/sigmoid_IDout/Sigmoid"},
0836                     &DNN_CEoutputsIdent);
0837 
0838     // Apply DNN for calibration
0839     std::vector<tensorflow::Tensor> DNN_CEoutputsCalib;
0840     tensorflow::run((globalCache()->DNNcalib_CEsession),
0841                     DNN_CEinputsList,
0842                     {"TauMinator_CE_calib/LIN_DNNout/Relu"},
0843                     &DNN_CEoutputsCalib);
0844 
0845     // Fill the output collection of L1 taus with the endcap candidates
0846     for (int clNxMIdx = 0; clNxMIdx < Nclusters_CE; clNxMIdx++) {
0847       l1t::Tau l1Tau =
0848           MakeTauCandidate(false, clNxMIdx, DNN_CEoutputsIdent, DNN_CEoutputsCalib, l1TowerClustersNxM_CE_pstn);
0849       if (l1Tau.pt() < 0) {
0850         continue;
0851       }
0852       L1NNCaloTauCollectionBXV->push_back(0, l1Tau);
0853     }
0854   }
0855 
0856   // Fill output
0857   iEvent.put(std::move(L1NNCaloTauCollectionBXV), "L1NNCaloTauCollectionBXV");
0858 
0859 }  // End of produce function
0860 
0861 template <class outPrecision, class inPrecision>
0862 outPrecision L1NNCaloTauEmulator::dPhi(inPrecision iPhi_1, inPrecision iPhi_2) {
0863   outPrecision dphi = iPhi_1 - iPhi_2;
0864 
0865   outPrecision dphi0 = dphi > outPrecision(INTPHI_PI) ? outPrecision(dphi - INTPHI_2PI) : dphi;
0866   outPrecision dphi1 = dphi <= outPrecision(-INTPHI_PI) ? outPrecision(dphi + INTPHI_2PI) : dphi;
0867 
0868   outPrecision result = dphi > outPrecision(0) ? dphi0 : dphi1;
0869 
0870   return result;
0871 }
0872 
0873 L1NNCaloTauEmulator::dIEtaPhi_t L1NNCaloTauEmulator::tower_dIEta(IEta_t iEta_1, IEta_t iEta_2) {
0874   ap_int<12> mult = iEta_1 * iEta_2;
0875   dIEtaPhi_t result = iEta_1 - iEta_2;
0876   if (mult < 0) {
0877     result = iEta_1 > 0 ? result - 1 : result + 1;
0878   }
0879 
0880   return result;
0881 }
0882 
0883 L1NNCaloTauEmulator::dEtaPhi_t L1NNCaloTauEmulator::tw2cl_dPhi(EtaPhi_t iPhi_1, IPhi_t iPhi_2) {
0884   EtaPhi_t shiftediPhi_2 = iPhi_2 <= IPhi_t(36) ? EtaPhi_t(iPhi_2) : EtaPhi_t(iPhi_2 - INTPHI_2PI + 1);
0885 
0886   EtaPhi_t fineiPhi_2 = shiftediPhi_2 * (IETAPHI_LSB / ETAPHI_LSB);
0887   // subrtaction of half rescaling corrects from edge to center of tower
0888   fineiPhi_2 = fineiPhi_2 > EtaPhi_t(0) ? EtaPhi_t(fineiPhi_2 - (IETAPHI_LSB / ETAPHI_LSB) / 2)
0889                                         : EtaPhi_t(fineiPhi_2 + (IETAPHI_LSB / ETAPHI_LSB) / 2);
0890 
0891   return dPhi<dEtaPhi_t, EtaPhi_t>(iPhi_1, fineiPhi_2);
0892 }
0893 
0894 L1NNCaloTauEmulator::dEtaPhi_t L1NNCaloTauEmulator::tw2cl_dEta(EtaPhi_t iEta_1, IEta_t iEta_2) {
0895   // change from hgcal frame to barrel-centered frame
0896   EtaPhi_t framechangeCl3d = 303;  // 303*pi/720 = 1.322
0897   iEta_1 = iEta_1 > EtaPhi_t(0) ? EtaPhi_t(iEta_1 + framechangeCl3d) : EtaPhi_t(iEta_1 - framechangeCl3d);
0898 
0899   // the actual depth is 330 but 329 corrects for 0.0808 tower
0900   EtaPhi_t barrelEtaDepth = 329;
0901   EtaPhi_t fineiEta_2 = barrelEtaDepth + (iEta_2 - IETAHGCAL_OFFSET) * (IETAHGCAL_LSB / ETAPHI_LSB);
0902 
0903   return iEta_1 - fineiEta_2;
0904 }
0905 
0906 L1NNCaloTauEmulator::IEta_t L1NNCaloTauEmulator::makeEndcapHwIEta(float eta) {
0907   IEta_t ieta = floor(eta / IETAHGCAL_LSB);
0908   // +1 because flooring gets it 1 unit lower when negative
0909   ieta = ieta < IEta_t(0) ? IEta_t(ieta + 1) : ieta;
0910 
0911   return ieta;
0912 }
0913 
0914 L1NNCaloTauEmulator::IPhi_t L1NNCaloTauEmulator::makeEndcapHwIPhi(float phi) {
0915   phi = phi < 0 ? phi + 2 * M_PI : phi;
0916 
0917   // +1 because tower 0 does not exist
0918   return floor(phi / IETAPHI_LSB) + 1;
0919 }
0920 
0921 template <int W>
0922 ap_int<W> L1NNCaloTauEmulator::ap_abs(ap_int<W> x) {
0923   ap_int<W> result;
0924   if (x < 0) {
0925     result = -x;
0926   } else {
0927     result = x;
0928   }
0929 
0930   return result;
0931 }
0932 
0933 template <int W, int I, ap_q_mode _AP_Q, ap_o_mode _AP_O>
0934 ap_ufixed<W, I> L1NNCaloTauEmulator::ap_abs(ap_fixed<W, I, _AP_Q, _AP_O> x) {
0935   ap_ufixed<W, I> result;
0936   if (x < 0) {
0937     result = -x;
0938   } else {
0939     result = x;
0940   }
0941 
0942   return result;
0943 }
0944 
0945 float L1NNCaloTauEmulator::apfixedQuantizer(float inputF, float LSB, int nbits) {
0946   return min(floor(inputF / LSB), float(pow(2, nbits) - 1)) * LSB;
0947 }
0948 
0949 int L1NNCaloTauEmulator::apintQuantizer(float inputF, float LSB, int nbits) {
0950   return min(floor(inputF / LSB), float(pow(2, nbits) - 1));
0951 }
0952 
0953 float L1NNCaloTauEmulator::inputScaler(float inputF, std::string feature) {
0954   float mean = (globalCache()->FeatScaler_CE).get_child(feature).get<float>("mean");
0955   float std = (globalCache()->FeatScaler_CE).get_child(feature).get<float>("std");
0956 
0957   return (inputF - mean) / std;
0958 }
0959 
0960 float L1NNCaloTauEmulator::correctInputEtaCl3d(float eta) {
0961   return eta > 0 ? eta - ETAHGCAL_OFFSET : eta + ETAHGCAL_OFFSET;
0962 }
0963 
0964 float L1NNCaloTauEmulator::correctInputMeanzCl3d(float meanz) { return CM2MM * (abs(meanz) - MEANZ_OFFSET); }
0965 
0966 float L1NNCaloTauEmulator::floatIEta(IEta_t eta) {
0967   // transform eta of towers from integer to float, correcting for different barrel/endcap LSB
0968   float feta;
0969   if (abs(eta) > IETAHGCAL_OFFSET) {
0970     if (eta > 0) {
0971       feta = IETAHGCAL_OFFSET * IETAPHI_LSB - (IETAHGCAL_LSB - IETAHGCAL_LSBp) +
0972              (eta.to_float() - IETAHGCAL_OFFSET) * IETAHGCAL_LSB;
0973     } else {
0974       feta = -IETAHGCAL_OFFSET * IETAPHI_LSB + (IETAHGCAL_LSB - IETAHGCAL_LSBp) +
0975              (eta.to_float() + IETAHGCAL_OFFSET) * IETAHGCAL_LSB;
0976     }
0977   } else {
0978     feta = eta.to_float() * IETAPHI_LSB;
0979   }
0980 
0981   // shift by half a tower to consider the tower center instead of the edge
0982   return feta > 0 ? feta - IETAPHI_LSB / 2 : feta + IETAPHI_LSB / 2;
0983 }
0984 
0985 float L1NNCaloTauEmulator::floatIPhi(IPhi_t phi) {
0986   float fphi = phi.to_float();
0987   // add 2pi + 1 because tower 0 does not exist
0988   fphi = fphi > INTPHI_PI ? fphi - INTPHI_2PI + 1 : fphi;
0989   fphi *= IETAPHI_LSB;
0990 
0991   // shift by half a tower to consider the tower center instead of the edge
0992   return fphi > 0 ? fphi - IETAPHI_LSB / 2 : fphi + IETAPHI_LSB / 2;
0993 }
0994 
0995 l1t::Tau L1NNCaloTauEmulator::MakeTauCandidate(
0996     bool isBarrel,
0997     int clNxMIdx,
0998     std::vector<tensorflow::Tensor> outputsIdent,
0999     std::vector<tensorflow::Tensor> outputsCalib,
1000     std::vector<L1NNCaloTauEmulator::InputTowerCluster_pstn> clustersNxM_pstn) {
1001   int seedIeta = clustersNxM_pstn[clNxMIdx].seedIeta;
1002   int seedIphi = clustersNxM_pstn[clNxMIdx].seedIphi;
1003 
1004   if (seedIeta > intEtaRestriction) {
1005     return l1t::Tau(reco::Candidate::PolarLorentzVector(-1, 0, 0, 0), -1, 0, 0, 0, 0);
1006     ;
1007   }
1008 
1009   float tau_IDscore = outputsIdent[0].matrix<float>()(0, clNxMIdx);
1010   float tau_calibPt = outputsCalib[0].matrix<float>()(0, clNxMIdx);
1011   float tau_eta = floatIEta(seedIeta);
1012   float tau_phi = floatIPhi(seedIphi);
1013 
1014   // Assign increasing quality to higher scoring candidates
1015   int quality = 0;
1016   if (isBarrel) {
1017     // 99% WP
1018     if (tau_IDscore > IdWp99_CB) {
1019       quality = 1;
1020     }
1021     // 95% WP
1022     if (tau_IDscore > IdWp95_CB) {
1023       quality = 2;
1024     }
1025     // 90% WP
1026     if (tau_IDscore > IdWp90_CB) {
1027       quality = 3;
1028     }
1029   } else {
1030     // 99% WP
1031     if (tau_IDscore > IdWp99_CE) {
1032       quality = 1;
1033     }
1034     // 95% WP
1035     if (tau_IDscore > IdWp95_CE) {
1036       quality = 2;
1037     }
1038     // 90% WP
1039     if (tau_IDscore > IdWp90_CE) {
1040       quality = 3;
1041     }
1042   }
1043 
1044   reco::Candidate::PolarLorentzVector tauP4 = reco::Candidate::PolarLorentzVector(tau_calibPt, tau_eta, tau_phi, 0);
1045 
1046   // store ID score multiplied by 10E4 to have good precision even using the Phase1 tau int iso format
1047   // (this is stored just in case for possible additional offline studies)
1048   // tau initialisation =  (p4,    pt,          eta,     phi,     qual,    iso)
1049   l1t::Tau l1Tau = l1t::Tau(tauP4, tau_calibPt, tau_eta, tau_phi, quality, tau_IDscore * 10E4);
1050   l1Tau.setTowerIEta(seedIeta);
1051   l1Tau.setTowerIPhi(seedIphi);
1052 
1053   return l1Tau;
1054 }
1055 
1056 void L1NNCaloTauEmulator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1057   edm::ParameterSetDescription desc;
1058 
1059   desc.add<edm::InputTag>("l1CaloTowers", edm::InputTag("l1tEGammaClusterEmuProducer", "L1CaloTowerCollection"));
1060   desc.add<edm::InputTag>("hgcalTowers", edm::InputTag("l1tHGCalTowerProducer", "HGCalTowerProcessor"));
1061   desc.add<edm::InputTag>("HgcalClusters",
1062                           edm::InputTag("l1tHGCalBackEndLayer2Producer", "HGCalBackendLayer2Processor3DClustering"));
1063 
1064   desc.add<std::string>("preEmId", "hOverE < 0.3 && hOverE >= 0");
1065   {
1066     edm::ParameterSetDescription psd0;
1067     psd0.add<bool>("isPUFilter", true);
1068     psd0.add<std::string>("preselection", "");
1069     psd0.add<std::string>("method", "BDT");
1070     {
1071       edm::ParameterSetDescription vpsd2;
1072       vpsd2.add<std::string>("name");
1073       vpsd2.add<std::string>("value");
1074       std::vector<edm::ParameterSet> temp2;
1075       temp2.reserve(5);
1076       {
1077         edm::ParameterSet temp3;
1078         temp3.addParameter<std::string>("name", "eMax");
1079         temp3.addParameter<std::string>("value", "eMax()");
1080         temp2.push_back(temp3);
1081       }
1082       {
1083         edm::ParameterSet temp3;
1084         temp3.addParameter<std::string>("name", "eMaxOverE");
1085         temp3.addParameter<std::string>("value", "eMax()/energy()");
1086         temp2.push_back(temp3);
1087       }
1088       {
1089         edm::ParameterSet temp3;
1090         temp3.addParameter<std::string>("name", "sigmaPhiPhiTot");
1091         temp3.addParameter<std::string>("value", "sigmaPhiPhiTot()");
1092         temp2.push_back(temp3);
1093       }
1094       {
1095         edm::ParameterSet temp3;
1096         temp3.addParameter<std::string>("name", "sigmaRRTot");
1097         temp3.addParameter<std::string>("value", "sigmaRRTot()");
1098         temp2.push_back(temp3);
1099       }
1100       {
1101         edm::ParameterSet temp3;
1102         temp3.addParameter<std::string>("name", "triggerCells90percent");
1103         temp3.addParameter<std::string>("value", "triggerCells90percent()");
1104         temp2.push_back(temp3);
1105       }
1106       psd0.addVPSet("variables", vpsd2, temp2);
1107     }
1108     psd0.add<std::string>(
1109         "weightsFile", "L1Trigger/Phase2L1ParticleFlow/data/hgcal_egID/Photon_Pion_vs_Neutrino_BDTweights_1116.xml.gz");
1110     psd0.add<std::string>("wp", "-0.10");
1111     desc.add<edm::ParameterSetDescription>("VsPuId", psd0);
1112   }
1113 
1114   desc.add<double>("EcalEtMinForClustering", 0.0);
1115   desc.add<double>("HcalEtMinForClustering", 0.0);
1116   desc.add<double>("EtMinForSeeding", 2.5);
1117   desc.add<double>("EtaRestriction", 2.4);
1118   desc.add<double>("CB_CE_split", 1.55);
1119   desc.add<double>("PuidThr", -0.1);
1120 
1121   desc.add<std::string>("CNNmodel_CB_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/CNNmodel_CB.pb");
1122   desc.add<std::string>("DNNident_CB_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNident_CB.pb");
1123   desc.add<std::string>("DNNcalib_CB_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNcalib_CB.pb");
1124   desc.add<std::string>("CNNmodel_CE_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/CNNmodel_CE.pb");
1125   desc.add<std::string>("DNNident_CE_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNident_CE.pb");
1126   desc.add<std::string>("DNNcalib_CE_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNcalib_CE.pb");
1127   desc.add<std::string>("FeatScaler_CE_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/Cl3dFeatScaler_CE.json");
1128 
1129   desc.add<double>("IdWp90_CB", 0.706);
1130   desc.add<double>("IdWp95_CB", 0.3432);
1131   desc.add<double>("IdWp99_CB", 0.0337);
1132   desc.add<double>("IdWp90_CE", 0.5711);
1133   desc.add<double>("IdWp95_CE", 0.2742);
1134   desc.add<double>("IdWp99_CE", 0.0394);
1135 
1136   desc.add<bool>("DEBUG", false);
1137 
1138   descriptions.add("l1tNNCaloTauEmulator", desc);
1139 }
1140 
1141 DEFINE_FWK_MODULE(L1NNCaloTauEmulator);