diff --git a/L1Trigger/Configuration/python/L1Trigger_EventContent_cff.py b/L1Trigger/Configuration/python/L1Trigger_EventContent_cff.py
index 6399f341f7984..bd35740e2123a 100644
--- a/L1Trigger/Configuration/python/L1Trigger_EventContent_cff.py
+++ b/L1Trigger/Configuration/python/L1Trigger_EventContent_cff.py
@@ -191,6 +191,8 @@ def _appendPhase2Digis(obj):
'keep *_l1tTowerCalibration_*_*',
'keep *_l1tCaloJet_*_*',
'keep *_l1tCaloJetHTT_*_*',
+ 'keep *_l1tNNCaloTauProducer_*_*',
+ 'keep *_l1tNNCaloTauEmulator_*_*',
'keep *_l1tPFClustersFromL1EGClusters_*_*',
'keep *_l1tPFClustersFromCombinedCaloHCal_*_*',
'keep *_l1tPFClustersFromCombinedCaloHF_*_*',
diff --git a/L1Trigger/Configuration/python/SimL1Emulator_cff.py b/L1Trigger/Configuration/python/SimL1Emulator_cff.py
index 99c9ded41520e..1b346c3e35072 100644
--- a/L1Trigger/Configuration/python/SimL1Emulator_cff.py
+++ b/L1Trigger/Configuration/python/SimL1Emulator_cff.py
@@ -94,7 +94,7 @@
from L1Trigger.L1CaloTrigger.l1tPhase2L1CaloEGammaEmulator_cfi import *
_phase2_siml1emulator.add(l1tPhase2L1CaloEGammaEmulator)
-# Barrel and EndCap CaloJet/HT
+# Barrel and EndCap CaloJet/HT/NNCaloTau
# ########################################################################
# ---- Produce the calibrated tower collection combining Barrel, HGCal, HF
from L1Trigger.L1CaloTrigger.l1tTowerCalibrationProducer_cfi import *
@@ -113,6 +113,12 @@
l1tCaloJetHTT = l1tCaloJetHTTProducer.clone(
BXVCaloJetsInputTag = ("L1CaloJet", "CaloJets")
)
+# ---- Produce the NNCaloTau
+from L1Trigger.L1CaloTrigger.l1tNNCaloTauProducer_cfi import *
+_phase2_siml1emulator.add(l1tNNCaloTauProducer)
+
+from L1Trigger.L1CaloTrigger.l1tNNCaloTauEmulator_cfi import *
+_phase2_siml1emulator.add(l1tNNCaloTauEmulator)
_phase2_siml1emulator.add(l1tTowerCalibration)
diff --git a/L1Trigger/L1CaloTrigger/plugins/BuildFile.xml b/L1Trigger/L1CaloTrigger/plugins/BuildFile.xml
index a732a86066576..356e44afd1fb2 100644
--- a/L1Trigger/L1CaloTrigger/plugins/BuildFile.xml
+++ b/L1Trigger/L1CaloTrigger/plugins/BuildFile.xml
@@ -12,7 +12,10 @@
+
+
+
diff --git a/L1Trigger/L1CaloTrigger/plugins/L1NNCaloTauEmulator.cc b/L1Trigger/L1CaloTrigger/plugins/L1NNCaloTauEmulator.cc
new file mode 100644
index 0000000000000..d4abb78f7b032
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/plugins/L1NNCaloTauEmulator.cc
@@ -0,0 +1,1141 @@
+/* -*- C++ -*-
+
+Package: L1CaloTrigger
+Class: L1NNCaloTauEmulator
+Frinedly name: The TauMinator
+
+\class L1NNCaloTauEmulator L1NNCaloTauEmulator.cc
+
+Description:
+Perform firmware-exact emulation of the l1tNNCaloTauProducer
+that implements the NN Calo Tau.
+(Perform reconstruction and identification of tau
+candidates at L1 Trigger with a CNN.)
+
+Implementation:
+The implementation is done forseeing the integration
+of the algorithm in the GCT Sum card. This means that
+the full detector information can be accessed at the same
+time (ECAL, HCAL, HGCAL full eta-phi coverage). This will
+come in the form of arrays of towers and clusters.
+Given that the emulators of the upstream algortihms are
+not fully determined yet, this emulator takes as input
+the simulation-based information, manipulates with software
+precision to pruduce the arrays of towers and clusters as
+they should be available in the GCT sum card.
+Only then the actual fixed point algorithm emulation arrives.
+
+** INFO : THE NNs ARE APPLIED USING THE TENSORFLOW SOFTWARE
+ the implementation of full emulation via hls4ml is ongoing
+ (it has already been shown in other contexts that tensorflow
+ softwrae and full emulation are very close to each other)
+
+Original Author: Jona Motta
+Created: Tue June 7th 2023
+
+*/
+
+#include
+#include
+#include
+
+#include "boost/property_tree/ptree.hpp"
+#include "boost/property_tree/json_parser.hpp"
+
+#include "ap_int.h"
+#include "ap_fixed.h"
+// #include "hls4ml/emulator.h"
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h"
+#include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
+#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
+#include "DataFormats/L1TParticleFlow/interface/PFCluster.h"
+#include "DataFormats/L1THGCal/interface/HGCalTower.h"
+#include "DataFormats/Math/interface/deltaPhi.h"
+#include "DataFormats/L1Trigger/interface/Tau.h"
+
+#include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h"
+#include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h"
+
+#include "L1Trigger/L1THGCal/interface/backend/HGCalTriggerClusterIdentificationBase.h"
+#include "L1Trigger/Phase2L1ParticleFlow/interface/HGC3DClusterEgID.h"
+#include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
+
+#include "PhysicsTools/TensorFlow/interface/TensorFlow.h"
+
+struct NNmodels_GlobalCache {
+ std::string CNNmodel_CB_path;
+ std::string DNNident_CB_path;
+ std::string DNNcalib_CB_path;
+
+ std::string CNNmodel_CE_path;
+ std::string DNNident_CE_path;
+ std::string DNNcalib_CE_path;
+ std::string FeatScaler_CE_path;
+ boost::property_tree::ptree FeatScaler_CE;
+
+ tensorflow::GraphDef* CNNmodel_CB;
+ tensorflow::GraphDef* DNNident_CB;
+ tensorflow::GraphDef* DNNcalib_CB;
+
+ tensorflow::Session* CNNmodel_CBsession;
+ tensorflow::Session* DNNident_CBsession;
+ tensorflow::Session* DNNcalib_CBsession;
+
+ tensorflow::GraphDef* CNNmodel_CE;
+ tensorflow::GraphDef* DNNident_CE;
+ tensorflow::GraphDef* DNNcalib_CE;
+
+ tensorflow::Session* CNNmodel_CEsession;
+ tensorflow::Session* DNNident_CEsession;
+ tensorflow::Session* DNNcalib_CEsession;
+};
+
+class L1NNCaloTauEmulator : public edm::stream::EDProducer> {
+public:
+ explicit L1NNCaloTauEmulator(const edm::ParameterSet&, const NNmodels_GlobalCache*);
+
+ static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+ static std::unique_ptr initializeGlobalCache(const edm::ParameterSet&);
+ static void globalEndJob(const NNmodels_GlobalCache*){/*do nothing*/};
+
+private:
+ // ----fixed LSBs, Nbits, scales, and types----
+ static constexpr int INTPHI_PI = 36;
+ static constexpr int INTPHI_2PI = 2 * INTPHI_PI;
+ static constexpr float IETAPHI_LSB = M_PI / INTPHI_PI;
+
+ static constexpr int FINEINTPHI_PI = 720;
+ static constexpr int FINEINTPHI_2PI = 2 * FINEINTPHI_PI;
+ static constexpr float ETAPHI_LSB = M_PI / FINEINTPHI_PI;
+
+ static constexpr float SHAPEFEAT_LSB = 0.0000153; // pow(2, -16)
+ static constexpr float SZZ_LSB = SHAPEFEAT_LSB * 100;
+ static constexpr float ETAHGCAL_OFFSET = 1.321; // inferred from hgcal schematics
+ static constexpr float IETAHGCAL_LSBp = 0.0808; // inferred from simulation
+ static constexpr float IETAHGCAL_LSB = 0.0845; // inferred from simulation
+ static constexpr float PUID_LSB = 0.00390625; // pow(2, -8)
+ static constexpr float MEANZ_OFFSET = 321.05; // inferred from hgcal schematics
+ static constexpr int IETAHGCAL_OFFSET = 17;
+ static constexpr float MEANZ_LSB = 0.5;
+ static constexpr float PTET_LSB = 0.25;
+ static constexpr float CM2MM = 10;
+ static constexpr int R2cone = 0.25 / ETAPHI_LSB / ETAPHI_LSB;
+
+ static constexpr int SHAPEFEAT_W = 16; // maximum forseen per shape
+ static constexpr int DETAPHI_W = 12;
+ static constexpr int DIETAPHI_W = 8;
+ static constexpr int IETAPHI_W = 7;
+ static constexpr int SHOWLEN_W = 6;
+ static constexpr int ETAPHI_W = 11; // precision going to correlator
+ static constexpr int MEANZ_W = 12;
+ static constexpr int PUID_W = 9;
+
+ static constexpr int PT_W = 14;
+ static constexpr int PT_I = 12;
+ static constexpr int ET_W = 10;
+ static constexpr int ET_I = 8;
+
+ // forseen precision of the HLS4ML emulation of the NNs
+ static constexpr int CALIBPT_W = 10;
+ static constexpr int CALIBPT_I = 9;
+ static constexpr int ID_W = 8;
+ static constexpr int ID_I = 1;
+
+ typedef ap_ufixed Pt_t;
+ typedef ap_ufixed Et_t;
+
+ typedef ap_ufixed CalibPt_t;
+ typedef ap_ufixed Id_t;
+
+ typedef ap_uint ShapeFeat_t;
+ typedef ap_int dIEtaPhi_t;
+ typedef ap_int dEtaPhi_t;
+ typedef ap_uint ShowLen_t;
+ typedef ap_int EtaPhi_t;
+ typedef ap_uint IPhi_t;
+ typedef ap_int IEta_t;
+ typedef ap_uint Meanz_t;
+ typedef ap_int PUid_t;
+
+ // ----fixed dimensions of tower clusters----
+ const int seedIdx = 22;
+ const int IEta_dim = 5;
+ const int IPhi_dim = 9;
+ const int Eta_limit = 33;
+
+ //----edm control---
+ void produce(edm::Event&, const edm::EventSetup&) override;
+
+ //----private functions----
+ template
+ outPrecision dPhi(inPrecision iPhi_1, inPrecision iPhi_2);
+ dIEtaPhi_t tower_dIEta(IEta_t iEta_1, IEta_t iEta_2);
+ dEtaPhi_t tw2cl_dPhi(EtaPhi_t iPhi_1, IPhi_t iPhi_2);
+ dEtaPhi_t tw2cl_dEta(EtaPhi_t iEta_1, IEta_t iEta_2);
+ IEta_t makeEndcapHwIEta(float eta);
+ IPhi_t makeEndcapHwIPhi(float phi);
+ float apfixedQuantizer(float inputF, float LSB, int nbits);
+ int apintQuantizer(float inputF, float LSB, int nbits);
+ float inputScaler(float inputF, std::string feature);
+ float correctInputEtaCl3d(float eta);
+ float correctInputMeanzCl3d(float meanz);
+
+ inline float floatPt(Pt_t pt) { return pt.to_float(); }
+ inline float floatEt(Et_t et) { return et.to_float(); }
+ inline float floatEta(EtaPhi_t eta) { return eta.to_float() * ETAPHI_LSB; }
+ inline float floatPhi(EtaPhi_t phi) { return phi.to_float() * ETAPHI_LSB; }
+ inline float floatShape(ShapeFeat_t shape) { return shape.to_float() * SHAPEFEAT_LSB; };
+ inline float floatSzz(ShapeFeat_t szz) { return szz.to_float() * SZZ_LSB; };
+ inline float floatMeanZ(Meanz_t meanz) { return meanz.to_float() * MEANZ_LSB + MEANZ_OFFSET; };
+ inline float floatMeanZHgcalCoord(Meanz_t meanz) { return meanz.to_float() * MEANZ_LSB; };
+ inline float floatPuId(PUid_t pu) { return pu.to_float() * PUID_LSB; };
+ float floatIEta(IEta_t eta);
+ float floatIPhi(IPhi_t phi);
+
+ template
+ ap_int ap_abs(ap_int x);
+ template
+ ap_ufixed ap_abs(ap_fixed x);
+
+ //----tokens and handles----
+ edm::EDGetTokenT l1TowersToken;
+ edm::Handle l1CaloTowerHandle;
+
+ edm::EDGetToken hgcalTowersToken;
+ edm::Handle hgcalTowersHandle;
+
+ edm::EDGetTokenT HGClusterToken;
+ edm::Handle HGClusterHandle;
+
+ //----private variables----
+ enum class UseEmInterp { No, EmOnly, AllKeepHad, AllKeepTot };
+ UseEmInterp scenario;
+ StringCutObjectSelector preEmId;
+ l1tpf::HGC3DClusterEgID VsPuId;
+
+ double EcalEtMinForClustering;
+ double HcalEtMinForClustering;
+ double EtMinForSeeding;
+ double EtaRestriction;
+ double CB_CE_split;
+ double PuidThr;
+
+ double IdWp90_CB;
+ double IdWp95_CB;
+ double IdWp99_CB;
+
+ double IdWp90_CE;
+ double IdWp95_CE;
+ double IdWp99_CE;
+
+ PUid_t intPuidThr;
+ IEta_t intEtaRestriction;
+ IEta_t intCB_CE_split;
+
+ // Class for the towers info as they should be in GCT
+ class SimpleTowerHit {
+ public:
+ IEta_t towerIeta = 0;
+ IPhi_t towerIphi = 0;
+ Et_t towerEm = 0.;
+ Et_t towerHad = 0.;
+ Et_t l1egTowerEt = 0.;
+ Et_t towerEt = 0.;
+ ap_uint<1> isBarrel = 0x1;
+ ap_uint<1> stale = 0x0;
+ ap_uint<1> stale4seed = 0x0;
+ };
+
+ // Class for the clusters info as they should arrive from HGCAL
+ class SimpleHGCluster {
+ public:
+ Pt_t pt;
+ EtaPhi_t eta;
+ EtaPhi_t phi;
+ ShowLen_t showerlength;
+ ShowLen_t coreshowerlength;
+ ShapeFeat_t spptot;
+ ShapeFeat_t szz;
+ ShapeFeat_t srrtot;
+ Meanz_t meanz;
+ PUid_t PUid;
+ ap_uint<1> stale = 0x0;
+ };
+
+ // Classes for the tower clusters
+ class SimplifiedTower {
+ public:
+ Et_t towerEm = 0.;
+ Et_t towerHad = 0.;
+ Et_t l1egTowerEt = 0.;
+
+ void fill(SimpleTowerHit Tower) {
+ towerEm = Tower.towerEm;
+ towerHad = Tower.towerHad;
+ l1egTowerEt = Tower.l1egTowerEt;
+ }
+ };
+
+ class InputTowerCluster {
+ public:
+ SimplifiedTower towerHits[45];
+ ap_uint<1> barrelSeeded = 0x0;
+ ap_uint<1> filled[45];
+
+ void fill(int idx, SimpleTowerHit Tower) {
+ towerHits[idx].fill(Tower);
+ filled[idx] = 0x1;
+ }
+
+ void init() {
+ SimplifiedTower emptyT;
+ std::fill(towerHits, towerHits + 44, emptyT);
+ std::fill(filled, filled + 44, 0x0);
+ }
+ };
+
+ class InputTowerCluster_pstn {
+ public:
+ IEta_t seedIeta = 0;
+ IPhi_t seedIphi = 0;
+
+ void fill(SimpleTowerHit Tower) {
+ seedIeta = Tower.towerIeta;
+ seedIphi = Tower.towerIphi;
+ }
+ };
+
+ // INFO : now variables are in GCT precision, they should be in NN precision
+ // after scaling, i.e. something like ap_fixed<16, 6, AP_TRN, AP_SAT>, when the
+ // full hls4ml emulation is available
+ class InputHGCluster {
+ public:
+ Pt_t pt;
+ EtaPhi_t eta;
+ ShowLen_t showerlength;
+ ShowLen_t coreshowerlength;
+ ShapeFeat_t spptot;
+ ShapeFeat_t szz;
+ ShapeFeat_t srrtot;
+ Meanz_t meanz;
+
+ void fill(SimpleHGCluster Cluster) {
+ pt = Cluster.pt;
+ eta = Cluster.eta;
+ showerlength = Cluster.showerlength;
+ coreshowerlength = Cluster.coreshowerlength;
+ spptot = Cluster.spptot;
+ szz = Cluster.szz;
+ srrtot = Cluster.srrtot;
+ meanz = Cluster.meanz;
+ }
+ };
+
+ l1t::Tau MakeTauCandidate(bool isBarrel,
+ int clNxMIdx,
+ std::vector outputsIdent,
+ std::vector outputsCalib,
+ std::vector clustersNxM_pstn);
+};
+
+/*
+████████ ██ ██ ██████ ████████ █████ ██ ██ ███ ███ ██ ███ ██ █████ ████████ ██████ ██████
+ ██ ██ ██ ██ ██ ██ ██ ██ ██ ████ ████ ██ ████ ██ ██ ██ ██ ██ ██ ██ ██
+ ██ ███████ █████ ██ ███████ ██ ██ ██ ████ ██ ██ ██ ██ ██ ███████ ██ ██ ██ ██████
+ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██
+ ██ ██ ██ ██████ ██ ██ ██ ███████ ██ ██ ██ ██ ████ ██ ██ ██ ██████ ██ ██
+*/
+
+std::unique_ptr L1NNCaloTauEmulator::initializeGlobalCache(const edm::ParameterSet& iConfig) {
+ edm::LogInfo("Initialization") << "Init NN models Global Cache " << std::endl;
+
+ std::unique_ptr GlobalCache(new NNmodels_GlobalCache);
+
+ GlobalCache->CNNmodel_CB_path = iConfig.getParameter("CNNmodel_CB_path");
+ GlobalCache->DNNident_CB_path = iConfig.getParameter("DNNident_CB_path");
+ GlobalCache->DNNcalib_CB_path = iConfig.getParameter("DNNcalib_CB_path");
+ GlobalCache->CNNmodel_CE_path = iConfig.getParameter("CNNmodel_CE_path");
+ GlobalCache->DNNident_CE_path = iConfig.getParameter("DNNident_CE_path");
+ GlobalCache->DNNcalib_CE_path = iConfig.getParameter("DNNcalib_CE_path");
+ GlobalCache->FeatScaler_CE_path = iConfig.getParameter("FeatScaler_CE_path");
+
+ // Create sessions for Tensorflow inferece
+ (GlobalCache->CNNmodel_CB) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->CNNmodel_CB_path)).fullPath());
+ (GlobalCache->CNNmodel_CBsession) = tensorflow::createSession((GlobalCache->CNNmodel_CB));
+
+ (GlobalCache->DNNident_CB) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNident_CB_path)).fullPath());
+ (GlobalCache->DNNident_CBsession) = tensorflow::createSession((GlobalCache->DNNident_CB));
+
+ (GlobalCache->DNNcalib_CB) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNcalib_CB_path)).fullPath());
+ (GlobalCache->DNNcalib_CBsession) = tensorflow::createSession((GlobalCache->DNNcalib_CB));
+
+ (GlobalCache->CNNmodel_CE) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->CNNmodel_CE_path)).fullPath());
+ (GlobalCache->CNNmodel_CEsession) = tensorflow::createSession((GlobalCache->CNNmodel_CE));
+
+ (GlobalCache->DNNident_CE) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNident_CE_path)).fullPath());
+ (GlobalCache->DNNident_CEsession) = tensorflow::createSession((GlobalCache->DNNident_CE));
+
+ (GlobalCache->DNNcalib_CE) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNcalib_CE_path)).fullPath());
+ (GlobalCache->DNNcalib_CEsession) = tensorflow::createSession((GlobalCache->DNNcalib_CE));
+
+ // Read features scaler
+ boost::property_tree::read_json(edm::FileInPath((GlobalCache->FeatScaler_CE_path)).fullPath(),
+ (GlobalCache->FeatScaler_CE));
+
+ return GlobalCache;
+}
+
+// ----Constructor and Destructor -----
+L1NNCaloTauEmulator::L1NNCaloTauEmulator(const edm::ParameterSet& iConfig, const NNmodels_GlobalCache* globalCache)
+ : l1TowersToken(consumes(iConfig.getParameter("l1CaloTowers"))),
+ hgcalTowersToken(consumes(iConfig.getParameter("hgcalTowers"))),
+
+ HGClusterToken(
+ consumes(iConfig.getParameter("HgcalClusters"))),
+ scenario(UseEmInterp::No),
+ preEmId(iConfig.getParameter("preEmId")),
+ VsPuId(iConfig.getParameter("VsPuId")),
+
+ EcalEtMinForClustering(iConfig.getParameter("EcalEtMinForClustering")),
+ HcalEtMinForClustering(iConfig.getParameter("HcalEtMinForClustering")),
+ EtMinForSeeding(iConfig.getParameter("EtMinForSeeding")),
+ EtaRestriction(iConfig.getParameter("EtaRestriction")),
+ CB_CE_split(iConfig.getParameter("CB_CE_split")),
+ PuidThr(iConfig.getParameter("PuidThr")),
+
+ IdWp90_CB(iConfig.getParameter("IdWp90_CB")),
+ IdWp95_CB(iConfig.getParameter("IdWp95_CB")),
+ IdWp99_CB(iConfig.getParameter("IdWp99_CB")),
+
+ IdWp90_CE(iConfig.getParameter("IdWp90_CE")),
+ IdWp95_CE(iConfig.getParameter("IdWp95_CE")),
+ IdWp99_CE(iConfig.getParameter("IdWp99_CE")) {
+ // Initialize HGCAL BDTs
+ if (!VsPuId.method().empty()) {
+ VsPuId.prepareTMVA();
+ }
+
+ intPuidThr = apintQuantizer(PuidThr, PUID_LSB, PUID_W);
+ intEtaRestriction = apintQuantizer(EtaRestriction, IETAPHI_LSB, IETAPHI_W);
+ intCB_CE_split = apintQuantizer(CB_CE_split, IETAPHI_LSB, IETAPHI_W) + 1;
+
+ // Create produced outputs
+ produces>("L1NNCaloTauCollectionBXV");
+
+ // Settings output
+ edm::LogInfo("Settings") << "EtaRestriction = " << EtaRestriction << " (" << intEtaRestriction << ")"
+ << " , CB_CE_split = " << CB_CE_split << "(" << intCB_CE_split
+ << ") , EtMinForSeeding = " << EtMinForSeeding
+ << " , HcalTpEtMin = " << HcalEtMinForClustering
+ << " , EcalTpEtMin = " << EcalEtMinForClustering << " , PuidThr = " << PuidThr << "("
+ << intPuidThr << ")" << std::endl;
+}
+
+void L1NNCaloTauEmulator::produce(edm::Event& iEvent, const edm::EventSetup& eSetup) {
+ // Output collection
+ std::unique_ptr> L1NNCaloTauCollectionBXV(new l1t::TauBxCollection);
+
+ // Create and Fill collection of all calotowers and their attributes
+ std::vector l1CaloTowers;
+
+ iEvent.getByToken(l1TowersToken, l1CaloTowerHandle);
+ int warnings = 0;
+ for (auto& hit : *l1CaloTowerHandle.product()) {
+ // Skip this weird towers and store warning
+ if (hit.towerIEta() == -1016 && hit.towerIPhi() == -962) {
+ warnings += 1;
+ continue;
+ }
+
+ SimpleTowerHit l1Hit;
+ l1Hit.isBarrel = 0x1;
+ l1Hit.l1egTowerEt = apfixedQuantizer(hit.l1egTowerEt(), PTET_LSB, ET_W);
+ l1Hit.towerEm = apfixedQuantizer(hit.ecalTowerEt(), PTET_LSB, ET_W);
+ l1Hit.towerHad = apfixedQuantizer(hit.hcalTowerEt(), PTET_LSB, ET_W);
+ l1Hit.towerEt = apfixedQuantizer(hit.ecalTowerEt() + hit.hcalTowerEt() + hit.l1egTowerEt(), PTET_LSB, ET_W);
+ l1Hit.towerIeta = hit.towerIEta();
+ l1Hit.towerIphi = hit.towerIPhi();
+
+ l1CaloTowers.push_back(l1Hit);
+ }
+ if (warnings != 0) {
+ edm::LogWarning("BrokenTowers") << " ** WARNING : FOUND " << warnings
+ << " TOWERS WITH towerIeta=-1016 AND towerIphi=-962" << std::endl;
+ }
+
+ iEvent.getByToken(hgcalTowersToken, hgcalTowersHandle);
+ for (auto& hit : *hgcalTowersHandle.product()) {
+ SimpleTowerHit l1Hit;
+ l1Hit.isBarrel = 0x0;
+ l1Hit.l1egTowerEt = 0.0;
+ l1Hit.towerEm = apfixedQuantizer(hit.etEm(), PTET_LSB, ET_W);
+ l1Hit.towerHad = apfixedQuantizer(hit.etHad(), PTET_LSB, ET_W);
+ l1Hit.towerEt = apfixedQuantizer(hit.etEm() + hit.etHad(), PTET_LSB, ET_W);
+ l1Hit.towerIeta = makeEndcapHwIEta(hit.eta());
+ l1Hit.towerIphi = makeEndcapHwIPhi(hit.phi());
+
+ l1CaloTowers.push_back(l1Hit);
+ }
+
+ // Sort the ECAL+HCAL+L1EGs tower sums based on total ET
+ std::sort(begin(l1CaloTowers), end(l1CaloTowers), [](const SimpleTowerHit& a, SimpleTowerHit& b) {
+ return a.towerEt > b.towerEt;
+ });
+
+ // Create and Fill the collection of 3D clusters and their attributes
+ std::vector AllHGClusters;
+ iEvent.getByToken(HGClusterToken, HGClusterHandle);
+
+ for (auto cl3dIt = HGClusterHandle->begin(0); cl3dIt != HGClusterHandle->end(0); ++cl3dIt) {
+ auto& cl3d = *cl3dIt;
+
+ // Implement cl3d PU ID as done in
+ // https://github.com/cms-sw/cmssw/blob/master/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc#L120
+ bool isEM = preEmId(*cl3dIt);
+ l1t::PFCluster cluster(cl3d.pt(), cl3d.eta(), cl3d.phi(), cl3d.hOverE());
+ if (scenario == UseEmInterp::EmOnly) // for emID objs, use EM interp as pT and set H = 0
+ {
+ if (isEM) {
+ float pt_new = cl3d.iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
+ float hoe_new = 0.;
+ cluster = l1t::PFCluster(pt_new, cl3d.eta(), cl3d.phi(), hoe_new, isEM);
+ }
+ } else if (scenario == UseEmInterp::AllKeepHad) // for all objs, replace EM part with EM interp, preserve H
+ {
+ float had_old = cl3d.pt() - cluster.emEt();
+ float em_new = cl3d.iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
+ float pt_new = had_old + em_new;
+ float hoe_new = em_new > 0 ? (had_old / em_new) : -1;
+ cluster = l1t::PFCluster(pt_new, cl3d.eta(), cl3d.phi(), hoe_new, isEM);
+ } else if (scenario == UseEmInterp::AllKeepTot) // for all objs, replace EM part with EM interp, preserve pT
+ {
+ float em_new = cl3d.iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
+ float hoe_new = em_new > 0 ? (cl3d.pt() / em_new - 1) : -1;
+ cluster = l1t::PFCluster(cl3d.pt(), cl3d.eta(), cl3d.phi(), hoe_new, isEM);
+ }
+
+ float idScore = -1.;
+ if (!VsPuId.method().empty()) {
+ idScore = VsPuId.passID(*cl3dIt, cluster);
+ idScore = cluster.egVsPUMVAOut();
+ }
+
+ float eta_hgcalCoord = correctInputEtaCl3d(cl3d.eta());
+ float meanz_hgcalCoord = correctInputMeanzCl3d(cl3d.zBarycenter());
+
+ SimpleHGCluster HGCluster;
+ HGCluster.pt = apfixedQuantizer(cl3d.pt(), PTET_LSB, PT_W);
+ HGCluster.eta = apintQuantizer(eta_hgcalCoord, ETAPHI_LSB, ETAPHI_W);
+ HGCluster.phi = apintQuantizer(cl3d.phi(), ETAPHI_LSB, ETAPHI_W);
+ HGCluster.showerlength = cl3d.showerLength();
+ HGCluster.coreshowerlength = cl3d.coreShowerLength();
+ HGCluster.spptot = apintQuantizer(cl3d.sigmaPhiPhiTot(), SHAPEFEAT_LSB, SHAPEFEAT_W);
+ HGCluster.szz = apintQuantizer(cl3d.sigmaZZ(), SZZ_LSB, SHAPEFEAT_W);
+ HGCluster.srrtot = apintQuantizer(cl3d.sigmaRRTot(), SHAPEFEAT_LSB, SHAPEFEAT_W);
+ HGCluster.meanz = apintQuantizer(meanz_hgcalCoord, MEANZ_LSB, MEANZ_W);
+ HGCluster.PUid = apintQuantizer(idScore, PUID_LSB, PUID_W);
+
+ AllHGClusters.push_back(HGCluster);
+ }
+
+ // Order the collection in pt (the input to the GCT will be pt ordered)
+ std::sort(begin(AllHGClusters), end(AllHGClusters), [](const SimpleHGCluster& a, SimpleHGCluster& b) {
+ return a.pt > b.pt;
+ });
+
+ /*
+ // END OF SOFTWARE PRECISION SECTION
+ // up to here treated inputs from simulation with SW precision
+ // to massage them into the HW precision varibales as they are
+ // forseen (roughly) to be available at the GCT Sum card level
+ // ------------------------------------------------------------- */
+
+ // Make NxM TowerClusters and HGClusters collections for TauMinator
+ std::vector l1TowerClustersNxM_CB;
+ std::vector l1TowerClustersNxM_CB_pstn;
+ std::vector l1TowerClustersNxM_CE;
+ std::vector l1TowerClustersNxM_CE_pstn;
+ std::vector HGClusters;
+
+ // Supporting collection of endcap clusters before cl3d matching
+ std::vector AllL1TowerClustersNxM_CE;
+ std::vector AllL1TowerClustersNxM_CE_pstn;
+
+ int Nclusters_CB = 0;
+ int AllNclusters_CE = 0;
+ bool caloTauSeedingFinished = false;
+ // Loop for seeding of clNxM objects
+ while (!caloTauSeedingFinished) {
+ InputTowerCluster clNxM;
+ clNxM.init();
+ InputTowerCluster_pstn clNxM_pstn;
+ bool seeded = false;
+
+ for (auto& l1CaloTower : l1CaloTowers) {
+ // Skip seeding in towers that would make the cluster extend in HF
+ // Skip l1CaloTowers which are already used by this clusters' mask
+ if (ap_abs(l1CaloTower.towerIeta) > Eta_limit || ap_abs(l1CaloTower.towerIeta) > intEtaRestriction ||
+ l1CaloTower.stale4seed) {
+ continue;
+ }
+
+ // If not seded do the seeding
+ if (!seeded) {
+ // The leading unused tower has ET < min, stop jet clustering
+ if (l1CaloTower.towerEt < EtMinForSeeding) {
+ caloTauSeedingFinished = true;
+ continue;
+ }
+
+ clNxM.fill(seedIdx, l1CaloTower);
+ clNxM_pstn.fill(l1CaloTower);
+ if (l1CaloTower.isBarrel) {
+ clNxM.barrelSeeded = 0x1;
+ }
+
+ l1CaloTower.stale4seed = 0x1;
+ l1CaloTower.stale = 0x1;
+ seeded = true;
+
+ continue;
+ }
+
+ dIEtaPhi_t d_iEta = tower_dIEta(l1CaloTower.towerIeta, clNxM_pstn.seedIeta);
+ dIEtaPhi_t d_iPhi = dPhi(l1CaloTower.towerIphi, clNxM_pstn.seedIphi);
+
+ // Stale tower for seeding if it would lead to overalp between clusters
+ if ((ap_abs(d_iEta) <= IEta_dim - 1 && ap_abs(d_iPhi) <= IPhi_dim - 1)) {
+ l1CaloTower.stale4seed = 0x1;
+ }
+
+ } // End for loop over TPs
+
+ // Pushback seeds split in barrel and endcap
+ if (seeded) {
+ if (ap_abs(clNxM_pstn.seedIeta) <= intCB_CE_split) {
+ l1TowerClustersNxM_CB.push_back(clNxM);
+ l1TowerClustersNxM_CB_pstn.push_back(clNxM_pstn);
+ Nclusters_CB++;
+ } else {
+ AllL1TowerClustersNxM_CE.push_back(clNxM);
+ AllL1TowerClustersNxM_CE_pstn.push_back(clNxM_pstn);
+ AllNclusters_CE++;
+ }
+ }
+
+ } // End while loop of TowerClusters seeding
+
+ // Loop for barrel NxM TowerClusters clustering starting from the seeds
+ for (int clNxMIdx = 0; clNxMIdx < Nclusters_CB; clNxMIdx++) {
+ for (auto& l1CaloTower : l1CaloTowers) {
+ // Skip l1CaloTowers which are already used
+ if (l1CaloTower.stale) {
+ continue;
+ }
+
+ dIEtaPhi_t d_iEta = tower_dIEta(l1CaloTower.towerIeta, l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIeta);
+ dIEtaPhi_t d_iPhi =
+ dPhi(l1CaloTower.towerIphi, l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIphi);
+ int hitIdx = d_iEta * 9 + d_iPhi + seedIdx;
+
+ // Cluster all towers in a NxM towers mask
+ if ((ap_abs(d_iEta) <= (IEta_dim - 1) / 2 && ap_abs(d_iPhi) <= (IPhi_dim - 1) / 2)) {
+ l1TowerClustersNxM_CB[clNxMIdx].fill(hitIdx, l1CaloTower);
+ l1CaloTower.stale = 0x1;
+ }
+
+ } // End for loop over TPs
+
+ } // End while loop of barrel TowerClusters creation
+
+ // In the endcap cross-loop over clNxM and cl3d to match them
+ // (we can do it before full clustering just using the seed info)
+ int Nclusters_CE = 0;
+ for (int clNxMIdx = 0; clNxMIdx < AllNclusters_CE; clNxMIdx++) {
+ bool matched = false;
+ for (auto& HGCluster : AllHGClusters) {
+ // In case the clNxM or HGCluster have already been matched just continue through the list to the end
+ // only use cl3ds above 4GeV and above -0.10 pu id
+ if (matched || HGCluster.stale || HGCluster.pt < Pt_t(4.) || HGCluster.PUid < intPuidThr) {
+ continue;
+ }
+
+ dEtaPhi_t d_iEta = tw2cl_dEta(HGCluster.eta, AllL1TowerClustersNxM_CE_pstn[clNxMIdx].seedIeta);
+ dEtaPhi_t d_iPhi = tw2cl_dPhi(HGCluster.phi, AllL1TowerClustersNxM_CE_pstn[clNxMIdx].seedIphi);
+ matched = d_iEta * d_iEta + d_iPhi * d_iPhi < R2cone;
+
+ if (matched) {
+ HGCluster.stale = 0x1;
+ InputHGCluster cl3d;
+ cl3d.fill(HGCluster);
+ HGClusters.push_back(cl3d);
+ l1TowerClustersNxM_CE.push_back(AllL1TowerClustersNxM_CE[clNxMIdx]);
+ l1TowerClustersNxM_CE_pstn.push_back(AllL1TowerClustersNxM_CE_pstn[clNxMIdx]);
+ Nclusters_CE++;
+ }
+
+ } // End for loop over cl3ds
+
+ } // End for loop over clNxM
+
+ // Loop for endcap matched NxM TowerClusters clustering starting from the seeds just found
+ for (int clNxMIdx = 0; clNxMIdx < Nclusters_CE; clNxMIdx++) {
+ for (auto& l1CaloTower : l1CaloTowers) {
+ // Skip l1CaloTowers which are already used
+ if (l1CaloTower.stale) {
+ continue;
+ }
+
+ dIEtaPhi_t d_iEta = tower_dIEta(l1CaloTower.towerIeta, l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIeta);
+ dIEtaPhi_t d_iPhi =
+ dPhi(l1CaloTower.towerIphi, l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIphi);
+ int hitIdx = d_iEta * 9 + d_iPhi + seedIdx;
+
+ // Cluster all towers in a NxM towers mask
+ if ((ap_abs(d_iEta) <= (IEta_dim - 1) / 2 && ap_abs(d_iPhi) <= (IPhi_dim - 1) / 2)) {
+ l1TowerClustersNxM_CE[clNxMIdx].fill(hitIdx, l1CaloTower);
+ l1CaloTower.stale = 0x1;
+ }
+
+ } // End for loop over TPs
+
+ } // End while loop of barrel TowerClusters creation
+
+ // Barrel TauMinator application
+ tensorflow::setLogging("2");
+ int batchSize_CB = (int)(Nclusters_CB);
+ tensorflow::TensorShape imageShape_CB({batchSize_CB, IEta_dim, IPhi_dim, 2});
+ tensorflow::TensorShape positionShape_CB({batchSize_CB, 2});
+ tensorflow::Tensor TowerClusterImage_CB(tensorflow::DT_FLOAT, imageShape_CB);
+ tensorflow::Tensor TowerClusterPosition_CB(tensorflow::DT_FLOAT, positionShape_CB);
+
+ for (int clNxMIdx = 0; clNxMIdx < Nclusters_CB; clNxMIdx++) {
+ // Fill inputs for Tensorflow inference
+ for (int eta = 0; eta < IEta_dim; ++eta) {
+ for (int phi = 0; phi < IPhi_dim; ++phi) {
+ int towerIdx = eta * IPhi_dim + phi;
+ TowerClusterImage_CB.tensor()(clNxMIdx, eta, phi, 0) =
+ (l1TowerClustersNxM_CB[clNxMIdx].towerHits[towerIdx].l1egTowerEt.to_float() +
+ l1TowerClustersNxM_CB[clNxMIdx].towerHits[towerIdx].towerEm.to_float());
+ TowerClusterImage_CB.tensor()(clNxMIdx, eta, phi, 1) =
+ (l1TowerClustersNxM_CB[clNxMIdx].towerHits[towerIdx].towerHad.to_float());
+ }
+ }
+
+ TowerClusterPosition_CB.tensor()(clNxMIdx, 0) = floatIEta(l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIeta);
+ TowerClusterPosition_CB.tensor()(clNxMIdx, 1) = floatIPhi(l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIphi);
+ }
+
+ if (batchSize_CB >
+ 0) // from CMSSW_14_0_X tensorflow does not seem to be able to deal with a tensor of dimension 0 anymore
+ {
+ // Apply CNN model
+ tensorflow::NamedTensorList CNNmodel_CBinputList = {{"TowerClusterImage", TowerClusterImage_CB},
+ {"TowerClusterPosition", TowerClusterPosition_CB}};
+ std::vector CNNmodel_CBoutputs;
+ tensorflow::run((globalCache()->CNNmodel_CBsession),
+ CNNmodel_CBinputList,
+ {"TauMinator_CB_conv/middleMan/concat"},
+ &CNNmodel_CBoutputs);
+ tensorflow::NamedTensorList DNN_CBinputsList = {{"middleMan", CNNmodel_CBoutputs[0]}};
+
+ // Apply DNN for identification
+ std::vector DNN_CBoutputsIdent;
+ tensorflow::run((globalCache()->DNNident_CBsession),
+ DNN_CBinputsList,
+ {"TauMinator_CB_ident/sigmoid_IDout/Sigmoid"},
+ &DNN_CBoutputsIdent);
+
+ // Apply DNN for calibration
+ std::vector DNN_CBoutputsCalib;
+ tensorflow::run((globalCache()->DNNcalib_CBsession),
+ DNN_CBinputsList,
+ {"TauMinator_CB_calib/DNNout/MatMul"},
+ &DNN_CBoutputsCalib);
+
+ // Fill the output collection of L1 taus with the barrel candidates
+ for (int clNxMIdx = 0; clNxMIdx < Nclusters_CB; clNxMIdx++) {
+ l1t::Tau l1Tau =
+ MakeTauCandidate(true, clNxMIdx, DNN_CBoutputsIdent, DNN_CBoutputsCalib, l1TowerClustersNxM_CB_pstn);
+ if (l1Tau.pt() < 0) {
+ continue;
+ }
+ L1NNCaloTauCollectionBXV->push_back(0, l1Tau);
+ }
+ }
+
+ // Endcap TauMinator application
+ int batchSize_CE = (int)(Nclusters_CE);
+ tensorflow::TensorShape imageShape_CE({batchSize_CE, IEta_dim, IPhi_dim, 2});
+ tensorflow::TensorShape positionShape_CE({batchSize_CE, 2});
+ tensorflow::TensorShape cl3dfeatShape_CE({batchSize_CE, 8});
+ tensorflow::Tensor TowerClusterImage_CE(tensorflow::DT_FLOAT, imageShape_CE);
+ tensorflow::Tensor TowerClusterPosition_CE(tensorflow::DT_FLOAT, positionShape_CE);
+ tensorflow::Tensor Cl3dShapeFeatures_CE(tensorflow::DT_FLOAT, cl3dfeatShape_CE);
+
+ for (int clNxMIdx = 0; clNxMIdx < Nclusters_CE; clNxMIdx++) {
+ // Indexing of cl3ds is the same as the one of clNxMs
+ InputHGCluster HGClu = HGClusters[clNxMIdx];
+
+ // Fill inputs for Tensorflow inference
+ for (int eta = 0; eta < IEta_dim; ++eta) {
+ for (int phi = 0; phi < IPhi_dim; ++phi) {
+ int towerIdx = eta * IPhi_dim + phi;
+ TowerClusterImage_CE.tensor()(clNxMIdx, eta, phi, 0) =
+ (l1TowerClustersNxM_CE[clNxMIdx].towerHits[towerIdx].l1egTowerEt.to_float() +
+ l1TowerClustersNxM_CE[clNxMIdx].towerHits[towerIdx].towerEm.to_float());
+ TowerClusterImage_CE.tensor()(clNxMIdx, eta, phi, 1) =
+ (l1TowerClustersNxM_CE[clNxMIdx].towerHits[towerIdx].towerHad.to_float());
+ }
+ }
+
+ TowerClusterPosition_CE.tensor()(clNxMIdx, 0) = floatIEta(l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIeta);
+ TowerClusterPosition_CE.tensor()(clNxMIdx, 1) = floatIPhi(l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIphi);
+
+ Cl3dShapeFeatures_CE.tensor()(clNxMIdx, 0) = inputScaler(HGClu.pt.to_float(), "pt");
+ Cl3dShapeFeatures_CE.tensor()(clNxMIdx, 1) = inputScaler(abs(floatEta(HGClu.eta)), "eta");
+ Cl3dShapeFeatures_CE.tensor()(clNxMIdx, 2) = inputScaler(HGClu.showerlength.to_float(), "showerlength");
+ Cl3dShapeFeatures_CE.tensor()(clNxMIdx, 3) =
+ inputScaler(HGClu.coreshowerlength.to_float(), "coreshowerlength");
+ Cl3dShapeFeatures_CE.tensor()(clNxMIdx, 4) = inputScaler(floatShape(HGClu.spptot), "spptot");
+ Cl3dShapeFeatures_CE.tensor()(clNxMIdx, 5) = inputScaler(floatSzz(HGClu.szz), "szz");
+ Cl3dShapeFeatures_CE.tensor()(clNxMIdx, 6) = inputScaler(floatShape(HGClu.srrtot), "srrtot");
+ Cl3dShapeFeatures_CE.tensor()(clNxMIdx, 7) = inputScaler(floatMeanZHgcalCoord(HGClu.meanz), "meanz");
+ }
+
+ if (batchSize_CE >
+ 0) // from CMSSW_14_0_X tensorflow does not seem to be able to deal with a tensor of dimension 0 anymore
+ {
+ // Apply CNN model
+ tensorflow::NamedTensorList CNNmodel_CEinputList = {{"TowerClusterImage", TowerClusterImage_CE},
+ {"TowerClusterPosition", TowerClusterPosition_CE},
+ {"AssociatedCl3dFeatures", Cl3dShapeFeatures_CE}};
+ std::vector CNNmodel_CEoutputs;
+ tensorflow::run((globalCache()->CNNmodel_CEsession),
+ CNNmodel_CEinputList,
+ {"TauMinator_CE_conv/middleMan/concat"},
+ &CNNmodel_CEoutputs);
+ tensorflow::NamedTensorList DNN_CEinputsList = {{"middleMan", CNNmodel_CEoutputs[0]}};
+
+ // Apply DNN for identification
+ std::vector DNN_CEoutputsIdent;
+ tensorflow::run((globalCache()->DNNident_CEsession),
+ DNN_CEinputsList,
+ {"TauMinator_CE_ident/sigmoid_IDout/Sigmoid"},
+ &DNN_CEoutputsIdent);
+
+ // Apply DNN for calibration
+ std::vector DNN_CEoutputsCalib;
+ tensorflow::run((globalCache()->DNNcalib_CEsession),
+ DNN_CEinputsList,
+ {"TauMinator_CE_calib/LIN_DNNout/Relu"},
+ &DNN_CEoutputsCalib);
+
+ // Fill the output collection of L1 taus with the endcap candidates
+ for (int clNxMIdx = 0; clNxMIdx < Nclusters_CE; clNxMIdx++) {
+ l1t::Tau l1Tau =
+ MakeTauCandidate(false, clNxMIdx, DNN_CEoutputsIdent, DNN_CEoutputsCalib, l1TowerClustersNxM_CE_pstn);
+ if (l1Tau.pt() < 0) {
+ continue;
+ }
+ L1NNCaloTauCollectionBXV->push_back(0, l1Tau);
+ }
+ }
+
+ // Fill output
+ iEvent.put(std::move(L1NNCaloTauCollectionBXV), "L1NNCaloTauCollectionBXV");
+
+} // End of produce function
+
+template
+outPrecision L1NNCaloTauEmulator::dPhi(inPrecision iPhi_1, inPrecision iPhi_2) {
+ outPrecision dphi = iPhi_1 - iPhi_2;
+
+ outPrecision dphi0 = dphi > outPrecision(INTPHI_PI) ? outPrecision(dphi - INTPHI_2PI) : dphi;
+ outPrecision dphi1 = dphi <= outPrecision(-INTPHI_PI) ? outPrecision(dphi + INTPHI_2PI) : dphi;
+
+ outPrecision result = dphi > outPrecision(0) ? dphi0 : dphi1;
+
+ return result;
+}
+
+L1NNCaloTauEmulator::dIEtaPhi_t L1NNCaloTauEmulator::tower_dIEta(IEta_t iEta_1, IEta_t iEta_2) {
+ ap_int<12> mult = iEta_1 * iEta_2;
+ dIEtaPhi_t result = iEta_1 - iEta_2;
+ if (mult < 0) {
+ result = iEta_1 > 0 ? result - 1 : result + 1;
+ }
+
+ return result;
+}
+
+L1NNCaloTauEmulator::dEtaPhi_t L1NNCaloTauEmulator::tw2cl_dPhi(EtaPhi_t iPhi_1, IPhi_t iPhi_2) {
+ EtaPhi_t shiftediPhi_2 = iPhi_2 <= IPhi_t(36) ? EtaPhi_t(iPhi_2) : EtaPhi_t(iPhi_2 - INTPHI_2PI + 1);
+
+ EtaPhi_t fineiPhi_2 = shiftediPhi_2 * (IETAPHI_LSB / ETAPHI_LSB);
+ // subrtaction of half rescaling corrects from edge to center of tower
+ fineiPhi_2 = fineiPhi_2 > EtaPhi_t(0) ? EtaPhi_t(fineiPhi_2 - (IETAPHI_LSB / ETAPHI_LSB) / 2)
+ : EtaPhi_t(fineiPhi_2 + (IETAPHI_LSB / ETAPHI_LSB) / 2);
+
+ return dPhi(iPhi_1, fineiPhi_2);
+}
+
+L1NNCaloTauEmulator::dEtaPhi_t L1NNCaloTauEmulator::tw2cl_dEta(EtaPhi_t iEta_1, IEta_t iEta_2) {
+ // change from hgcal frame to barrel-centered frame
+ EtaPhi_t framechangeCl3d = 303; // 303*pi/720 = 1.322
+ iEta_1 = iEta_1 > EtaPhi_t(0) ? EtaPhi_t(iEta_1 + framechangeCl3d) : EtaPhi_t(iEta_1 - framechangeCl3d);
+
+ // the actual depth is 330 but 329 corrects for 0.0808 tower
+ EtaPhi_t barrelEtaDepth = 329;
+ EtaPhi_t fineiEta_2 = barrelEtaDepth + (iEta_2 - IETAHGCAL_OFFSET) * (IETAHGCAL_LSB / ETAPHI_LSB);
+
+ return iEta_1 - fineiEta_2;
+}
+
+L1NNCaloTauEmulator::IEta_t L1NNCaloTauEmulator::makeEndcapHwIEta(float eta) {
+ IEta_t ieta = floor(eta / IETAHGCAL_LSB);
+ // +1 because flooring gets it 1 unit lower when negative
+ ieta = ieta < IEta_t(0) ? IEta_t(ieta + 1) : ieta;
+
+ return ieta;
+}
+
+L1NNCaloTauEmulator::IPhi_t L1NNCaloTauEmulator::makeEndcapHwIPhi(float phi) {
+ phi = phi < 0 ? phi + 2 * M_PI : phi;
+
+ // +1 because tower 0 does not exist
+ return floor(phi / IETAPHI_LSB) + 1;
+}
+
+template
+ap_int L1NNCaloTauEmulator::ap_abs(ap_int x) {
+ ap_int result;
+ if (x < 0) {
+ result = -x;
+ } else {
+ result = x;
+ }
+
+ return result;
+}
+
+template
+ap_ufixed L1NNCaloTauEmulator::ap_abs(ap_fixed x) {
+ ap_ufixed result;
+ if (x < 0) {
+ result = -x;
+ } else {
+ result = x;
+ }
+
+ return result;
+}
+
+float L1NNCaloTauEmulator::apfixedQuantizer(float inputF, float LSB, int nbits) {
+ return min(floor(inputF / LSB), float(pow(2, nbits) - 1)) * LSB;
+}
+
+int L1NNCaloTauEmulator::apintQuantizer(float inputF, float LSB, int nbits) {
+ return min(floor(inputF / LSB), float(pow(2, nbits) - 1));
+}
+
+float L1NNCaloTauEmulator::inputScaler(float inputF, std::string feature) {
+ float mean = (globalCache()->FeatScaler_CE).get_child(feature).get("mean");
+ float std = (globalCache()->FeatScaler_CE).get_child(feature).get("std");
+
+ return (inputF - mean) / std;
+}
+
+float L1NNCaloTauEmulator::correctInputEtaCl3d(float eta) {
+ return eta > 0 ? eta - ETAHGCAL_OFFSET : eta + ETAHGCAL_OFFSET;
+}
+
+float L1NNCaloTauEmulator::correctInputMeanzCl3d(float meanz) { return CM2MM * (abs(meanz) - MEANZ_OFFSET); }
+
+float L1NNCaloTauEmulator::floatIEta(IEta_t eta) {
+ // transform eta of towers from integer to float, correcting for different barrel/endcap LSB
+ float feta;
+ if (abs(eta) > IETAHGCAL_OFFSET) {
+ if (eta > 0) {
+ feta = IETAHGCAL_OFFSET * IETAPHI_LSB - (IETAHGCAL_LSB - IETAHGCAL_LSBp) +
+ (eta.to_float() - IETAHGCAL_OFFSET) * IETAHGCAL_LSB;
+ } else {
+ feta = -IETAHGCAL_OFFSET * IETAPHI_LSB + (IETAHGCAL_LSB - IETAHGCAL_LSBp) +
+ (eta.to_float() + IETAHGCAL_OFFSET) * IETAHGCAL_LSB;
+ }
+ } else {
+ feta = eta.to_float() * IETAPHI_LSB;
+ }
+
+ // shift by half a tower to consider the tower center instead of the edge
+ return feta > 0 ? feta - IETAPHI_LSB / 2 : feta + IETAPHI_LSB / 2;
+}
+
+float L1NNCaloTauEmulator::floatIPhi(IPhi_t phi) {
+ float fphi = phi.to_float();
+ // add 2pi + 1 because tower 0 does not exist
+ fphi = fphi > INTPHI_PI ? fphi - INTPHI_2PI + 1 : fphi;
+ fphi *= IETAPHI_LSB;
+
+ // shift by half a tower to consider the tower center instead of the edge
+ return fphi > 0 ? fphi - IETAPHI_LSB / 2 : fphi + IETAPHI_LSB / 2;
+}
+
+l1t::Tau L1NNCaloTauEmulator::MakeTauCandidate(
+ bool isBarrel,
+ int clNxMIdx,
+ std::vector outputsIdent,
+ std::vector outputsCalib,
+ std::vector clustersNxM_pstn) {
+ int seedIeta = clustersNxM_pstn[clNxMIdx].seedIeta;
+ int seedIphi = clustersNxM_pstn[clNxMIdx].seedIphi;
+
+ if (seedIeta > intEtaRestriction) {
+ return l1t::Tau(reco::Candidate::PolarLorentzVector(-1, 0, 0, 0), -1, 0, 0, 0, 0);
+ ;
+ }
+
+ float tau_IDscore = outputsIdent[0].matrix()(0, clNxMIdx);
+ float tau_calibPt = outputsCalib[0].matrix()(0, clNxMIdx);
+ float tau_eta = floatIEta(seedIeta);
+ float tau_phi = floatIPhi(seedIphi);
+
+ // Assign increasing quality to higher scoring candidates
+ int quality = 0;
+ if (isBarrel) {
+ // 99% WP
+ if (tau_IDscore > IdWp99_CB) {
+ quality = 1;
+ }
+ // 95% WP
+ if (tau_IDscore > IdWp95_CB) {
+ quality = 2;
+ }
+ // 90% WP
+ if (tau_IDscore > IdWp90_CB) {
+ quality = 3;
+ }
+ } else {
+ // 99% WP
+ if (tau_IDscore > IdWp99_CE) {
+ quality = 1;
+ }
+ // 95% WP
+ if (tau_IDscore > IdWp95_CE) {
+ quality = 2;
+ }
+ // 90% WP
+ if (tau_IDscore > IdWp90_CE) {
+ quality = 3;
+ }
+ }
+
+ reco::Candidate::PolarLorentzVector tauP4 = reco::Candidate::PolarLorentzVector(tau_calibPt, tau_eta, tau_phi, 0);
+
+ // store ID score multiplied by 10E4 to have good precision even using the Phase1 tau int iso format
+ // (this is stored just in case for possible additional offline studies)
+ // tau initialisation = (p4, pt, eta, phi, qual, iso)
+ l1t::Tau l1Tau = l1t::Tau(tauP4, tau_calibPt, tau_eta, tau_phi, quality, tau_IDscore * 10E4);
+ l1Tau.setTowerIEta(seedIeta);
+ l1Tau.setTowerIPhi(seedIphi);
+
+ return l1Tau;
+}
+
+void L1NNCaloTauEmulator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+ edm::ParameterSetDescription desc;
+
+ desc.add("l1CaloTowers", edm::InputTag("l1tEGammaClusterEmuProducer", "L1CaloTowerCollection"));
+ desc.add("hgcalTowers", edm::InputTag("l1tHGCalTowerProducer", "HGCalTowerProcessor"));
+ desc.add("HgcalClusters",
+ edm::InputTag("l1tHGCalBackEndLayer2Producer", "HGCalBackendLayer2Processor3DClustering"));
+
+ desc.add("preEmId", "hOverE < 0.3 && hOverE >= 0");
+ {
+ edm::ParameterSetDescription psd0;
+ psd0.add("isPUFilter", true);
+ psd0.add("preselection", "");
+ psd0.add("method", "BDT");
+ {
+ edm::ParameterSetDescription vpsd2;
+ vpsd2.add("name");
+ vpsd2.add("value");
+ std::vector temp2;
+ temp2.reserve(5);
+ {
+ edm::ParameterSet temp3;
+ temp3.addParameter("name", "eMax");
+ temp3.addParameter("value", "eMax()");
+ temp2.push_back(temp3);
+ }
+ {
+ edm::ParameterSet temp3;
+ temp3.addParameter("name", "eMaxOverE");
+ temp3.addParameter("value", "eMax()/energy()");
+ temp2.push_back(temp3);
+ }
+ {
+ edm::ParameterSet temp3;
+ temp3.addParameter("name", "sigmaPhiPhiTot");
+ temp3.addParameter("value", "sigmaPhiPhiTot()");
+ temp2.push_back(temp3);
+ }
+ {
+ edm::ParameterSet temp3;
+ temp3.addParameter("name", "sigmaRRTot");
+ temp3.addParameter("value", "sigmaRRTot()");
+ temp2.push_back(temp3);
+ }
+ {
+ edm::ParameterSet temp3;
+ temp3.addParameter("name", "triggerCells90percent");
+ temp3.addParameter("value", "triggerCells90percent()");
+ temp2.push_back(temp3);
+ }
+ psd0.addVPSet("variables", vpsd2, temp2);
+ }
+ psd0.add(
+ "weightsFile", "L1Trigger/Phase2L1ParticleFlow/data/hgcal_egID/Photon_Pion_vs_Neutrino_BDTweights_1116.xml.gz");
+ psd0.add("wp", "-0.10");
+ desc.add("VsPuId", psd0);
+ }
+
+ desc.add("EcalEtMinForClustering", 0.0);
+ desc.add("HcalEtMinForClustering", 0.0);
+ desc.add("EtMinForSeeding", 2.5);
+ desc.add("EtaRestriction", 2.4);
+ desc.add("CB_CE_split", 1.55);
+ desc.add("PuidThr", -0.1);
+
+ desc.add("CNNmodel_CB_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/CNNmodel_CB.pb");
+ desc.add("DNNident_CB_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNident_CB.pb");
+ desc.add("DNNcalib_CB_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNcalib_CB.pb");
+ desc.add("CNNmodel_CE_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/CNNmodel_CE.pb");
+ desc.add("DNNident_CE_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNident_CE.pb");
+ desc.add("DNNcalib_CE_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNcalib_CE.pb");
+ desc.add("FeatScaler_CE_path", "L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/Cl3dFeatScaler_CE.json");
+
+ desc.add("IdWp90_CB", 0.706);
+ desc.add("IdWp95_CB", 0.3432);
+ desc.add("IdWp99_CB", 0.0337);
+ desc.add("IdWp90_CE", 0.5711);
+ desc.add("IdWp95_CE", 0.2742);
+ desc.add("IdWp99_CE", 0.0394);
+
+ desc.add("DEBUG", false);
+
+ descriptions.add("l1tNNCaloTauEmulator", desc);
+}
+
+DEFINE_FWK_MODULE(L1NNCaloTauEmulator);
diff --git a/L1Trigger/L1CaloTrigger/plugins/L1NNCaloTauProducer.cc b/L1Trigger/L1CaloTrigger/plugins/L1NNCaloTauProducer.cc
new file mode 100644
index 0000000000000..2a17bdd84915d
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/plugins/L1NNCaloTauProducer.cc
@@ -0,0 +1,933 @@
+/* -*- C++ -*-
+
+Package: L1CaloTrigger
+Class: L1NNCaloTauProducer
+Frinedly name: The TauMinator
+
+\class L1NNCaloTauProducer L1NNCaloTauProducer.cc
+
+Description:
+Perform reconstruction and identification of tau
+candidates at L1 Trigger with a CNN.
+
+Implementation:
+Take as input the HCAL TPs, the ECAL TPs from
+l1tEGammaClusterEmuProducer, and the HGCAL TPs
+from l1tHGCalTowerProducer and l1tHGCalBackEndLayer2Producer.
+Proceed to clustering of trigger towers in NxM
+clusters, match to HGcal 3D clusters in the endcap.
+Finally apply the CNNs.
+
+Original Author: Jona Motta
+Created: Tue May 30th 2023
+
+*/
+
+#include
+#include
+#include
+
+#include "boost/property_tree/ptree.hpp"
+#include "boost/property_tree/json_parser.hpp"
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h"
+#include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
+#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
+#include "DataFormats/L1TParticleFlow/interface/PFCluster.h"
+#include "DataFormats/L1THGCal/interface/HGCalTower.h"
+#include "DataFormats/Math/interface/deltaPhi.h"
+#include "DataFormats/L1Trigger/interface/Tau.h"
+
+#include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h"
+#include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h"
+
+#include "L1Trigger/L1THGCal/interface/backend/HGCalTriggerClusterIdentificationBase.h"
+#include "L1Trigger/Phase2L1ParticleFlow/interface/HGC3DClusterEgID.h"
+#include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
+
+#include "PhysicsTools/TensorFlow/interface/TensorFlow.h"
+
+struct NNmodels_GlobalCache {
+ std::string CNNmodel_CB_path;
+ std::string DNNident_CB_path;
+ std::string DNNcalib_CB_path;
+
+ std::string CNNmodel_CE_path;
+ std::string DNNident_CE_path;
+ std::string DNNcalib_CE_path;
+ std::string FeatScaler_CE_path;
+ boost::property_tree::ptree FeatScaler_CE;
+
+ tensorflow::GraphDef* CNNmodel_CB;
+ tensorflow::GraphDef* DNNident_CB;
+ tensorflow::GraphDef* DNNcalib_CB;
+
+ tensorflow::Session* CNNmodel_CBsession;
+ tensorflow::Session* DNNident_CBsession;
+ tensorflow::Session* DNNcalib_CBsession;
+
+ tensorflow::GraphDef* CNNmodel_CE;
+ tensorflow::GraphDef* DNNident_CE;
+ tensorflow::GraphDef* DNNcalib_CE;
+
+ tensorflow::Session* CNNmodel_CEsession;
+ tensorflow::Session* DNNident_CEsession;
+ tensorflow::Session* DNNcalib_CEsession;
+};
+
+class L1NNCaloTauProducer : public edm::stream::EDProducer> {
+public:
+ explicit L1NNCaloTauProducer(const edm::ParameterSet&, const NNmodels_GlobalCache*);
+
+ static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+ static std::unique_ptr initializeGlobalCache(const edm::ParameterSet&);
+ static void globalEndJob(const NNmodels_GlobalCache*){/*do nothing*/};
+
+private:
+ //----edm control---
+ void produce(edm::Event&, const edm::EventSetup&) override;
+
+ //----private functions----
+ int tower_dIPhi(int& iPhi_1, int& iPhi_2) const;
+ int tower_dIEta(int& iEta_1, int& iEta_2) const;
+ int endcap_iphi(float& phi) const;
+ int endcap_ieta(float& eta) const;
+ float inputQuantizer(float inputF, float LSB, int nbits);
+ float inputScaler(float inputF, std::string feature);
+
+ //----tokens and handles----
+ edm::EDGetTokenT l1TowersToken;
+ edm::Handle l1CaloTowerHandle;
+
+ edm::EDGetToken hgcalTowersToken;
+ edm::Handle hgcalTowersHandle;
+
+ edm::EDGetTokenT HGClusterToken;
+ edm::Handle HGClusterHandle;
+
+ //----private variables----
+ enum class UseEmInterp { No, EmOnly, AllKeepHad, AllKeepTot };
+ UseEmInterp scenario;
+ StringCutObjectSelector preEmId;
+ l1tpf::HGC3DClusterEgID VsPuId;
+
+ double EcalEtMinForClustering;
+ double HcalEtMinForClustering;
+ double EtMinForSeeding;
+ double EtaRestriction;
+ double CB_CE_split;
+
+ double IdWp90_CB;
+ double IdWp95_CB;
+ double IdWp99_CB;
+
+ double IdWp90_CE;
+ double IdWp95_CE;
+ double IdWp99_CE;
+
+ // hardoced dimensions of the tower clusters
+ const int seedIdx = 22;
+ const int IEta_dim = 5;
+ const int IPhi_dim = 9;
+ const float Eta_dim = 0.2;
+ const float Phi_dim = 0.4;
+ const float Eta_dim_seed = 0.35;
+ const float Phi_dim_seed = 0.7;
+ const float Eta_limit = 2.83;
+
+ // classes of objects used only in this producer
+ class SimpleTowerHit {
+ public:
+ float towerEta = -99.;
+ float towerPhi = -99.;
+ float towerEm = 0.;
+ float towerHad = 0.;
+ float l1egTowerEt = 0.;
+ float towerEt = 0.;
+ int towerIeta = -99;
+ int towerIphi = -99;
+ bool isBarrel = true;
+ bool stale = false;
+ bool stale4seed = false;
+ };
+
+ class SimpleTowerCluster {
+ public:
+ bool barrelSeeded = false;
+ int seedIeta = -99;
+ int seedIphi = -99;
+ float seedEta = -99.;
+ float seedPhi = -99.;
+ float rawEt = 0.;
+ float IDscore = -99.;
+ float calibPt = -99.;
+
+ std::vector towerHits;
+
+ void InitHits(int N, int M) { towerHits.resize(N * M); }
+ };
+
+ class SimpleHGCluster {
+ public:
+ float pt = -99.;
+ float eta = -99.;
+ float phi = -99.;
+ float showerlength = -99.;
+ float coreshowerlength = -99.;
+ float spptot = -99.;
+ float szz = -99.;
+ float srrtot = -99.;
+ float meanz = -99.;
+ bool stale = false;
+ };
+};
+
+/*
+████████ ██ ██ ██████ ████████ █████ ██ ██ ███ ███ ██ ███ ██ █████ ████████ ██████ ██████
+ ██ ██ ██ ██ ██ ██ ██ ██ ██ ████ ████ ██ ████ ██ ██ ██ ██ ██ ██ ██ ██
+ ██ ███████ █████ ██ ███████ ██ ██ ██ ████ ██ ██ ██ ██ ██ ███████ ██ ██ ██ ██████
+ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██
+ ██ ██ ██ ██████ ██ ██ ██ ███████ ██ ██ ██ ██ ████ ██ ██ ██ ██████ ██ ██
+*/
+
+std::unique_ptr L1NNCaloTauProducer::initializeGlobalCache(const edm::ParameterSet& iConfig) {
+ edm::LogInfo("Initialization") << "Init NN models Global Cache " << std::endl;
+
+ std::unique_ptr GlobalCache(new NNmodels_GlobalCache);
+
+ GlobalCache->CNNmodel_CB_path = iConfig.getParameter("CNNmodel_CB_path");
+ GlobalCache->DNNident_CB_path = iConfig.getParameter("DNNident_CB_path");
+ GlobalCache->DNNcalib_CB_path = iConfig.getParameter("DNNcalib_CB_path");
+ GlobalCache->CNNmodel_CE_path = iConfig.getParameter("CNNmodel_CE_path");
+ GlobalCache->DNNident_CE_path = iConfig.getParameter("DNNident_CE_path");
+ GlobalCache->DNNcalib_CE_path = iConfig.getParameter("DNNcalib_CE_path");
+ GlobalCache->FeatScaler_CE_path = iConfig.getParameter("FeatScaler_CE_path");
+
+ // Create sessions for Tensorflow inferece
+ (GlobalCache->CNNmodel_CB) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->CNNmodel_CB_path)).fullPath());
+ (GlobalCache->CNNmodel_CBsession) = tensorflow::createSession((GlobalCache->CNNmodel_CB));
+
+ (GlobalCache->DNNident_CB) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNident_CB_path)).fullPath());
+ (GlobalCache->DNNident_CBsession) = tensorflow::createSession((GlobalCache->DNNident_CB));
+
+ (GlobalCache->DNNcalib_CB) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNcalib_CB_path)).fullPath());
+ (GlobalCache->DNNcalib_CBsession) = tensorflow::createSession((GlobalCache->DNNcalib_CB));
+
+ (GlobalCache->CNNmodel_CE) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->CNNmodel_CE_path)).fullPath());
+ (GlobalCache->CNNmodel_CEsession) = tensorflow::createSession((GlobalCache->CNNmodel_CE));
+
+ (GlobalCache->DNNident_CE) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNident_CE_path)).fullPath());
+ (GlobalCache->DNNident_CEsession) = tensorflow::createSession((GlobalCache->DNNident_CE));
+
+ (GlobalCache->DNNcalib_CE) = tensorflow::loadGraphDef(edm::FileInPath((GlobalCache->DNNcalib_CE_path)).fullPath());
+ (GlobalCache->DNNcalib_CEsession) = tensorflow::createSession((GlobalCache->DNNcalib_CE));
+
+ // Read features scaler
+ boost::property_tree::read_json(edm::FileInPath((GlobalCache->FeatScaler_CE_path)).fullPath(),
+ (GlobalCache->FeatScaler_CE));
+
+ return GlobalCache;
+}
+
+// ----Constructor and Destructor -----
+L1NNCaloTauProducer::L1NNCaloTauProducer(const edm::ParameterSet& iConfig, const NNmodels_GlobalCache* globalCache)
+ : l1TowersToken(consumes(iConfig.getParameter("l1CaloTowers"))),
+ hgcalTowersToken(consumes(iConfig.getParameter("hgcalTowers"))),
+
+ HGClusterToken(
+ consumes(iConfig.getParameter("HgcalClusters"))),
+ scenario(UseEmInterp::No),
+ preEmId(iConfig.getParameter("preEmId")),
+ VsPuId(iConfig.getParameter("VsPuId")),
+
+ EcalEtMinForClustering(iConfig.getParameter("EcalEtMinForClustering")),
+ HcalEtMinForClustering(iConfig.getParameter("HcalEtMinForClustering")),
+ EtMinForSeeding(iConfig.getParameter("EtMinForSeeding")),
+ EtaRestriction(iConfig.getParameter("EtaRestriction")),
+ CB_CE_split(iConfig.getParameter("CB_CE_split")),
+
+ IdWp90_CB(iConfig.getParameter("IdWp90_CB")),
+ IdWp95_CB(iConfig.getParameter("IdWp95_CB")),
+ IdWp99_CB(iConfig.getParameter("IdWp99_CB")),
+
+ IdWp90_CE(iConfig.getParameter("IdWp90_CE")),
+ IdWp95_CE(iConfig.getParameter("IdWp95_CE")),
+ IdWp99_CE(iConfig.getParameter("IdWp99_CE")) {
+ // Initialize HGCAL BDTs
+ if (!VsPuId.method().empty()) {
+ VsPuId.prepareTMVA();
+ }
+
+ // Create produced outputs
+ produces>("L1NNCaloTauCollectionBXV");
+
+ // Settings output
+ edm::LogInfo("Settings") << "EtaRestriction = " << EtaRestriction << " , CB_CE_split = " << CB_CE_split
+ << " , EtMinForSeeding = " << EtMinForSeeding
+ << " , HcalTpEtMin = " << HcalEtMinForClustering
+ << " , EcalTpEtMin = " << EcalEtMinForClustering << std::endl;
+}
+
+void L1NNCaloTauProducer::produce(edm::Event& iEvent, const edm::EventSetup& eSetup) {
+ // Output collection
+ std::unique_ptr> L1NNCaloTauCollectionBXV(new l1t::TauBxCollection);
+
+ // Create and Fill collection of all calotowers and their attributes
+ std::vector l1CaloTowers;
+
+ iEvent.getByToken(l1TowersToken, l1CaloTowerHandle);
+ int warnings = 0;
+ for (auto& hit : *l1CaloTowerHandle.product()) {
+ // Skip this weird towers and store warning
+ if (hit.towerIEta() == -1016 && hit.towerIPhi() == -962) {
+ warnings += 1;
+ continue;
+ }
+
+ SimpleTowerHit l1Hit;
+ l1Hit.isBarrel = true;
+ l1Hit.l1egTowerEt = hit.l1egTowerEt();
+ l1Hit.towerEta = hit.towerEta();
+ l1Hit.towerPhi = hit.towerPhi();
+ l1Hit.towerEm = hit.ecalTowerEt();
+ l1Hit.towerHad = hit.hcalTowerEt();
+ l1Hit.towerEt = l1Hit.towerEm + l1Hit.towerHad + l1Hit.l1egTowerEt;
+ l1Hit.towerIeta = hit.towerIEta();
+ l1Hit.towerIphi = hit.towerIPhi();
+
+ l1CaloTowers.push_back(l1Hit);
+ }
+ if (warnings != 0) {
+ edm::LogWarning("BrokenTowers") << " ** WARNING : FOUND " << warnings
+ << " TOWERS WITH towerIeta=-1016 AND towerIphi=-962" << std::endl;
+ }
+
+ iEvent.getByToken(hgcalTowersToken, hgcalTowersHandle);
+ for (auto& hit : *hgcalTowersHandle.product()) {
+ SimpleTowerHit l1Hit;
+ l1Hit.isBarrel = false;
+ l1Hit.l1egTowerEt = 0.0;
+ l1Hit.towerEta = hit.eta();
+ l1Hit.towerPhi = hit.phi();
+ l1Hit.towerEm = hit.etEm();
+ l1Hit.towerHad = hit.etHad();
+ l1Hit.towerEt = l1Hit.towerEm + l1Hit.towerHad;
+ l1Hit.towerIeta = endcap_ieta(l1Hit.towerEta); // computed and filled but not used
+ l1Hit.towerIphi = endcap_iphi(l1Hit.towerPhi); // computed and filled but not used
+
+ l1CaloTowers.push_back(l1Hit);
+ }
+
+ // Sort the ECAL+HCAL+L1EGs tower sums based on total ET
+ std::sort(begin(l1CaloTowers), end(l1CaloTowers), [](const SimpleTowerHit& a, SimpleTowerHit& b) {
+ return a.towerEt > b.towerEt;
+ });
+
+ // Create and Fill the collection of 3D clusters and their attributes
+ std::vector AllHGClusters;
+ iEvent.getByToken(HGClusterToken, HGClusterHandle);
+
+ for (auto cl3dIt = HGClusterHandle->begin(0); cl3dIt != HGClusterHandle->end(0); ++cl3dIt) {
+ auto& cl3d = *cl3dIt;
+
+ // Implement cl3d PU ID as done in
+ // https://github.com/cms-sw/cmssw/blob/master/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc#L120
+ bool isEM = preEmId(*cl3dIt);
+ l1t::PFCluster cluster(cl3d.pt(), cl3d.eta(), cl3d.phi(), cl3d.hOverE());
+ if (scenario == UseEmInterp::EmOnly) // for emID objs, use EM interp as pT and set H = 0
+ {
+ if (isEM) {
+ float pt_new = cl3d.iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
+ float hoe_new = 0.;
+ cluster = l1t::PFCluster(pt_new, cl3d.eta(), cl3d.phi(), hoe_new, isEM);
+ }
+ } else if (scenario == UseEmInterp::AllKeepHad) // for all objs, replace EM part with EM interp, preserve H
+ {
+ float had_old = cl3d.pt() - cluster.emEt();
+ float em_new = cl3d.iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
+ float pt_new = had_old + em_new;
+ float hoe_new = em_new > 0 ? (had_old / em_new) : -1;
+ cluster = l1t::PFCluster(pt_new, cl3d.eta(), cl3d.phi(), hoe_new, isEM);
+ } else if (scenario == UseEmInterp::AllKeepTot) // for all objs, replace EM part with EM interp, preserve pT
+ {
+ float em_new = cl3d.iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
+ float hoe_new = em_new > 0 ? (cl3d.pt() / em_new - 1) : -1;
+ cluster = l1t::PFCluster(cl3d.pt(), cl3d.eta(), cl3d.phi(), hoe_new, isEM);
+ }
+
+ if (!VsPuId.method().empty()) {
+ int id = VsPuId.passID(*cl3dIt, cluster);
+ if (!id) {
+ continue;
+ } // skip cl3d if it does not pass puid
+ }
+
+ SimpleHGCluster HGCluster;
+ HGCluster.pt = cl3d.pt();
+ HGCluster.eta = cl3d.eta();
+ HGCluster.phi = cl3d.phi();
+ HGCluster.showerlength = cl3d.showerLength();
+ HGCluster.coreshowerlength = cl3d.coreShowerLength();
+ HGCluster.spptot = cl3d.sigmaPhiPhiTot();
+ HGCluster.szz = cl3d.sigmaZZ();
+ HGCluster.srrtot = cl3d.sigmaRRTot();
+ HGCluster.meanz = cl3d.zBarycenter();
+
+ AllHGClusters.push_back(HGCluster);
+ }
+
+ // Order the collection in pt (the input to the GCT will be pt ordered)
+ std::sort(begin(AllHGClusters), end(AllHGClusters), [](const SimpleHGCluster& a, SimpleHGCluster& b) {
+ return a.pt > b.pt;
+ });
+
+ // Make NxM TowerClusters and HGClusters collections for TauMinator
+ std::vector l1TowerClustersNxM_CB;
+ std::vector l1TowerClustersNxM_CE;
+ std::vector HGClusters;
+
+ // Supporting collection of endcap clusters before cl3d matching
+ std::vector AllL1TowerClustersNxM_CE;
+
+ bool caloTauSeedingFinished = false;
+ // Loop for seeding of clNxM objects
+ while (!caloTauSeedingFinished) {
+ SimpleTowerCluster clNxM;
+ clNxM.InitHits(IEta_dim, IPhi_dim);
+ bool seeded = false;
+
+ for (auto& l1CaloTower : l1CaloTowers) {
+ // Skip seeding in towers that would make the cluster extend in HF
+ // Skip l1CaloTowers which are already used by this clusters' mask
+ if (abs(l1CaloTower.towerEta) > Eta_limit || abs(l1CaloTower.towerEta) > EtaRestriction ||
+ l1CaloTower.stale4seed) {
+ continue;
+ }
+
+ // If not seded do the seeding
+ if (!seeded) {
+ // The leading unused tower has ET < min, stop jet clustering
+ if (l1CaloTower.towerEt < EtMinForSeeding) {
+ caloTauSeedingFinished = true;
+ continue;
+ }
+
+ clNxM.seedIeta = l1CaloTower.towerIeta;
+ clNxM.seedIphi = l1CaloTower.towerIphi;
+ clNxM.seedEta = l1CaloTower.towerEta;
+ clNxM.seedPhi = l1CaloTower.towerPhi;
+ if (l1CaloTower.isBarrel) {
+ clNxM.barrelSeeded = true;
+ }
+
+ clNxM.rawEt += l1CaloTower.towerEt;
+ clNxM.towerHits[seedIdx] = l1CaloTower;
+ l1CaloTower.stale4seed = true;
+ l1CaloTower.stale = true;
+ seeded = true;
+
+ continue;
+ }
+
+ int d_iEta = 99;
+ int d_iPhi = 99;
+ float d_Eta = 99.;
+ float d_Phi = 99.;
+ // Ese iEta/iPhi comparisons in the barrel and eta/phi in HGCal
+ if (clNxM.barrelSeeded && l1CaloTower.isBarrel) {
+ d_iEta = tower_dIEta(l1CaloTower.towerIeta, clNxM.seedIeta);
+ d_iPhi = tower_dIPhi(l1CaloTower.towerIphi, clNxM.seedIphi);
+ } else {
+ d_Eta = l1CaloTower.towerEta - clNxM.seedEta;
+ d_Phi = reco::deltaPhi(l1CaloTower.towerPhi, clNxM.seedPhi);
+ }
+
+ // Stale tower for seeding if it would lead to overalp between clusters
+ if ((abs(d_iEta) <= IEta_dim - 1 && abs(d_iPhi) <= IPhi_dim - 1) ||
+ (abs(d_Eta) < Eta_dim_seed && abs(d_Phi) < Phi_dim_seed)) {
+ l1CaloTower.stale4seed = true;
+ }
+
+ } // End for loop over TPs
+
+ // Pushback seeds split in barrel and endcap
+ if (seeded) {
+ if (abs(clNxM.seedEta) < CB_CE_split) {
+ l1TowerClustersNxM_CB.push_back(clNxM);
+ } else {
+ AllL1TowerClustersNxM_CE.push_back(clNxM);
+ }
+ }
+
+ } // End while loop of TowerClusters seeding
+
+ // Loop for barrel NxM TowerClusters clustering starting from the seeds
+ for (auto& clNxM : l1TowerClustersNxM_CB) {
+ for (auto& l1CaloTower : l1CaloTowers) {
+ // Skip l1CaloTowers which are already used
+ if (l1CaloTower.stale) {
+ continue;
+ }
+
+ int d_iEta = 99;
+ int d_iPhi = 99;
+ float d_Eta = 99.;
+ float d_Phi = 99.;
+ int hitIdx = 99.;
+ // Use iEta/iPhi comparisons in the barrel and use eta/phi in HGCal
+ if (l1CaloTower.isBarrel) {
+ d_iEta = tower_dIEta(l1CaloTower.towerIeta, clNxM.seedIeta);
+ d_iPhi = tower_dIPhi(l1CaloTower.towerIphi, clNxM.seedIphi);
+
+ hitIdx = d_iEta * IPhi_dim + d_iPhi + seedIdx;
+ } else {
+ d_Eta = l1CaloTower.towerEta - clNxM.seedEta;
+ d_Phi = reco::deltaPhi(l1CaloTower.towerPhi, clNxM.seedPhi);
+
+ int dieta = d_Eta / 0.0807; // minimal difference in endcap is 0.0808
+ int diphi = d_Phi / 0.0872;
+ hitIdx = dieta * IPhi_dim + diphi + seedIdx;
+ }
+
+ // Cluster all towers in a NxM towers mask
+ if ((abs(d_iEta) <= (IEta_dim - 1) / 2 && abs(d_iPhi) <= (IPhi_dim - 1) / 2) ||
+ (abs(d_Eta) < Eta_dim && abs(d_Phi) < Phi_dim)) {
+ clNxM.rawEt += l1CaloTower.towerEt;
+ clNxM.towerHits[hitIdx] = l1CaloTower;
+ l1CaloTower.stale = true;
+ }
+
+ } // End for loop over TPs
+
+ } // End while loop of barrel TowerClusters creation
+
+ // In the endcap cross-loop over clNxM and cl3d to match them
+ // (we can do it before full clustering just using the seed info)
+ for (auto& clNxM : AllL1TowerClustersNxM_CE) {
+ bool matched = false;
+ for (auto& HGCluster : AllHGClusters) {
+ // In case the clNxM or HGCluster have already been matched just continue through the list to the end
+ // only use cl3ds above 4GeV
+ if (matched || HGCluster.stale || HGCluster.pt < 4) {
+ continue;
+ }
+
+ float d_Eta = HGCluster.eta - clNxM.seedEta;
+ float d_Phi = reco::deltaPhi(HGCluster.phi, clNxM.seedPhi);
+ float d_R2 = pow(d_Eta, 2) + pow(d_Phi, 2);
+
+ if (d_R2 < 0.25) {
+ HGCluster.stale = true;
+ HGClusters.push_back(HGCluster);
+ l1TowerClustersNxM_CE.push_back(clNxM);
+ matched = true;
+ }
+
+ } // End for loop over cl3ds
+
+ } // End for loop over clNxM
+
+ // Loop for endcap matched NxM TowerClusters clustering starting from the seeds just found
+ for (auto& clNxM : l1TowerClustersNxM_CE) {
+ for (auto& l1CaloTower : l1CaloTowers) {
+ // Skip l1CaloTowers which are already used
+ if (l1CaloTower.stale) {
+ continue;
+ }
+
+ int d_iEta = 99;
+ int d_iPhi = 99;
+ float d_Eta = 99.;
+ float d_Phi = 99.;
+ int hitIdx = 99.;
+ // Use iEta/iPhi comparisons in the endcap and use eta/phi in HGCal
+ if (l1CaloTower.isBarrel) {
+ d_iEta = tower_dIEta(l1CaloTower.towerIeta, clNxM.seedIeta);
+ d_iPhi = tower_dIPhi(l1CaloTower.towerIphi, clNxM.seedIphi);
+
+ hitIdx = d_iEta * IPhi_dim + d_iPhi + seedIdx;
+ } else {
+ d_Eta = l1CaloTower.towerEta - clNxM.seedEta;
+ d_Phi = reco::deltaPhi(l1CaloTower.towerPhi, clNxM.seedPhi);
+
+ int dieta = d_Eta / 0.0807; // minimal difference in endcap is 0.0808
+ int diphi = d_Phi / 0.0872;
+ hitIdx = dieta * IPhi_dim + diphi + seedIdx;
+ }
+
+ // Cluster all towers in a NxM towers mask
+ if ((abs(d_iEta) <= (IEta_dim - 1) / 2 && abs(d_iPhi) <= (IPhi_dim - 1) / 2) ||
+ (abs(d_Eta) < Eta_dim && abs(d_Phi) < Phi_dim)) {
+ clNxM.rawEt += l1CaloTower.towerEt;
+ clNxM.towerHits[hitIdx] = l1CaloTower;
+ l1CaloTower.stale = true;
+ }
+
+ } // End for loop over TPs
+
+ } // End while loop of endcap TowerClusters creation
+
+ // Barrel TauMinator application
+ tensorflow::setLogging("2");
+ int batchSize_CB = (int)(l1TowerClustersNxM_CB.size());
+ tensorflow::TensorShape imageShape_CB({batchSize_CB, IEta_dim, IPhi_dim, 2});
+ tensorflow::TensorShape positionShape_CB({batchSize_CB, 2});
+ tensorflow::Tensor TowerClusterImage_CB(tensorflow::DT_FLOAT, imageShape_CB);
+ tensorflow::Tensor TowerClusterPosition_CB(tensorflow::DT_FLOAT, positionShape_CB);
+
+ int clIdx = 0;
+ for (auto& clNxM : l1TowerClustersNxM_CB) {
+ // Fill inputs for Tensorflow inference
+ for (int eta = 0; eta < IEta_dim; ++eta) {
+ for (int phi = 0; phi < IPhi_dim; ++phi) {
+ int towerIdx = eta * IPhi_dim + phi;
+ TowerClusterImage_CB.tensor()(clIdx, eta, phi, 0) =
+ inputQuantizer(clNxM.towerHits[towerIdx].l1egTowerEt + clNxM.towerHits[towerIdx].towerEm, 0.25, 10);
+ TowerClusterImage_CB.tensor()(clIdx, eta, phi, 1) =
+ inputQuantizer(clNxM.towerHits[towerIdx].towerHad, 0.25, 10);
+ }
+ }
+
+ TowerClusterPosition_CB.tensor()(clIdx, 0) = clNxM.seedEta;
+ TowerClusterPosition_CB.tensor()(clIdx, 1) = clNxM.seedPhi;
+
+ clIdx++; // Increase batch index
+ }
+
+ if (batchSize_CB >
+ 0) // from CMSSW_14_0_X tensorflow does not seem to be able to deal with a tensor of dimension 0 anymore
+ {
+ // Apply CNN model
+ tensorflow::NamedTensorList CNNmodel_CBinputList = {{"TowerClusterImage", TowerClusterImage_CB},
+ {"TowerClusterPosition", TowerClusterPosition_CB}};
+ std::vector CNNmodel_CBoutputs;
+ tensorflow::run((globalCache()->CNNmodel_CBsession),
+ CNNmodel_CBinputList,
+ {"TauMinator_CB_conv/middleMan/concat"},
+ &CNNmodel_CBoutputs);
+ tensorflow::NamedTensorList DNN_CBinputsList = {{"middleMan", CNNmodel_CBoutputs[0]}};
+
+ // Apply DNN for identification
+ std::vector DNN_CBoutputsIdent;
+ tensorflow::run((globalCache()->DNNident_CBsession),
+ DNN_CBinputsList,
+ {"TauMinator_CB_ident/sigmoid_IDout/Sigmoid"},
+ &DNN_CBoutputsIdent);
+
+ // Apply DNN for calibration
+ std::vector DNN_CBoutputsCalib;
+ tensorflow::run((globalCache()->DNNcalib_CBsession),
+ DNN_CBinputsList,
+ {"TauMinator_CB_calib/LIN_DNNout/Relu"},
+ &DNN_CBoutputsCalib);
+
+ // Fill TauMinator output variables of TowerClusters
+ clIdx = 0;
+ for (auto& clNxM : l1TowerClustersNxM_CB) {
+ clNxM.IDscore = DNN_CBoutputsIdent[0].matrix()(0, clIdx);
+ clNxM.calibPt = DNN_CBoutputsCalib[0].matrix()(0, clIdx);
+ clIdx++; // Increase batch index
+ }
+ }
+
+ // Endcap TauMinator application
+ int batchSize_CE = (int)(l1TowerClustersNxM_CE.size());
+ tensorflow::TensorShape imageShape_CE({batchSize_CE, IEta_dim, IPhi_dim, 2});
+ tensorflow::TensorShape positionShape_CE({batchSize_CE, 2});
+ tensorflow::TensorShape cl3dfeatShape_CE({batchSize_CE, 8});
+ tensorflow::Tensor TowerClusterImage_CE(tensorflow::DT_FLOAT, imageShape_CE);
+ tensorflow::Tensor TowerClusterPosition_CE(tensorflow::DT_FLOAT, positionShape_CE);
+ tensorflow::Tensor Cl3dShapeFeatures_CE(tensorflow::DT_FLOAT, cl3dfeatShape_CE);
+
+ clIdx = 0;
+ for (auto& clNxM : l1TowerClustersNxM_CE) {
+ // Indexing of cl3ds is the same as the one of clNxMs
+ SimpleHGCluster HGClu = HGClusters[clIdx];
+
+ // Fill inputs for Tensorflow inference
+ for (int eta = 0; eta < IEta_dim; ++eta) {
+ for (int phi = 0; phi < IPhi_dim; ++phi) {
+ int towerIdx = eta * IPhi_dim + phi;
+ TowerClusterImage_CE.tensor()(clIdx, eta, phi, 0) =
+ inputQuantizer(clNxM.towerHits[towerIdx].l1egTowerEt + clNxM.towerHits[towerIdx].towerEm, 0.25, 10);
+ TowerClusterImage_CE.tensor()(clIdx, eta, phi, 1) =
+ inputQuantizer(clNxM.towerHits[towerIdx].towerHad, 0.25, 10);
+ }
+ }
+
+ TowerClusterPosition_CE.tensor()(clIdx, 0) = clNxM.seedEta;
+ TowerClusterPosition_CE.tensor()(clIdx, 1) = clNxM.seedPhi;
+
+ Cl3dShapeFeatures_CE.tensor()(clIdx, 0) = inputScaler(inputQuantizer(HGClu.pt, 0.25, 14), "pt");
+ Cl3dShapeFeatures_CE.tensor()(clIdx, 1) =
+ inputScaler(inputQuantizer(abs(HGClu.eta) - 1.321, 0.004, 9), "eta");
+ Cl3dShapeFeatures_CE.tensor()(clIdx, 2) = inputScaler(HGClu.showerlength, "showerlength");
+ Cl3dShapeFeatures_CE.tensor()(clIdx, 3) = inputScaler(HGClu.coreshowerlength, "coreshowerlength");
+ Cl3dShapeFeatures_CE.tensor()(clIdx, 4) =
+ inputScaler(inputQuantizer(HGClu.spptot, 0.0000153, 16), "spptot");
+ Cl3dShapeFeatures_CE.tensor()(clIdx, 5) = inputScaler(inputQuantizer(HGClu.szz, 0.00153, 16), "szz");
+ Cl3dShapeFeatures_CE.tensor()(clIdx, 6) =
+ inputScaler(inputQuantizer(HGClu.srrtot, 0.0000153, 16), "srrtot");
+ Cl3dShapeFeatures_CE.tensor()(clIdx, 7) =
+ inputScaler(inputQuantizer(10 * (abs(HGClu.meanz) - 321.05), 0.5, 12), "meanz");
+
+ clIdx++; // Increase batch index
+ }
+
+ if (batchSize_CE >
+ 0) // from CMSSW_14_0_X tensorflow does not seem to be able to deal with a tensor of dimension 0 anymore
+ {
+ // Apply CNN model
+ tensorflow::NamedTensorList CNNmodel_CEinputList = {{"TowerClusterImage", TowerClusterImage_CE},
+ {"TowerClusterPosition", TowerClusterPosition_CE},
+ {"AssociatedCl3dFeatures", Cl3dShapeFeatures_CE}};
+ std::vector CNNmodel_CEoutputs;
+ tensorflow::run((globalCache()->CNNmodel_CEsession),
+ CNNmodel_CEinputList,
+ {"TauMinator_CE_conv/middleMan/concat"},
+ &CNNmodel_CEoutputs);
+ tensorflow::NamedTensorList DNN_CEinputsList = {{"middleMan", CNNmodel_CEoutputs[0]}};
+
+ // Apply DNN for identification
+ std::vector DNN_CEoutputsIdent;
+ tensorflow::run((globalCache()->DNNident_CEsession),
+ DNN_CEinputsList,
+ {"TauMinator_CE_ident/sigmoid_IDout/Sigmoid"},
+ &DNN_CEoutputsIdent);
+
+ // Apply DNN for calibration
+ std::vector DNN_CEoutputsCalib;
+ tensorflow::run((globalCache()->DNNcalib_CEsession),
+ DNN_CEinputsList,
+ {"TauMinator_CE_calib/LIN_DNNout/Relu"},
+ &DNN_CEoutputsCalib);
+
+ // Fill TauMinator output variables of TowerClusters
+ clIdx = 0;
+ for (auto& clNxM : l1TowerClustersNxM_CE) {
+ clNxM.IDscore = DNN_CEoutputsIdent[0].matrix()(0, clIdx);
+ clNxM.calibPt = DNN_CEoutputsCalib[0].matrix()(0, clIdx);
+ clIdx++; // Increase batch index
+ }
+ }
+
+ // Fill the output collection of L1 taus
+ for (auto& clNxM : l1TowerClustersNxM_CB) {
+ // Apply eta restriction
+ if (abs(clNxM.seedEta) > EtaRestriction) {
+ continue;
+ }
+
+ // Assign increasing quality to higher scoring candidates
+ int quality = 0;
+ // 99% WP
+ if (clNxM.IDscore > IdWp99_CB) {
+ quality = 1;
+ }
+ // 95% WP
+ if (clNxM.IDscore > IdWp95_CB) {
+ quality = 2;
+ }
+ // 90% WP
+ if (clNxM.IDscore > IdWp90_CB) {
+ quality = 3;
+ }
+
+ reco::Candidate::PolarLorentzVector tauP4 =
+ reco::Candidate::PolarLorentzVector(clNxM.calibPt, clNxM.seedEta, clNxM.seedPhi, 0);
+
+ // store ID score multiplied by 10E4 to have good precision even using the Phase1 tau int iso format
+ // (this is stored just in case for possible additional offline studies)
+ // tau initialisation = (p4, pt, eta, phi, qual, iso)
+ l1t::Tau l1Tau = l1t::Tau(tauP4, clNxM.calibPt, clNxM.seedEta, clNxM.seedPhi, quality, clNxM.IDscore * 10E4);
+ l1Tau.setTowerIEta(clNxM.seedIeta);
+ l1Tau.setTowerIPhi(clNxM.seedIphi);
+ l1Tau.setRawEt(clNxM.rawEt);
+
+ L1NNCaloTauCollectionBXV->push_back(0, l1Tau);
+ }
+
+ for (auto& clNxM : l1TowerClustersNxM_CE) {
+ // Apply eta restriction
+ if (abs(clNxM.seedEta) > EtaRestriction) {
+ continue;
+ }
+
+ // Assign increasing quality to higher scoring candidates
+ int quality = 0;
+ // 99% WP
+ if (clNxM.IDscore > IdWp99_CE) {
+ quality = 1;
+ }
+ // 95% WP
+ if (clNxM.IDscore > IdWp95_CE) {
+ quality = 2;
+ }
+ // 90% WP
+ if (clNxM.IDscore > IdWp90_CE) {
+ quality = 3;
+ }
+
+ reco::Candidate::PolarLorentzVector tauP4 =
+ reco::Candidate::PolarLorentzVector(clNxM.calibPt, clNxM.seedEta, clNxM.seedPhi, 0);
+
+ // store ID score multiplied by 10E4 to have good precision even using the Phase1 tau int iso format
+ // (this is stored just in case for possible additional offline studies)
+ // tau initialisation = (p4, pt, eta, phi, qual, iso)
+ l1t::Tau l1Tau = l1t::Tau(tauP4, clNxM.calibPt, clNxM.seedEta, clNxM.seedPhi, quality, clNxM.IDscore * 10E4);
+ l1Tau.setTowerIEta(clNxM.seedIeta);
+ l1Tau.setTowerIPhi(clNxM.seedIphi);
+ l1Tau.setRawEt(clNxM.rawEt);
+
+ L1NNCaloTauCollectionBXV->push_back(0, l1Tau);
+ }
+
+ // Fill output
+ iEvent.put(std::move(L1NNCaloTauCollectionBXV), "L1NNCaloTauCollectionBXV");
+
+} // End of produce function
+
+int L1NNCaloTauProducer::tower_dIPhi(int& iPhi_1, int& iPhi_2) const {
+ const int PI = 36;
+ int result = iPhi_1 - iPhi_2;
+ if (result > PI) {
+ result -= 2 * PI;
+ }
+ if (result <= -PI) {
+ result += 2 * PI;
+ }
+ return result;
+}
+
+int L1NNCaloTauProducer::tower_dIEta(int& iEta_1, int& iEta_2) const {
+ if (iEta_1 * iEta_2 > 0) {
+ return iEta_1 - iEta_2;
+ } else {
+ if (iEta_1 > 0) {
+ return iEta_1 - iEta_2 - 1;
+ } else {
+ return iEta_1 - iEta_2 + 1;
+ }
+ }
+}
+
+int L1NNCaloTauProducer::endcap_iphi(float& phi) const {
+ const float phi_step = 0.0872664;
+ if (phi > 0) {
+ return floor(phi / phi_step) + 1;
+ } else {
+ return floor(phi / phi_step) + 73;
+ }
+}
+
+int L1NNCaloTauProducer::endcap_ieta(float& eta) const {
+ const float eta_step = 0.0845;
+ return floor(abs(eta) / eta_step) * std::copysign(1, eta);
+}
+
+float L1NNCaloTauProducer::inputQuantizer(float inputF, float LSB, int nbits) {
+ return min(floor(inputF / LSB), float(pow(2, nbits) - 1)) * LSB;
+}
+
+float L1NNCaloTauProducer::inputScaler(float inputF, std::string feature) {
+ float mean = (globalCache()->FeatScaler_CE).get_child(feature).get("mean");
+ float std = (globalCache()->FeatScaler_CE).get_child(feature).get("std");
+
+ return (inputF - mean) / std;
+}
+
+void L1NNCaloTauProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+ edm::ParameterSetDescription desc;
+
+ desc.add("l1CaloTowers", edm::InputTag("l1tEGammaClusterEmuProducer", "L1CaloTowerCollection"));
+ desc.add("hgcalTowers", edm::InputTag("l1tHGCalTowerProducer", "HGCalTowerProcessor"));
+ desc.add("HgcalClusters",
+ edm::InputTag("l1tHGCalBackEndLayer2Producer", "HGCalBackendLayer2Processor3DClustering"));
+
+ desc.add("preEmId", "hOverE < 0.3 && hOverE >= 0");
+ {
+ edm::ParameterSetDescription psd0;
+ psd0.add("isPUFilter", true);
+ psd0.add("preselection", "");
+ psd0.add("method", "BDT");
+ {
+ edm::ParameterSetDescription vpsd2;
+ vpsd2.add("name");
+ vpsd2.add("value");
+ std::vector temp2;
+ temp2.reserve(5);
+ {
+ edm::ParameterSet temp3;
+ temp3.addParameter("name", "eMax");
+ temp3.addParameter("value", "eMax()");
+ temp2.push_back(temp3);
+ }
+ {
+ edm::ParameterSet temp3;
+ temp3.addParameter("name", "eMaxOverE");
+ temp3.addParameter("value", "eMax()/energy()");
+ temp2.push_back(temp3);
+ }
+ {
+ edm::ParameterSet temp3;
+ temp3.addParameter("name", "sigmaPhiPhiTot");
+ temp3.addParameter("value", "sigmaPhiPhiTot()");
+ temp2.push_back(temp3);
+ }
+ {
+ edm::ParameterSet temp3;
+ temp3.addParameter("name", "sigmaRRTot");
+ temp3.addParameter("value", "sigmaRRTot()");
+ temp2.push_back(temp3);
+ }
+ {
+ edm::ParameterSet temp3;
+ temp3.addParameter("name", "triggerCells90percent");
+ temp3.addParameter