From 6ba8576f3fe2f47587ecb19c7b06b37a2de8bc69 Mon Sep 17 00:00:00 2001 From: Jona Motta Date: Thu, 21 Sep 2023 10:00:27 +0200 Subject: [PATCH] upload TauMinator (NN Calo Tau) code --- .../python/L1Trigger_EventContent_cff.py | 2 + .../Configuration/python/SimL1Emulator_cff.py | 8 +- L1Trigger/L1CaloTrigger/plugins/BuildFile.xml | 3 + .../plugins/L1NNCaloTauEmulator.cc | 1065 +++++++++++++++++ .../plugins/L1NNCaloTauProducer.cc | 809 +++++++++++++ .../python/l1tNNCaloTauEmulator_cff.py | 5 + .../python/l1tNNCaloTauEmulator_cfi.py | 50 + .../python/l1tNNCaloTauProducer_cff.py | 5 + .../python/l1tNNCaloTauProducer_cfi.py | 49 + 9 files changed, 1995 insertions(+), 1 deletion(-) create mode 100644 L1Trigger/L1CaloTrigger/plugins/L1NNCaloTauEmulator.cc create mode 100644 L1Trigger/L1CaloTrigger/plugins/L1NNCaloTauProducer.cc create mode 100644 L1Trigger/L1CaloTrigger/python/l1tNNCaloTauEmulator_cff.py create mode 100644 L1Trigger/L1CaloTrigger/python/l1tNNCaloTauEmulator_cfi.py create mode 100644 L1Trigger/L1CaloTrigger/python/l1tNNCaloTauProducer_cff.py create mode 100644 L1Trigger/L1CaloTrigger/python/l1tNNCaloTauProducer_cfi.py 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..55d010d5c8b73 --- /dev/null +++ b/L1Trigger/L1CaloTrigger/plugins/L1NNCaloTauEmulator.cc @@ -0,0 +1,1065 @@ +/* -*- C++ -*- + +Package: L1CaloTrigger +Class: l1tNNCaloTauEmulator +Frinedly name: The TauMinator + +\class l1tNNCaloTauEmulator l1tNNCaloTauEmulator.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/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" + +class l1tNNCaloTauEmulator : public edm::stream::EDProducer<> { +public: + explicit l1tNNCaloTauEmulator(const edm::ParameterSet&); + +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---- + dIEtaPhi_t tower_dIPhi(IPhi_t iPhi_1, IPhi_t iPhi_2); + dIEtaPhi_t tower_dIEta(IEta_t iEta_1, IEta_t iEta_2); + dEtaPhi_t dPhi(EtaPhi_t iPhi_1, EtaPhi_t iPhi_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; + + 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; + + 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.; + // IEta_t ieta = 0; // DEBUG ONLY + // IPhi_t iphi = 0; // DEBUG ONLY + + void fill(SimpleTowerHit Tower) { + towerEm = Tower.towerEm; + towerHad = Tower.towerHad; + l1egTowerEt = Tower.l1egTowerEt; + // ieta = Tower.towerIeta; // DEBUG ONLY + // iphi = Tower.towerIphi; // DEBUG ONLY + } + }; + + 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; + // EtaPhi_t phi; // DEBUG ONLY + 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; + // phi = Cluster.phi; // DEBUG ONLY + showerlength = Cluster.showerlength; + coreshowerlength = Cluster.coreshowerlength; + spptot = Cluster.spptot; + szz = Cluster.szz; + srrtot = Cluster.srrtot; + meanz = Cluster.meanz; + } + }; + + /* // output structure of taus + struct NNCaloTau { + Et_t hwPt; + Et_t hwSeedPt; + IEta_t hwEta; + IPhi_t hwPhi; + + ap_uint<2> quality; + ap_uint<9> IDscore; + + inline CALOtoGT_pt(Et_t x) { return (l1gt::pt_t)x; } + inline CALOtoGT_eta(IEta_t x) { return x * (IETAPHI_LSB / l1gt::Scales::ETAPHI_LSB); } + inline CALOtoGT_phi(IPhi_t x) { return x * (IETAPHI_LSB / l1gt::Scales::ETAPHI_LSB); } + + l1gt::Tau toGT() const { + l1gt::Tau tau; + tau.valid = hwPt != 0; + tau.v3.pt = CALOtoGT_pt(hwPt); + tau.v3.phi = CALOtoGT_phi(hwPhi); + tau.v3.eta = CALOtoGT_eta(hwEta); + tau.seed_pt = CALOtoGT_pt(hwSeedPt); + tau.seed_z0 = 0; + tau.charge = 0; + tau.type = 0; + tau.isolation = 0; + tau.id0 = quality; + tau.id1 = 0; + } + + }; */ +}; + +/* +████████ ██ ██ ██████ ████████ █████ ██ ██ ███ ███ ██ ███ ██ █████ ████████ ██████ ██████ + ██ ██ ██ ██ ██ ██ ██ ██ ██ ████ ████ ██ ████ ██ ██ ██ ██ ██ ██ ██ ██ + ██ ███████ █████ ██ ███████ ██ ██ ██ ████ ██ ██ ██ ██ ██ ███████ ██ ██ ██ ██████ + ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ + ██ ██ ██ ██████ ██ ██ ██ ███████ ██ ██ ██ ██ ████ ██ ██ ██ ██████ ██ ██ +*/ + +// ----Constructor and Destructor ----- +l1tNNCaloTauEmulator::l1tNNCaloTauEmulator(const edm::ParameterSet& iConfig) + : 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")), + + CNNmodel_CB_path(iConfig.getParameter("CNNmodel_CB_path")), + DNNident_CB_path(iConfig.getParameter("DNNident_CB_path")), + DNNcalib_CB_path(iConfig.getParameter("DNNcalib_CB_path")), + CNNmodel_CE_path(iConfig.getParameter("CNNmodel_CE_path")), + DNNident_CE_path(iConfig.getParameter("DNNident_CE_path")), + DNNcalib_CE_path(iConfig.getParameter("DNNcalib_CE_path")), + FeatScaler_CE_path(iConfig.getParameter("FeatScaler_CE_path")), + + 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")) { + // Create sessions for Tensorflow inferece + CNNmodel_CB = tensorflow::loadGraphDef(edm::FileInPath(CNNmodel_CB_path).fullPath()); + CNNmodel_CBsession = tensorflow::createSession(CNNmodel_CB); + + DNNident_CB = tensorflow::loadGraphDef(edm::FileInPath(DNNident_CB_path).fullPath()); + DNNident_CBsession = tensorflow::createSession(DNNident_CB); + + DNNcalib_CB = tensorflow::loadGraphDef(edm::FileInPath(DNNcalib_CB_path).fullPath()); + DNNcalib_CBsession = tensorflow::createSession(DNNcalib_CB); + + CNNmodel_CE = tensorflow::loadGraphDef(edm::FileInPath(CNNmodel_CE_path).fullPath()); + CNNmodel_CEsession = tensorflow::createSession(CNNmodel_CE); + + DNNident_CE = tensorflow::loadGraphDef(edm::FileInPath(DNNident_CE_path).fullPath()); + DNNident_CEsession = tensorflow::createSession(DNNident_CE); + + DNNcalib_CE = tensorflow::loadGraphDef(edm::FileInPath(DNNcalib_CE_path).fullPath()); + DNNcalib_CEsession = tensorflow::createSession(DNNcalib_CE); + + // Read features scaler + boost::property_tree::read_json(edm::FileInPath(FeatScaler_CE_path).fullPath(), FeatScaler_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 + std::cout << "EtaRestriction = " << EtaRestriction << " (" << intEtaRestriction << ")" + << " , CB_CE_split = " << CB_CE_split << "(" << intCB_CE_split + << ") , EtMinForSeeding = " << EtMinForSeeding << " , HcalTpEtMin = " << HcalEtMinForClustering + << " , EcalTpEtMin = " << EcalEtMinForClustering << " , PuidThr = " << PuidThr << "(" << intPuidThr << ")" + << std::endl; +} + +void l1tNNCaloTauEmulator::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) { + std::cout << " ** 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()) { + int id = 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 = tower_dIPhi(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 = tower_dIPhi(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 = tower_dIPhi(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); + } + + // Apply CNN model + tensorflow::NamedTensorList CNNmodel_CBinputList = {{"TowerClusterImage", TowerClusterImage_CB}, + {"TowerClusterPosition", TowerClusterPosition_CB}}; + std::vector CNNmodel_CBoutputs; + tensorflow::run( + 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( + DNNident_CBsession, DNN_CBinputsList, {"TauMinator_CB_ident/sigmoid_IDout/Sigmoid"}, &DNN_CBoutputsIdent); + + // Apply DNN for calibration + std::vector DNN_CBoutputsCalib; + tensorflow::run(DNNcalib_CBsession, DNN_CBinputsList, {"TauMinator_CB_calib/DNNout/MatMul"}, &DNN_CBoutputsCalib); + + // 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"); + } + + // Apply CNN model + tensorflow::NamedTensorList CNNmodel_CEinputList = {{"TowerClusterImage", TowerClusterImage_CE}, + {"TowerClusterPosition", TowerClusterPosition_CE}, + {"AssociatedCl3dFeatures", Cl3dShapeFeatures_CE}}; + std::vector CNNmodel_CEoutputs; + tensorflow::run( + 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( + DNNident_CEsession, DNN_CEinputsList, {"TauMinator_CE_ident/sigmoid_IDout/Sigmoid"}, &DNN_CEoutputsIdent); + + // Apply DNN for calibration + std::vector DNN_CEoutputsCalib; + tensorflow::run(DNNcalib_CEsession, DNN_CEinputsList, {"TauMinator_CE_calib/LIN_DNNout/Relu"}, &DNN_CEoutputsCalib); + + // ------------------------------------------------------------- */ + // RESTART OF SOFTWARE PRECISION SECTION + // from here on we go back to floating point precision to + // produce the output for Ntuplization and further work, + // and the output for the GT. + // * + + // Fill the output collection of L1 taus with the barrel candidates + for (int clNxMIdx = 0; clNxMIdx < Nclusters_CB; clNxMIdx++) { + int seedIeta = l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIeta; + int seedIphi = l1TowerClustersNxM_CB_pstn[clNxMIdx].seedIphi; + + if (seedIeta > intEtaRestriction) { + continue; + } + + float tau_IDscore = DNN_CBoutputsIdent[0].matrix()(0, clNxMIdx); + float tau_calibPt = DNN_CBoutputsCalib[0].matrix()(0, clNxMIdx); + float tau_eta = floatIEta(seedIeta); + float tau_phi = floatIPhi(seedIphi); + + // Assign increasing quality to higher scoring candidates + int quality = 0; + // 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; + } + + 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); + // l1Tau.setRawEt(clNxM.rawEt); // FIXME : not really needed tbh + + L1NNCaloTauCollectionBXV->push_back(0, l1Tau); + } + + // Fill the output collection of L1 taus with the endcap candidates + for (int clNxMIdx = 0; clNxMIdx < Nclusters_CE; clNxMIdx++) { + int seedIeta = l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIeta; + int seedIphi = l1TowerClustersNxM_CE_pstn[clNxMIdx].seedIphi; + + if (seedIeta > intEtaRestriction) { + continue; + } + + float tau_IDscore = DNN_CEoutputsIdent[0].matrix()(0, clNxMIdx); + float tau_calibPt = DNN_CEoutputsCalib[0].matrix()(0, clNxMIdx); + float tau_eta = floatIEta(seedIeta); + float tau_phi = floatIPhi(seedIphi); + + // Assign increasing quality to higher scoring candidates + int quality = 0; + // 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); + // l1Tau.setRawEt(clNxM.rawEt); // FIXME : not really needed tbh + + L1NNCaloTauCollectionBXV->push_back(0, l1Tau); + } + + // Fill output + iEvent.put(std::move(L1NNCaloTauCollectionBXV), "L1NNCaloTauCollectionBXV"); + +} // End of produce function + +l1tNNCaloTauEmulator::dIEtaPhi_t l1tNNCaloTauEmulator::tower_dIPhi(IPhi_t iPhi_1, IPhi_t iPhi_2) { + dIEtaPhi_t dphi = iPhi_1 - iPhi_2; + + dIEtaPhi_t dphi0 = dphi > dIEtaPhi_t(INTPHI_PI) ? dIEtaPhi_t(dphi - INTPHI_2PI) : dphi; + dIEtaPhi_t dphi1 = dphi <= dIEtaPhi_t(-INTPHI_PI) ? dIEtaPhi_t(dphi + INTPHI_2PI) : dphi; + + dIEtaPhi_t result = dphi > dIEtaPhi_t(0) ? dphi0 : dphi1; + + return result; +} + +l1tNNCaloTauEmulator::dIEtaPhi_t l1tNNCaloTauEmulator::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; +} + +l1tNNCaloTauEmulator::dEtaPhi_t l1tNNCaloTauEmulator::dPhi(EtaPhi_t iPhi_1, EtaPhi_t iPhi_2) { + dEtaPhi_t dphi = iPhi_1 - iPhi_2; + + dEtaPhi_t dphi0 = dphi > dEtaPhi_t(FINEINTPHI_PI) ? dEtaPhi_t(FINEINTPHI_2PI - dphi) : dphi; + dEtaPhi_t dphi1 = dphi < dEtaPhi_t(-FINEINTPHI_PI) ? dEtaPhi_t(FINEINTPHI_2PI + dphi) : dphi; + + dEtaPhi_t result = dphi > dEtaPhi_t(0) ? dphi0 : dphi1; + + return result; +} + +l1tNNCaloTauEmulator::dEtaPhi_t l1tNNCaloTauEmulator::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); +} + +l1tNNCaloTauEmulator::dEtaPhi_t l1tNNCaloTauEmulator::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; +} + +l1tNNCaloTauEmulator::IEta_t l1tNNCaloTauEmulator::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; +} + +l1tNNCaloTauEmulator::IPhi_t l1tNNCaloTauEmulator::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 l1tNNCaloTauEmulator::ap_abs(ap_int x) { + ap_int result; + if (x < 0) { + result = -x; + } else { + result = x; + } + + return result; +} + +template +ap_ufixed l1tNNCaloTauEmulator::ap_abs(ap_fixed x) { + ap_ufixed result; + if (x < 0) { + result = -x; + } else { + result = x; + } + + return result; +} + +float l1tNNCaloTauEmulator::apfixedQuantizer(float inputF, float LSB, int nbits) { + return min(floor(inputF / LSB), float(pow(2, nbits) - 1)) * LSB; +} + +int l1tNNCaloTauEmulator::apintQuantizer(float inputF, float LSB, int nbits) { + return min(floor(inputF / LSB), float(pow(2, nbits) - 1)); +} + +float l1tNNCaloTauEmulator::inputScaler(float inputF, std::string feature) { + float mean = FeatScaler_CE.get_child(feature).get("mean"); + float std = FeatScaler_CE.get_child(feature).get("std"); + + return (inputF - mean) / std; +} + +float l1tNNCaloTauEmulator::correctInputEtaCl3d(float eta) { + return eta > 0 ? eta - ETAHGCAL_OFFSET : eta + ETAHGCAL_OFFSET; +} + +float l1tNNCaloTauEmulator::correctInputMeanzCl3d(float meanz) { return CM2MM * (abs(meanz) - MEANZ_OFFSET); } + +float l1tNNCaloTauEmulator::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 l1tNNCaloTauEmulator::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; +} + +DEFINE_FWK_MODULE(l1tNNCaloTauEmulator); diff --git a/L1Trigger/L1CaloTrigger/plugins/L1NNCaloTauProducer.cc b/L1Trigger/L1CaloTrigger/plugins/L1NNCaloTauProducer.cc new file mode 100644 index 0000000000000..27c79ffedbc94 --- /dev/null +++ b/L1Trigger/L1CaloTrigger/plugins/L1NNCaloTauProducer.cc @@ -0,0 +1,809 @@ +/* -*- C++ -*- + +Package: L1CaloTrigger +Class: l1tNNCaloTauProducer +Frinedly name: The TauMinator + +\class l1tNNCaloTauProducer l1tNNCaloTauProducer.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/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" + +class l1tNNCaloTauProducer : public edm::stream::EDProducer<> { +public: + explicit l1tNNCaloTauProducer(const edm::ParameterSet&); + +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; + + 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; + + 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; + }; +}; + +/* +████████ ██ ██ ██████ ████████ █████ ██ ██ ███ ███ ██ ███ ██ █████ ████████ ██████ ██████ + ██ ██ ██ ██ ██ ██ ██ ██ ██ ████ ████ ██ ████ ██ ██ ██ ██ ██ ██ ██ ██ + ██ ███████ █████ ██ ███████ ██ ██ ██ ████ ██ ██ ██ ██ ██ ███████ ██ ██ ██ ██████ + ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ ██ + ██ ██ ██ ██████ ██ ██ ██ ███████ ██ ██ ██ ██ ████ ██ ██ ██ ██████ ██ ██ +*/ + +// ----Constructor and Destructor ----- +l1tNNCaloTauProducer::l1tNNCaloTauProducer(const edm::ParameterSet& iConfig) + : 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")), + + CNNmodel_CB_path(iConfig.getParameter("CNNmodel_CB_path")), + DNNident_CB_path(iConfig.getParameter("DNNident_CB_path")), + DNNcalib_CB_path(iConfig.getParameter("DNNcalib_CB_path")), + CNNmodel_CE_path(iConfig.getParameter("CNNmodel_CE_path")), + DNNident_CE_path(iConfig.getParameter("DNNident_CE_path")), + DNNcalib_CE_path(iConfig.getParameter("DNNcalib_CE_path")), + FeatScaler_CE_path(iConfig.getParameter("FeatScaler_CE_path")), + + 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")) { + // Create sessions for Tensorflow inferece + CNNmodel_CB = tensorflow::loadGraphDef(edm::FileInPath(CNNmodel_CB_path).fullPath()); + CNNmodel_CBsession = tensorflow::createSession(CNNmodel_CB); + + DNNident_CB = tensorflow::loadGraphDef(edm::FileInPath(DNNident_CB_path).fullPath()); + DNNident_CBsession = tensorflow::createSession(DNNident_CB); + + DNNcalib_CB = tensorflow::loadGraphDef(edm::FileInPath(DNNcalib_CB_path).fullPath()); + DNNcalib_CBsession = tensorflow::createSession(DNNcalib_CB); + + CNNmodel_CE = tensorflow::loadGraphDef(edm::FileInPath(CNNmodel_CE_path).fullPath()); + CNNmodel_CEsession = tensorflow::createSession(CNNmodel_CE); + + DNNident_CE = tensorflow::loadGraphDef(edm::FileInPath(DNNident_CE_path).fullPath()); + DNNident_CEsession = tensorflow::createSession(DNNident_CE); + + DNNcalib_CE = tensorflow::loadGraphDef(edm::FileInPath(DNNcalib_CE_path).fullPath()); + DNNcalib_CEsession = tensorflow::createSession(DNNcalib_CE); + + // Read features scaler + boost::property_tree::read_json(edm::FileInPath(FeatScaler_CE_path).fullPath(), FeatScaler_CE); + + // Initialize HGCAL BDTs + if (!VsPuId.method().empty()) { + VsPuId.prepareTMVA(); + } + + // Create produced outputs + produces>("L1NNCaloTauCollectionBXV"); + + // Settings output + std::cout << "EtaRestriction = " << EtaRestriction << " , CB_CE_split = " << CB_CE_split + << " , EtMinForSeeding = " << EtMinForSeeding << " , HcalTpEtMin = " << HcalEtMinForClustering + << " , EcalTpEtMin = " << EcalEtMinForClustering << std::endl; +} + +void l1tNNCaloTauProducer::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) { + std::cout << " ** 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 + } + + // Apply CNN model + tensorflow::NamedTensorList CNNmodel_CBinputList = {{"TowerClusterImage", TowerClusterImage_CB}, + {"TowerClusterPosition", TowerClusterPosition_CB}}; + std::vector CNNmodel_CBoutputs; + tensorflow::run( + 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( + DNNident_CBsession, DNN_CBinputsList, {"TauMinator_CB_ident/sigmoid_IDout/Sigmoid"}, &DNN_CBoutputsIdent); + + // Apply DNN for calibration + std::vector DNN_CBoutputsCalib; + tensorflow::run(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 + } + + // Apply CNN model + tensorflow::NamedTensorList CNNmodel_CEinputList = {{"TowerClusterImage", TowerClusterImage_CE}, + {"TowerClusterPosition", TowerClusterPosition_CE}, + {"AssociatedCl3dFeatures", Cl3dShapeFeatures_CE}}; + std::vector CNNmodel_CEoutputs; + tensorflow::run( + 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( + DNNident_CEsession, DNN_CEinputsList, {"TauMinator_CE_ident/sigmoid_IDout/Sigmoid"}, &DNN_CEoutputsIdent); + + // Apply DNN for calibration + std::vector DNN_CEoutputsCalib; + tensorflow::run(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 l1tNNCaloTauProducer::tower_dIPhi(int& iPhi_1, int& iPhi_2) 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 l1tNNCaloTauProducer::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 l1tNNCaloTauProducer::endcap_iphi(float& phi) const { + float phi_step = 0.0872664; + if (phi > 0) { + return floor(phi / phi_step) + 1; + } else { + return floor(phi / phi_step) + 73; + } +} + +int l1tNNCaloTauProducer::endcap_ieta(float& eta) const { + float eta_step = 0.0845; + return floor(abs(eta) / eta_step) * std::copysign(1, eta); +} + +float l1tNNCaloTauProducer::inputQuantizer(float inputF, float LSB, int nbits) { + return min(floor(inputF / LSB), float(pow(2, nbits) - 1)) * LSB; +} + +float l1tNNCaloTauProducer::inputScaler(float inputF, std::string feature) { + float mean = FeatScaler_CE.get_child(feature).get("mean"); + float std = FeatScaler_CE.get_child(feature).get("std"); + + return (inputF - mean) / std; +} + +DEFINE_FWK_MODULE(l1tNNCaloTauProducer); \ No newline at end of file diff --git a/L1Trigger/L1CaloTrigger/python/l1tNNCaloTauEmulator_cff.py b/L1Trigger/L1CaloTrigger/python/l1tNNCaloTauEmulator_cff.py new file mode 100644 index 0000000000000..32fc89da2fcd4 --- /dev/null +++ b/L1Trigger/L1CaloTrigger/python/l1tNNCaloTauEmulator_cff.py @@ -0,0 +1,5 @@ +import FWCore.ParameterSet.Config as cms + +from L1Trigger.L1CaloTrigger.l1tNNCaloTauEmulator_cfi import * + +l1tNNCaloTauEmulator_seq = cms.Sequence(l1tNNCaloTauEmulator) \ No newline at end of file diff --git a/L1Trigger/L1CaloTrigger/python/l1tNNCaloTauEmulator_cfi.py b/L1Trigger/L1CaloTrigger/python/l1tNNCaloTauEmulator_cfi.py new file mode 100644 index 0000000000000..8369345833f2a --- /dev/null +++ b/L1Trigger/L1CaloTrigger/python/l1tNNCaloTauEmulator_cfi.py @@ -0,0 +1,50 @@ +import FWCore.ParameterSet.Config as cms + +l1tNNCaloTauEmulator = cms.EDProducer("l1tNNCaloTauEmulator", + l1CaloTowers = cms.InputTag("l1tEGammaClusterEmuProducer","L1CaloTowerCollection",""), # uncalibrated towers (same input as L1TowerCalibrator) + hgcalTowers = cms.InputTag("l1tHGCalTowerProducer","HGCalTowerProcessor"), + + HgcalClusters = cms.InputTag("l1tHGCalBackEndLayer2Producer","HGCalBackendLayer2Processor3DClustering"), + preEmId = cms.string("hOverE < 0.3 && hOverE >= 0"), + VsPuId = cms.PSet( + isPUFilter = cms.bool(True), + preselection = cms.string(""), + method = cms.string("BDT"), # "" to be disabled, "BDT" to be enabled + variables = cms.VPSet( + cms.PSet(name = cms.string("eMax"), value = cms.string("eMax()")), + cms.PSet(name = cms.string("eMaxOverE"), value = cms.string("eMax()/energy()")), + cms.PSet(name = cms.string("sigmaPhiPhiTot"), value = cms.string("sigmaPhiPhiTot()")), + cms.PSet(name = cms.string("sigmaRRTot"), value = cms.string("sigmaRRTot()")), + cms.PSet(name = cms.string("triggerCells90percent"), value = cms.string("triggerCells90percent()")), + ), + weightsFile = cms.string("L1Trigger/Phase2L1ParticleFlow/data/hgcal_egID/Photon_Pion_vs_Neutrino_BDTweights_1116.xml.gz"), + wp = cms.string("-0.10") + ), + + EcalEtMinForClustering = cms.double(0.), + HcalEtMinForClustering = cms.double(0.), + EtMinForSeeding = cms.double(2.5), + EtaRestriction = cms.double(2.4), # the TauMinator algo can go up to 3.0, but 2.4 makes it directly comparable to "Square Calo Tau" algo + CB_CE_split = cms.double(1.55), + PuidThr = cms.double(-0.10), + + CNNmodel_CB_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/CNNmodel_CB.pb"), + DNNident_CB_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNident_CB.pb"), + DNNcalib_CB_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNcalib_CB.pb"), + + CNNmodel_CE_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/CNNmodel_CE.pb"), + DNNident_CE_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNident_CE.pb"), + DNNcalib_CE_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNcalib_CE.pb"), + FeatScaler_CE_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/Cl3dFeatScaler_CE.json"), + + # v22 + IdWp90_CB = cms.double(0.7060), + IdWp95_CB = cms.double(0.3432), + IdWp99_CB = cms.double(0.0337), + # v22 + IdWp90_CE = cms.double(0.5711), + IdWp95_CE = cms.double(0.2742), + IdWp99_CE = cms.double(0.0394), + + DEBUG = cms.bool(False) +) diff --git a/L1Trigger/L1CaloTrigger/python/l1tNNCaloTauProducer_cff.py b/L1Trigger/L1CaloTrigger/python/l1tNNCaloTauProducer_cff.py new file mode 100644 index 0000000000000..79aef7f3d15ef --- /dev/null +++ b/L1Trigger/L1CaloTrigger/python/l1tNNCaloTauProducer_cff.py @@ -0,0 +1,5 @@ +import FWCore.ParameterSet.Config as cms + +from L1Trigger.L1CaloTrigger.l1tNNCaloTauProducer_cfi import * + +l1tNNCaloTauProducer_seq = cms.Sequence(l1tNNCaloTauProducer) \ No newline at end of file diff --git a/L1Trigger/L1CaloTrigger/python/l1tNNCaloTauProducer_cfi.py b/L1Trigger/L1CaloTrigger/python/l1tNNCaloTauProducer_cfi.py new file mode 100644 index 0000000000000..0ff2d52098023 --- /dev/null +++ b/L1Trigger/L1CaloTrigger/python/l1tNNCaloTauProducer_cfi.py @@ -0,0 +1,49 @@ +import FWCore.ParameterSet.Config as cms + +l1tNNCaloTauProducer = cms.EDProducer("l1tNNCaloTauProducer", + l1CaloTowers = cms.InputTag("l1tEGammaClusterEmuProducer","L1CaloTowerCollection",""), # uncalibrated towers (same input as L1TowerCalibrator) + hgcalTowers = cms.InputTag("l1tHGCalTowerProducer","HGCalTowerProcessor"), + + HgcalClusters = cms.InputTag("l1tHGCalBackEndLayer2Producer","HGCalBackendLayer2Processor3DClustering"), + preEmId = cms.string("hOverE < 0.3 && hOverE >= 0"), + VsPuId = cms.PSet( + isPUFilter = cms.bool(True), + preselection = cms.string(""), + method = cms.string("BDT"), # "" to be disabled, "BDT" to be enabled + variables = cms.VPSet( + cms.PSet(name = cms.string("eMax"), value = cms.string("eMax()")), + cms.PSet(name = cms.string("eMaxOverE"), value = cms.string("eMax()/energy()")), + cms.PSet(name = cms.string("sigmaPhiPhiTot"), value = cms.string("sigmaPhiPhiTot()")), + cms.PSet(name = cms.string("sigmaRRTot"), value = cms.string("sigmaRRTot()")), + cms.PSet(name = cms.string("triggerCells90percent"), value = cms.string("triggerCells90percent()")), + ), + weightsFile = cms.string("L1Trigger/Phase2L1ParticleFlow/data/hgcal_egID/Photon_Pion_vs_Neutrino_BDTweights_1116.xml.gz"), + wp = cms.string("-0.10") + ), + + EcalEtMinForClustering = cms.double(0.), + HcalEtMinForClustering = cms.double(0.), + EtMinForSeeding = cms.double(2.5), + EtaRestriction = cms.double(2.4), # the TauMinator algo can go up to 3.0, but 2.4 makes it directly comparable to "Square Calo Tau" algo + CB_CE_split = cms.double(1.55), + + CNNmodel_CB_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/CNNmodel_CB.pb"), + DNNident_CB_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNident_CB.pb"), + DNNcalib_CB_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNcalib_CB.pb"), + + CNNmodel_CE_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/CNNmodel_CE.pb"), + DNNident_CE_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNident_CE.pb"), + DNNcalib_CE_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/v22/DNNcalib_CE.pb"), + FeatScaler_CE_path = cms.string("L1Trigger/L1CaloTrigger/data/Phase2_NNCaloTaus/Cl3dFeatScaler_CE.json"), + + # v22 + IdWp90_CB = cms.double(0.7060), + IdWp95_CB = cms.double(0.3432), + IdWp99_CB = cms.double(0.0337), + # v22 + IdWp90_CE = cms.double(0.5711), + IdWp95_CE = cms.double(0.2742), + IdWp99_CE = cms.double(0.0394), + + DEBUG = cms.bool(False) +)