diff --git a/L1Trigger/L1CaloTrigger/BuildFile.xml b/L1Trigger/L1CaloTrigger/BuildFile.xml
new file mode 100644
index 0000000000000..ea6148c865666
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/BuildFile.xml
@@ -0,0 +1,7 @@
+
+
+
+
+
+
+
diff --git a/L1Trigger/L1CaloTrigger/interface/L1EGammaEECalibrator.h b/L1Trigger/L1CaloTrigger/interface/L1EGammaEECalibrator.h
new file mode 100644
index 0000000000000..51a8efe024e57
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/interface/L1EGammaEECalibrator.h
@@ -0,0 +1,25 @@
+#ifndef L1Trigger_L1CaloTrigger_L1EGammaEECalibrator_h
+#define L1Trigger_L1CaloTrigger_L1EGammaEECalibrator_h
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include
+#include
+#include
+
+class L1EGammaEECalibrator {
+public:
+ explicit L1EGammaEECalibrator(const edm::ParameterSet&);
+
+ float calibrationFactor(const float& pt, const float& eta) const;
+
+private:
+ int etaBin(float eta) const { return bin(eta_bins, std::abs(eta)); }
+ int ptBin(float pt) const { return bin(pt_bins, pt); }
+ int bin(const std::set& container, float value) const;
+
+ std::set eta_bins;
+ std::set pt_bins;
+ std::vector calib_factors;
+};
+
+#endif
diff --git a/L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h b/L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h
new file mode 100644
index 0000000000000..feefed6b752fe
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h
@@ -0,0 +1,26 @@
+#ifndef L1Trigger_L1CaloTrigger_ParametricCalibration_h
+#define L1Trigger_L1CaloTrigger_ParametricCalibration_h
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/Exception.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include
+#include
+#include
+
+namespace l1tp2 {
+ class ParametricCalibration {
+ public:
+ ParametricCalibration() {}
+ ParametricCalibration(const edm::ParameterSet& cpset);
+ static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+ float operator()(const float pt, const float abseta) const;
+
+ protected:
+ std::vector etas, pts, scales;
+ };
+
+}; // namespace l1tp2
+
+#endif
diff --git a/L1Trigger/L1CaloTrigger/plugins/BuildFile.xml b/L1Trigger/L1CaloTrigger/plugins/BuildFile.xml
new file mode 100644
index 0000000000000..041ad226b7800
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/plugins/BuildFile.xml
@@ -0,0 +1,32 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc b/L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc
new file mode 100644
index 0000000000000..2781e7b15d7bd
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc
@@ -0,0 +1,1204 @@
+// -*- C++ -*-
+//
+// Package: L1CaloTrigger
+// Class: L1EGammaCrystalsEmulatorProducer
+//
+/**\class L1EGammaCrystalsEmulatorProducer L1EGammaCrystalsEmulatorProducer.cc L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc
+Description: Produces crystal clusters using crystal-level information and hardware constraints
+Implementation:
+[Notes on implementation]
+*/
+//
+// Original Author: Cecile Caillol
+// Created: Tue Aug 10 2018
+//
+// Redesign, Calibration: Vladimir Rekovic
+// Date: Tue May 12 2020
+//
+// $Id$
+//
+//
+
+// system include files
+#include
+#include
+#include
+#include
+
+// user include files
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h"
+#include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h"
+#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
+#include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
+#include "Geometry/HcalTowerAlgo/interface/HcalTrigTowerGeometry.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
+#include "DataFormats/HcalDetId/interface/HcalDetId.h"
+
+// ECAL TPs
+#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
+
+// HCAL TPs
+#include "DataFormats/HcalDigi/interface/HcalTriggerPrimitiveDigi.h"
+
+// Output tower collection
+#include "DataFormats/L1TCalorimeterPhase2/interface/CaloCrystalCluster.h"
+#include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h"
+#include "DataFormats/L1Trigger/interface/EGamma.h"
+
+#include "L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h"
+#include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+static constexpr bool do_brem = true;
+
+static constexpr int n_eta_bins = 2;
+static constexpr int n_borders_phi = 18;
+static constexpr int n_borders_eta = 18;
+static constexpr int n_clusters_max = 5;
+static constexpr int n_clusters_link = 3;
+static constexpr int n_clusters_4link = 4 * 3;
+static constexpr int n_crystals_towerEta = 5;
+static constexpr int n_crystals_towerPhi = 5;
+static constexpr int n_crystals_3towers = 3 * 5;
+static constexpr int n_towers_per_link = 17;
+static constexpr int n_clusters_per_link = 2;
+static constexpr int n_towers_Eta = 34;
+static constexpr int n_towers_Phi = 72;
+static constexpr int n_towers_halfPhi = 36;
+static constexpr int n_links_card = 4;
+static constexpr int n_links_GCTcard = 48;
+static constexpr int n_GCTcards = 3;
+static constexpr float ECAL_eta_range = 1.4841;
+static constexpr float half_crystal_size = 0.00873;
+static constexpr float slideIsoPtThreshold = 80;
+static constexpr float a0_80 = 0.85, a1_80 = 0.0080, a0 = 0.21; // passes_iso
+static constexpr float b0 = 0.38, b1 = 1.9, b2 = 0.05; //passes_looseTkiso
+static constexpr float c0_ss = 0.94, c1_ss = 0.052, c2_ss = 0.044; //passes_ss
+static constexpr float d0 = 0.96, d1 = 0.0003; //passes_photon
+static constexpr float e0_looseTkss = 0.944, e1_looseTkss = 0.65, e2_looseTkss = 0.4; //passes_looseTkss
+static constexpr float cut_500_MeV = 0.5;
+
+// absolue IDs range from 0-33
+// 0-16 are iEta -17 to -1
+// 17 to 33 are iEta 1 to 17
+static constexpr int toweriEta_fromAbsoluteID_shift = 16;
+
+// absolue IDs range from 0-71.
+// To align with detector tower IDs (1 - n_towers_Phi)
+// shift all indices by 37 and loop over after 72
+static constexpr int toweriPhi_fromAbsoluteID_shift_lowerHalf = 37;
+static constexpr int toweriPhi_fromAbsoluteID_shift_upperHalf = 35;
+
+float getEta_fromL2LinkCardTowerCrystal(int link, int card, int tower, int crystal) {
+ int etaID = n_crystals_towerEta * (n_towers_per_link * ((link / n_links_card) % 2) + (tower % n_towers_per_link)) +
+ crystal % n_crystals_towerEta;
+ float size_cell = 2 * ECAL_eta_range / (n_crystals_towerEta * n_towers_Eta);
+ return etaID * size_cell - ECAL_eta_range + half_crystal_size;
+}
+
+float getPhi_fromL2LinkCardTowerCrystal(int link, int card, int tower, int crystal) {
+ int phiID = n_crystals_towerPhi * ((card * 24) + (n_links_card * (link / 8)) + (tower / n_towers_per_link)) +
+ crystal / n_crystals_towerPhi;
+ float size_cell = 2 * M_PI / (n_crystals_towerPhi * n_towers_Phi);
+ return phiID * size_cell - M_PI + half_crystal_size;
+}
+
+int getCrystal_etaID(float eta) {
+ float size_cell = 2 * ECAL_eta_range / (n_crystals_towerEta * n_towers_Eta);
+ int etaID = int((eta + ECAL_eta_range) / size_cell);
+ return etaID;
+}
+
+int convert_L2toL1_link(int link) { return link % n_links_card; }
+
+int convert_L2toL1_tower(int tower) { return tower; }
+
+int convert_L2toL1_card(int card, int link) { return card * n_clusters_4link + link / n_links_card; }
+
+int getCrystal_phiID(float phi) {
+ float size_cell = 2 * M_PI / (n_crystals_towerPhi * n_towers_Phi);
+ int phiID = int((phi + M_PI) / size_cell);
+ return phiID;
+}
+
+int getTower_absoluteEtaID(float eta) {
+ float size_cell = 2 * ECAL_eta_range / n_towers_Eta;
+ int etaID = int((eta + ECAL_eta_range) / size_cell);
+ return etaID;
+}
+
+int getTower_absolutePhiID(float phi) {
+ float size_cell = 2 * M_PI / n_towers_Phi;
+ int phiID = int((phi + M_PI) / size_cell);
+ return phiID;
+}
+
+int getToweriEta_fromAbsoluteID(int id) {
+ if (id < n_towers_per_link)
+ return id - n_towers_per_link;
+ else
+ return id - toweriEta_fromAbsoluteID_shift;
+}
+
+int getToweriPhi_fromAbsoluteID(int id) {
+ if (id < n_towers_Phi / 2)
+ return id + toweriPhi_fromAbsoluteID_shift_lowerHalf;
+ else
+ return id - toweriPhi_fromAbsoluteID_shift_upperHalf;
+}
+
+float getTowerEta_fromAbsoluteID(int id) {
+ float size_cell = 2 * ECAL_eta_range / n_towers_Eta;
+ float eta = (id * size_cell) - ECAL_eta_range + 0.5 * size_cell;
+ return eta;
+}
+
+float getTowerPhi_fromAbsoluteID(int id) {
+ float size_cell = 2 * M_PI / n_towers_Phi;
+ float phi = (id * size_cell) - M_PI + 0.5 * size_cell;
+ return phi;
+}
+
+int getCrystalIDInTower(int etaID, int phiID) {
+ return int(n_crystals_towerPhi * (phiID % n_crystals_towerPhi) + (etaID % n_crystals_towerEta));
+}
+
+int get_towerEta_fromCardTowerInCard(int card, int towerincard) {
+ return n_towers_per_link * (card % 2) + towerincard % n_towers_per_link;
+}
+
+int get_towerPhi_fromCardTowerInCard(int card, int towerincard) {
+ return 4 * (card / 2) + towerincard / n_towers_per_link;
+}
+
+int get_towerEta_fromCardLinkTower(int card, int link, int tower) { return n_towers_per_link * (card % 2) + tower; }
+
+int get_towerPhi_fromCardLinkTower(int card, int link, int tower) { return 4 * (card / 2) + link; }
+
+int getTowerID(int etaID, int phiID) {
+ return int(n_towers_per_link * ((phiID / n_crystals_towerPhi) % 4) +
+ (etaID / n_crystals_towerEta) % n_towers_per_link);
+}
+
+int getTower_phiID(int cluster_phiID) { // Tower ID in card given crystal ID in total detector
+ return int((cluster_phiID / n_crystals_towerPhi) % 4);
+}
+
+int getTower_etaID(int cluster_etaID) { // Tower ID in card given crystal ID in total detector
+ return int((cluster_etaID / n_crystals_towerEta) % n_towers_per_link);
+}
+
+int getEtaMax_card(int card) {
+ int etamax = 0;
+ if (card % 2 == 0)
+ etamax = n_towers_per_link * n_crystals_towerEta - 1; // First eta half. 5 crystals in eta in 1 tower.
+ else
+ etamax = n_towers_Eta * n_crystals_towerEta - 1;
+ return etamax;
+}
+
+int getEtaMin_card(int card) {
+ int etamin = 0;
+ if (card % 2 == 0)
+ etamin = 0 * n_crystals_towerEta; // First eta half. 5 crystals in eta in 1 tower.
+ else
+ etamin = n_towers_per_link * n_crystals_towerEta;
+ return etamin;
+}
+
+int getPhiMax_card(int card) {
+ int phimax = ((card / 2) + 1) * 4 * n_crystals_towerPhi - 1;
+ return phimax;
+}
+
+int getPhiMin_card(int card) {
+ int phimin = (card / 2) * 4 * n_crystals_towerPhi;
+ return phimin;
+}
+
+class L1EGCrystalClusterEmulatorProducer : public edm::stream::EDProducer<> {
+public:
+ explicit L1EGCrystalClusterEmulatorProducer(const edm::ParameterSet&);
+ ~L1EGCrystalClusterEmulatorProducer() override;
+
+private:
+ void produce(edm::Event&, const edm::EventSetup&) override;
+ bool passes_ss(float pt, float ss);
+ bool passes_photon(float pt, float pss);
+ bool passes_iso(float pt, float iso);
+ bool passes_looseTkss(float pt, float ss);
+ bool passes_looseTkiso(float pt, float iso);
+ float get_calibrate(float uncorr);
+
+ edm::EDGetTokenT ecalTPEBToken_;
+ edm::EDGetTokenT > hcalTPToken_;
+ edm::ESHandle decoder_;
+
+ l1tp2::ParametricCalibration calib_;
+
+ edm::ESHandle caloGeometry_;
+ const CaloSubdetectorGeometry* ebGeometry;
+ const CaloSubdetectorGeometry* hbGeometry;
+ edm::ESHandle hbTopology;
+ const HcalTopology* hcTopology_;
+
+ struct mycluster {
+ float c2x2_;
+ float c2x5_;
+ float c5x5_;
+ int cshowershape_;
+ int cshowershapeloosetk_;
+ float cvalueshowershape_;
+ int cphotonshowershape_;
+ float cpt; // ECAL pt
+ int cbrem_; // if brem corrections were applied
+ float cWeightedEta_;
+ float cWeightedPhi_;
+ float ciso_; // pt of cluster divided by 7x7 ECAL towers
+ float chovere_; // 5x5 HCAL towers divided by the ECAL cluster pt
+ float craweta_; // coordinates between -1.44 and 1.44
+ float crawphi_; // coordinates between -pi and pi
+ float chcal_; // 5x5 HCAL towers
+ float ceta_; // eta ID in the whole detector (between 0 and 5*34-1)
+ float cphi_; // phi ID in the whole detector (between 0 and 5*72-1)
+ int ccrystalid_; // crystal ID inside tower (between 0 and 24)
+ int cinsidecrystalid_;
+ int ctowerid_; // tower ID inside card (between 0 and 4*n_towers_per_link-1)
+ };
+
+ class SimpleCaloHit {
+ private:
+ float pt_ = 0;
+ float energy_ = 0.;
+ bool isEndcapHit_ = false; // If using endcap, we won't be using integer crystal indices
+ bool stale_ = false; // Hits become stale once used in clustering algorithm to prevent overlap in clusters
+ bool used_ = false;
+ GlobalVector position_; // As opposed to GlobalPoint, so we can add them (for weighted average)
+ HcalDetId id_hcal_;
+ EBDetId id_;
+
+ public:
+ // tool functions
+ inline void setPt() { pt_ = (position_.mag2() > 0) ? energy_ * sin(position_.theta()) : 0; };
+ inline void setEnergy(float et) { energy_ = et / sin(position_.theta()); };
+ inline void setIsEndcapHit(bool isEC) { isEndcapHit_ = isEC; };
+ inline void setUsed(bool isUsed) { used_ = isUsed; };
+ inline void setPosition(const GlobalVector& pos) { position_ = pos; };
+ inline void setIdHcal(const HcalDetId& idhcal) { id_hcal_ = idhcal; };
+ inline void setId(const EBDetId& id) { id_ = id; };
+
+ inline float pt() const { return pt_; };
+ inline float energy() const { return energy_; };
+ inline bool isEndcapHit() const { return isEndcapHit_; };
+ inline bool used() const { return used_; };
+ inline const GlobalVector& position() const { return position_; };
+ inline const EBDetId& id() const { return id_; };
+
+ inline float deta(SimpleCaloHit& other) const { return position_.eta() - other.position().eta(); };
+ int dieta(SimpleCaloHit& other) const {
+ if (isEndcapHit_ || other.isEndcapHit())
+ return 9999; // We shouldn't compare integer indices in endcap, the map is not linear
+ if (id_.ieta() * other.id().ieta() > 0)
+ return id_.ieta() - other.id().ieta();
+ return id_.ieta() - other.id().ieta() - 1;
+ };
+ inline float dphi(SimpleCaloHit& other) const {
+ return reco::deltaPhi(static_cast(position_.phi()), static_cast(other.position().phi()));
+ };
+ int diphi(SimpleCaloHit& other) const {
+ if (isEndcapHit_ || other.isEndcapHit())
+ return 9999; // We shouldn't compare integer indices in endcap, the map is not linear
+ // Logic from EBDetId::distancePhi() without the abs()
+ static constexpr int PI = 180;
+ int result = id().iphi() - other.id().iphi();
+ while (result > PI)
+ result -= 2 * PI;
+ while (result <= -PI)
+ result += 2 * PI;
+ return result;
+ };
+ float distanceTo(SimpleCaloHit& other) const {
+ // Treat position as a point, measure 3D distance
+ // This is used for endcap hits, where we don't have a rectangular mapping
+ return (position() - other.position()).mag();
+ };
+ bool operator==(SimpleCaloHit& other) const {
+ return (id_ == other.id() && position() == other.position() && energy_ == other.energy() &&
+ isEndcapHit_ == other.isEndcapHit());
+ };
+ };
+};
+
+L1EGCrystalClusterEmulatorProducer::L1EGCrystalClusterEmulatorProducer(const edm::ParameterSet& iConfig)
+ : ecalTPEBToken_(consumes(iConfig.getParameter("ecalTPEB"))),
+ hcalTPToken_(
+ consumes >(iConfig.getParameter("hcalTP"))),
+ calib_(iConfig.getParameter("calib")) {
+ produces();
+ produces >();
+ produces();
+}
+
+L1EGCrystalClusterEmulatorProducer::~L1EGCrystalClusterEmulatorProducer() {}
+
+void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+ using namespace edm;
+
+ edm::Handle pcalohits;
+ iEvent.getByToken(ecalTPEBToken_, pcalohits);
+
+ iSetup.get().get(caloGeometry_);
+ ebGeometry = caloGeometry_->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
+ hbGeometry = caloGeometry_->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
+ iSetup.get().get(hbTopology);
+ hcTopology_ = hbTopology.product();
+ HcalTrigTowerGeometry theTrigTowerGeometry(hcTopology_);
+ iSetup.get().get(decoder_);
+
+ //****************************************************************
+ //******************* Get all the hits ***************************
+ //****************************************************************
+
+ // Get all the ECAL hits
+ iEvent.getByToken(ecalTPEBToken_, pcalohits);
+ std::vector ecalhits;
+
+ for (const auto& hit : *pcalohits.product()) {
+ if (hit.encodedEt() > 0) // hit.encodedEt() returns an int corresponding to 2x the crystal Et
+ {
+ // Et is 10 bit, by keeping the ADC saturation Et at 120 GeV it means that you have to divide by 8
+ float et = hit.encodedEt() / 8.;
+ if (et < cut_500_MeV)
+ continue; // keep the 500 MeV ET Cut
+
+ auto cell = ebGeometry->getGeometry(hit.id());
+
+ SimpleCaloHit ehit;
+ ehit.setId(hit.id());
+ ehit.setPosition(GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z()));
+ ehit.setEnergy(et);
+ ehit.setPt();
+ ecalhits.push_back(ehit);
+ }
+ }
+
+ // Get all the HCAL hits
+ std::vector hcalhits;
+ edm::Handle > hbhecoll;
+ iEvent.getByToken(hcalTPToken_, hbhecoll);
+ for (const auto& hit : *hbhecoll.product()) {
+ float et = decoder_->hcaletValue(hit.id(), hit.t0());
+ if (et <= 0)
+ continue;
+ if (!(hcTopology_->validHT(hit.id()))) {
+ LogError("L1EGCrystalClusterEmulatorProducer")
+ << " -- Hcal hit DetID not present in HCAL Geom: " << hit.id() << std::endl;
+ throw cms::Exception("L1EGCrystalClusterEmulatorProducer");
+ continue;
+ }
+ const std::vector& hcId = theTrigTowerGeometry.detIds(hit.id());
+ if (hcId.empty()) {
+ LogError("L1EGCrystalClusterEmulatorProducer")
+ << "Cannot find any HCalDetId corresponding to " << hit.id() << std::endl;
+ throw cms::Exception("L1EGCrystalClusterEmulatorProducer");
+ continue;
+ }
+ if (hcId[0].subdetId() > 1)
+ continue;
+ GlobalVector hcal_tp_position = GlobalVector(0., 0., 0.);
+ for (const auto& hcId_i : hcId) {
+ if (hcId_i.subdetId() > 1)
+ continue;
+ auto cell = hbGeometry->getGeometry(hcId_i);
+ if (cell == nullptr)
+ continue;
+ GlobalVector tmpVector = GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z());
+ hcal_tp_position = tmpVector;
+ break;
+ }
+ SimpleCaloHit hhit;
+ hhit.setId(hit.id());
+ hhit.setIdHcal(hit.id());
+ hhit.setPosition(hcal_tp_position);
+ hhit.setEnergy(et);
+ hhit.setPt();
+ hcalhits.push_back(hhit);
+ }
+
+ //*******************************************************************
+ //********************** Do layer 1 *********************************
+ //*******************************************************************
+
+ // Definition of L1 outputs
+ // 36 L1 cards send each 4 links with 17 towers
+ float ECAL_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
+ float HCAL_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
+ int iEta_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
+ int iPhi_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
+ // 36 L1 cards send each 4 links with 3 clusters
+ float energy_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
+ // 36 L1 cards send each 4 links with 3 clusters
+ int brem_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
+ int towerID_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
+ int crystalID_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
+ int showerShape_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
+ int showerShapeLooseTk_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
+ int photonShowerShape_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
+
+ for (int ii = 0; ii < n_links_card; ++ii) {
+ for (int jj = 0; jj < n_towers_per_link; ++jj) {
+ for (int ll = 0; ll < n_towers_halfPhi; ++ll) {
+ ECAL_tower_L1Card[ii][jj][ll] = 0;
+ HCAL_tower_L1Card[ii][jj][ll] = 0;
+ iPhi_tower_L1Card[ii][jj][ll] = -999;
+ iEta_tower_L1Card[ii][jj][ll] = -999;
+ }
+ }
+ }
+ for (int ii = 0; ii < n_links_card; ++ii) {
+ for (int jj = 0; jj < n_clusters_link; ++jj) {
+ for (int ll = 0; ll < n_towers_halfPhi; ++ll) {
+ energy_cluster_L1Card[ii][jj][ll] = 0;
+ brem_cluster_L1Card[ii][jj][ll] = 0;
+ towerID_cluster_L1Card[ii][jj][ll] = 0;
+ crystalID_cluster_L1Card[ii][jj][ll] = 0;
+ }
+ }
+ }
+
+ // There is one list of clusters per card. We take the 12 highest pt per card
+ vector cluster_list[n_towers_halfPhi];
+ // After merging the clusters in different regions of a single L1 card
+ vector cluster_list_merged[n_towers_halfPhi];
+
+ for (int cc = 0; cc < n_towers_halfPhi; ++cc) { // Loop over 36 L1 cards
+ // Loop over 3x4 etaxphi regions to search for max 5 clusters
+ for (int nregion = 0; nregion <= n_clusters_max; ++nregion) {
+ int nclusters = 0;
+ bool build_cluster = true;
+
+ // Continue until 5 clusters have been built or there is no cluster left
+ while (nclusters < n_clusters_max && build_cluster) {
+ build_cluster = false;
+ SimpleCaloHit centerhit;
+
+ for (const auto& hit : ecalhits) {
+ if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
+ getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
+ getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
+ getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) &&
+ // Check that the hit is in the good card
+ getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) &&
+ getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion &&
+ !hit.used() && hit.pt() >= 1.0 && hit.pt() > centerhit.pt()) // 3 towers x 5 crystals
+ {
+ // Highest hit in good region with pt>1 and not used in any other cluster
+ centerhit = hit;
+ build_cluster = true;
+ }
+ }
+ if (build_cluster)
+ nclusters++;
+
+ // Use only the 5 most energetic clusters
+ if (build_cluster && nclusters > 0 && nclusters <= n_clusters_max) {
+ mycluster mc1;
+ mc1.cpt = 0.0;
+ mc1.cWeightedEta_ = 0.0;
+ mc1.cWeightedPhi_ = 0.0;
+ float leftlobe = 0;
+ float rightlobe = 0;
+ float e5x5 = 0;
+ float n5x5 = 0;
+ float e2x5_1 = 0;
+ float n2x5_1 = 0;
+ float e2x5_2 = 0;
+ float n2x5_2 = 0;
+ float e2x2_1 = 0;
+ float n2x2_1 = 0;
+ float e2x2_2 = 0;
+ float n2x2_2 = 0;
+ float e2x2_3 = 0;
+ float n2x2_3 = 0;
+ float e2x2_4 = 0;
+ float n2x2_4 = 0;
+ for (auto& hit : ecalhits) {
+ if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
+ getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
+ getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
+ getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && hit.pt() > 0 &&
+ getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) &&
+ getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion) {
+ if (abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) > 2 && hit.diphi(centerhit) <= 7) {
+ rightlobe += hit.pt();
+ }
+ if (abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) < -2 && hit.diphi(centerhit) >= -7) {
+ leftlobe += hit.pt();
+ }
+ if (abs(hit.dieta(centerhit)) <= 2 && abs(hit.diphi(centerhit)) <= 2) {
+ e5x5 += hit.energy();
+ n5x5++;
+ }
+ if ((hit.dieta(centerhit) == 1 or hit.dieta(centerhit) == 0) &&
+ (hit.diphi(centerhit) == 1 or hit.diphi(centerhit) == 0)) {
+ e2x2_1 += hit.energy();
+ n2x2_1++;
+ }
+ if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) &&
+ (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == 1)) {
+ e2x2_2 += hit.energy();
+ n2x2_2++;
+ }
+ if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == 1) &&
+ (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == -1)) {
+ e2x2_3 += hit.energy();
+ n2x2_3++;
+ }
+ if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) &&
+ (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == -1)) {
+ e2x2_4 += hit.energy();
+ n2x2_4++;
+ }
+ if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == 1) && abs(hit.diphi(centerhit)) <= 2) {
+ e2x5_1 += hit.energy();
+ n2x5_1++;
+ }
+ if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) && abs(hit.diphi(centerhit)) <= 2) {
+ e2x5_2 += hit.energy();
+ n2x5_2++;
+ }
+ }
+ if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
+ getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
+ getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
+ getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && !hit.used() && hit.pt() > 0 &&
+ abs(hit.dieta(centerhit)) <= 1 && abs(hit.diphi(centerhit)) <= 2 &&
+ getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) &&
+ getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion) {
+ // clusters 3x5 in etaxphi using only the hits in the corresponding card and in the corresponding 3x4 region
+ hit.setUsed(true);
+ mc1.cpt += hit.pt();
+ mc1.cWeightedEta_ += float(hit.pt()) * float(hit.position().eta());
+ mc1.cWeightedPhi_ = mc1.cWeightedPhi_ + (float(hit.pt()) * float(hit.position().phi()));
+ }
+ }
+ if (do_brem && (rightlobe > 0.10 * mc1.cpt or leftlobe > 0.10 * mc1.cpt)) {
+ for (auto& hit : ecalhits) {
+ if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
+ getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
+ getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
+ getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && hit.pt() > 0 &&
+ getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) &&
+ getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion &&
+ !hit.used()) {
+ if (rightlobe > 0.10 * mc1.cpt && (leftlobe < 0.10 * mc1.cpt or rightlobe > leftlobe) &&
+ abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) > 2 && hit.diphi(centerhit) <= 7) {
+ mc1.cpt += hit.pt();
+ hit.setUsed(true);
+ mc1.cbrem_ = 1;
+ }
+ if (leftlobe > 0.10 * mc1.cpt && (rightlobe < 0.10 * mc1.cpt or leftlobe >= rightlobe) &&
+ abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) < -2 && hit.diphi(centerhit) >= -7) {
+ mc1.cpt += hit.pt();
+ hit.setUsed(true);
+ mc1.cbrem_ = 1;
+ }
+ }
+ }
+ }
+ mc1.c5x5_ = e5x5;
+ mc1.c2x5_ = max(e2x5_1, e2x5_2);
+ mc1.c2x2_ = e2x2_1;
+ if (e2x2_2 > mc1.c2x2_)
+ mc1.c2x2_ = e2x2_2;
+ if (e2x2_3 > mc1.c2x2_)
+ mc1.c2x2_ = e2x2_3;
+ if (e2x2_4 > mc1.c2x2_)
+ mc1.c2x2_ = e2x2_4;
+ mc1.cWeightedEta_ = mc1.cWeightedEta_ / mc1.cpt;
+ mc1.cWeightedPhi_ = mc1.cWeightedPhi_ / mc1.cpt;
+ mc1.ceta_ = getCrystal_etaID(centerhit.position().eta());
+ mc1.cphi_ = getCrystal_phiID(centerhit.position().phi());
+ mc1.crawphi_ = centerhit.position().phi();
+ mc1.craweta_ = centerhit.position().eta();
+ cluster_list[cc].push_back(mc1);
+ } // End if 5 clusters per region
+ } // End while to find the 5 clusters
+ } // End loop over regions to search for clusters
+ std::sort(begin(cluster_list[cc]), end(cluster_list[cc]), [](mycluster a, mycluster b) { return a.cpt > b.cpt; });
+
+ // Merge clusters from different regions
+ for (unsigned int jj = 0; jj < unsigned(cluster_list[cc].size()); ++jj) {
+ for (unsigned int kk = jj + 1; kk < unsigned(cluster_list[cc].size()); ++kk) {
+ if (std::abs(cluster_list[cc][jj].ceta_ - cluster_list[cc][kk].ceta_) < 2 &&
+ std::abs(cluster_list[cc][jj].cphi_ - cluster_list[cc][kk].cphi_) < 2) { //Diagonale + exact neighbors
+ if (cluster_list[cc][kk].cpt > cluster_list[cc][jj].cpt) {
+ cluster_list[cc][kk].cpt += cluster_list[cc][jj].cpt;
+ cluster_list[cc][kk].c5x5_ += cluster_list[cc][jj].c5x5_;
+ cluster_list[cc][kk].c2x5_ += cluster_list[cc][jj].c2x5_;
+ cluster_list[cc][jj].cpt = 0;
+ cluster_list[cc][jj].c5x5_ = 0;
+ cluster_list[cc][jj].c2x5_ = 0;
+ cluster_list[cc][jj].c2x2_ = 0;
+ } else {
+ cluster_list[cc][jj].cpt += cluster_list[cc][kk].cpt;
+ cluster_list[cc][jj].c5x5_ += cluster_list[cc][kk].c5x5_;
+ cluster_list[cc][jj].c2x5_ += cluster_list[cc][kk].c2x5_;
+ cluster_list[cc][kk].cpt = 0;
+ cluster_list[cc][kk].c2x2_ = 0;
+ cluster_list[cc][kk].c2x5_ = 0;
+ cluster_list[cc][kk].c5x5_ = 0;
+ }
+ }
+ }
+ if (cluster_list[cc][jj].cpt > 0) {
+ cluster_list[cc][jj].cpt =
+ cluster_list[cc][jj].cpt *
+ calib_(cluster_list[cc][jj].cpt,
+ std::abs(cluster_list[cc][jj].craweta_)); //Mark's calibration as a function of eta and pt
+ cluster_list_merged[cc].push_back(cluster_list[cc][jj]);
+ }
+ }
+ std::sort(begin(cluster_list_merged[cc]), end(cluster_list_merged[cc]), [](mycluster a, mycluster b) {
+ return a.cpt > b.cpt;
+ });
+
+ // Fill cluster information in the arrays. We keep max 12 clusters (distributed between 4 links)
+ for (unsigned int jj = 0; jj < unsigned(cluster_list_merged[cc].size()) && jj < n_clusters_4link; ++jj) {
+ crystalID_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] =
+ getCrystalIDInTower(cluster_list_merged[cc][jj].ceta_, cluster_list_merged[cc][jj].cphi_);
+ towerID_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] =
+ getTowerID(cluster_list_merged[cc][jj].ceta_, cluster_list_merged[cc][jj].cphi_);
+ energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = cluster_list_merged[cc][jj].cpt;
+ brem_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = cluster_list_merged[cc][jj].cbrem_;
+ if (passes_ss(cluster_list_merged[cc][jj].cpt,
+ cluster_list_merged[cc][jj].c2x5_ / cluster_list_merged[cc][jj].c5x5_))
+ showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1;
+ else
+ showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0;
+ if (passes_looseTkss(cluster_list_merged[cc][jj].cpt,
+ cluster_list_merged[cc][jj].c2x5_ / cluster_list_merged[cc][jj].c5x5_))
+ showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1;
+ else
+ showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0;
+ if (passes_photon(cluster_list_merged[cc][jj].cpt,
+ cluster_list_merged[cc][jj].c2x2_ / cluster_list_merged[cc][jj].c2x5_))
+ photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1;
+ else
+ photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0;
+ }
+
+ // Loop over calo ecal hits to get the ECAL towers. Take only hits that have not been used to make clusters
+ for (const auto& hit : ecalhits) {
+ if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
+ getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
+ getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
+ getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) &&
+ !hit.used()) { // Take all the hits inside the card that have not been used yet
+ for (int jj = 0; jj < n_links_card; ++jj) { // loop over 4 links per card
+ if ((getCrystal_phiID(hit.position().phi()) / n_crystals_towerPhi) % 4 == jj) { // Go to ID tower modulo 4
+ for (int ii = 0; ii < n_towers_per_link; ++ii) {
+ //Apply Mark's calibration at the same time (row of the lowest pT, as a function of eta)
+ if ((getCrystal_etaID(hit.position().eta()) / n_crystals_towerEta) % n_towers_per_link == ii) {
+ ECAL_tower_L1Card[jj][ii][cc] += hit.pt() * calib_(0, std::abs(hit.position().eta()));
+ iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(hit.position().eta()); //hit.id().ieta();
+ iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(hit.position().phi()); //hit.id().iphi();
+ }
+ } // end of loop over eta towers
+ }
+ } // end of loop over phi links
+ // Make sure towers with 0 ET are initialized with proper iEta, iPhi coordinates
+ static constexpr float tower_width = 0.0873;
+ for (int jj = 0; jj < n_links_card; ++jj) {
+ for (int ii = 0; ii < n_towers_per_link; ++ii) {
+ float phi = getPhiMin_card(cc) * tower_width / n_crystals_towerPhi - M_PI + (jj + 0.5) * tower_width;
+ float eta = getEtaMin_card(cc) * tower_width / n_crystals_towerEta - n_towers_per_link * tower_width +
+ (ii + 0.5) * tower_width;
+ iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(eta);
+ iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(phi);
+ }
+ }
+
+ } // end of check if inside card
+ } // end of loop over hits to build towers
+
+ // Loop over hcal hits to get the HCAL towers.
+ for (const auto& hit : hcalhits) {
+ if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
+ getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
+ getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
+ getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && hit.pt() > 0) {
+ for (int jj = 0; jj < n_links_card; ++jj) {
+ if ((getCrystal_phiID(hit.position().phi()) / n_crystals_towerPhi) % n_links_card == jj) {
+ for (int ii = 0; ii < n_towers_per_link; ++ii) {
+ if ((getCrystal_etaID(hit.position().eta()) / n_crystals_towerEta) % n_towers_per_link == ii) {
+ HCAL_tower_L1Card[jj][ii][cc] += hit.pt();
+ iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(hit.position().eta()); //hit.id().ieta();
+ iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(hit.position().phi()); //hit.id().iphi();
+ }
+ } // end of loop over eta towers
+ }
+ } // end of loop over phi links
+ } // end of check if inside card
+ } // end of loop over hits to build towers
+
+ // Give back energy of not used clusters to the towers (if there are more than 12 clusters)
+ for (unsigned int kk = n_clusters_4link; kk < cluster_list_merged[cc].size(); ++kk) {
+ if (cluster_list_merged[cc][kk].cpt > 0) {
+ ECAL_tower_L1Card[getTower_phiID(cluster_list_merged[cc][kk].cphi_)]
+ [getTower_etaID(cluster_list_merged[cc][kk].ceta_)][cc] += cluster_list_merged[cc][kk].cpt;
+ }
+ }
+ } //End of loop over cards
+
+ //*********************************************************
+ //******************** Do layer 2 *************************
+ //*********************************************************
+
+ // Definition of L2 outputs
+ float HCAL_tower_L2Card[n_links_GCTcard][n_towers_per_link]
+ [n_GCTcards]; // 3 L2 cards send each 48 links with 17 towers
+ float ECAL_tower_L2Card[n_links_GCTcard][n_towers_per_link][n_GCTcards];
+ int iEta_tower_L2Card[n_links_GCTcard][n_towers_per_link][n_GCTcards];
+ int iPhi_tower_L2Card[n_links_GCTcard][n_towers_per_link][n_GCTcards];
+ float energy_cluster_L2Card[n_links_GCTcard][n_clusters_per_link]
+ [n_GCTcards]; // 3 L2 cards send each 48 links with 2 clusters
+ float brem_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
+ int towerID_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
+ int crystalID_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
+ float isolation_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
+ float HE_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
+ int showerShape_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
+ int showerShapeLooseTk_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
+ int photonShowerShape_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
+
+ for (int ii = 0; ii < n_links_GCTcard; ++ii) {
+ for (int jj = 0; jj < n_towers_per_link; ++jj) {
+ for (int ll = 0; ll < n_GCTcards; ++ll) {
+ HCAL_tower_L2Card[ii][jj][ll] = 0;
+ ECAL_tower_L2Card[ii][jj][ll] = 0;
+ iEta_tower_L2Card[ii][jj][ll] = 0;
+ iPhi_tower_L2Card[ii][jj][ll] = 0;
+ }
+ }
+ }
+ for (int ii = 0; ii < n_links_GCTcard; ++ii) {
+ for (int jj = 0; jj < n_clusters_per_link; ++jj) {
+ for (int ll = 0; ll < n_GCTcards; ++ll) {
+ energy_cluster_L2Card[ii][jj][ll] = 0;
+ brem_cluster_L2Card[ii][jj][ll] = 0;
+ towerID_cluster_L2Card[ii][jj][ll] = 0;
+ crystalID_cluster_L2Card[ii][jj][ll] = 0;
+ isolation_cluster_L2Card[ii][jj][ll] = 0;
+ HE_cluster_L2Card[ii][jj][ll] = 0;
+ photonShowerShape_cluster_L2Card[ii][jj][ll] = 0;
+ showerShape_cluster_L2Card[ii][jj][ll] = 0;
+ showerShapeLooseTk_cluster_L2Card[ii][jj][ll] = 0;
+ }
+ }
+ }
+
+ // There is one list of clusters per equivalent of L1 card. We take the 8 highest pt.
+ vector cluster_list_L2[n_towers_halfPhi];
+
+ // Merge clusters on the phi edges
+ for (int ii = 0; ii < n_borders_phi; ++ii) { // 18 borders in phi
+ for (int jj = 0; jj < n_eta_bins; ++jj) { // 2 eta bins
+ int card_left = 2 * ii + jj;
+ int card_right = 2 * ii + jj + 2;
+ if (card_right > 35)
+ card_right = card_right - n_towers_halfPhi;
+ for (int kk = 0; kk < n_clusters_4link; ++kk) { // 12 clusters in the first card. We check the right side
+ if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 50 &&
+ crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 19 &&
+ energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 0) {
+ for (int ll = 0; ll < n_clusters_4link; ++ll) { // We check the 12 clusters in the card on the right
+ if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_per_link &&
+ crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < 5 &&
+ std::abs(
+ 5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) % n_towers_per_link +
+ crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] % 5 -
+ 5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) % n_towers_per_link -
+ crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] % 5) < 2) {
+ if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] >
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) {
+ energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] +=
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right];
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 0;
+ } // The least energetic cluster is merged to the most energetic one
+ else {
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] +=
+ energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left];
+ energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 0;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Bremsstrahlung corrections. Merge clusters on the phi edges depending on pT (pt more than 10 percent, dphi leq 5, deta leq 1)
+ if (do_brem) {
+ for (int ii = 0; ii < n_borders_phi; ++ii) { // 18 borders in phi
+ for (int jj = 0; jj < n_eta_bins; ++jj) { // 2 eta bins
+ int card_left = 2 * ii + jj;
+ int card_right = 2 * ii + jj + 2;
+ if (card_right > 35)
+ card_right = card_right - n_towers_halfPhi;
+ // 12 clusters in the first card. We check the right side crystalID_cluster_L1Card[kk%4][kk/4][card_left]
+ for (int kk = 0; kk < n_clusters_4link; ++kk) {
+ // if the tower is at the edge there might be brem corrections, whatever the crystal ID
+ if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 50 &&
+ energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 0) {
+ for (int ll = 0; ll < n_clusters_4link; ++ll) { // We check the 12 clusters in the card on the right
+ //Distance of 1 max in eta
+ if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_per_link &&
+ std::abs(5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) %
+ n_towers_per_link +
+ crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] % 5 -
+ 5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) %
+ n_towers_per_link -
+ crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] % 5) <= 1) {
+ //Distance of 5 max in phi
+ if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_per_link &&
+ std::abs(5 + crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] / 5 -
+ (crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] / 5)) <= 5) {
+ if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] >
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] &&
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] >
+ 0.10 * energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) {
+ energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] +=
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right];
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 0;
+ brem_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 1;
+ }
+ // The least energetic cluster is merged to the most energetic one, if its pt is at least ten percent
+ else if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_right] >=
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_left] &&
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_left] >
+ 0.10 * energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_right]) {
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] +=
+ energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left];
+ energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 0;
+ brem_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 1;
+ }
+ } //max distance eta
+ } //max distance phi
+ } //max distance phi
+ }
+ }
+ }
+ }
+ }
+
+ // Merge clusters on the eta edges
+ for (int ii = 0; ii < n_borders_eta; ++ii) { // 18 borders in eta
+ int card_bottom = 2 * ii;
+ int card_top = 2 * ii + 1;
+ for (int kk = 0; kk < n_clusters_4link; ++kk) { // 12 clusters in the first card. We check the top side
+ if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] % n_towers_per_link == 16 &&
+ crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] % 5 == n_links_card &&
+ energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] >
+ 0) { // If there is one cluster on the right side of the first card
+ for (int ll = 0; ll < n_clusters_4link; ++ll) { // We check the card on the right
+ if (std::abs(
+ 5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] / n_towers_per_link) +
+ crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] / 5 -
+ 5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] / n_towers_per_link) -
+ crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] / 5) < 2) {
+ if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] >
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_bottom]) {
+ energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] +=
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top];
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] = 0;
+ } else {
+ energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] +=
+ energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom];
+ energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] = 0;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Regroup the new clusters per equivalent of L1 card geometry
+ for (int ii = 0; ii < n_towers_halfPhi; ++ii) {
+ for (int jj = 0; jj < n_clusters_4link; ++jj) {
+ if (energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii] > 0) {
+ mycluster mc1;
+ mc1.cpt = energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
+ mc1.cbrem_ = brem_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
+ mc1.ctowerid_ = towerID_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
+ mc1.ccrystalid_ = crystalID_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
+ mc1.cshowershape_ = showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
+ mc1.cshowershapeloosetk_ = showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
+ mc1.cphotonshowershape_ = photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
+ cluster_list_L2[ii].push_back(mc1);
+ }
+ }
+ std::sort(
+ begin(cluster_list_L2[ii]), end(cluster_list_L2[ii]), [](mycluster a, mycluster b) { return a.cpt > b.cpt; });
+ }
+
+ // If there are more than 8 clusters per equivalent of L1 card we need to put them back in the towers
+ for (int ii = 0; ii < n_towers_halfPhi; ++ii) {
+ for (unsigned int jj = 8; jj < n_clusters_4link && jj < cluster_list_L2[ii].size(); ++jj) {
+ if (cluster_list_L2[ii][jj].cpt > 0) {
+ ECAL_tower_L1Card[cluster_list_L2[ii][jj].ctowerid_ / n_towers_per_link]
+ [cluster_list_L2[ii][jj].ctowerid_ % n_towers_per_link][ii] += cluster_list_L2[ii][jj].cpt;
+ cluster_list_L2[ii][jj].cpt = 0;
+ cluster_list_L2[ii][jj].ctowerid_ = 0;
+ cluster_list_L2[ii][jj].ccrystalid_ = 0;
+ }
+ }
+ }
+
+ // Compute isolation (7*7 ECAL towers) and HCAL energy (5x5 HCAL towers)
+ for (int ii = 0; ii < n_towers_halfPhi; ++ii) { // Loop over the new cluster list (stored in 36x8 format)
+ for (unsigned int jj = 0; jj < 8 && jj < cluster_list_L2[ii].size(); ++jj) {
+ int cluster_etaOfTower_fullDetector = get_towerEta_fromCardTowerInCard(ii, cluster_list_L2[ii][jj].ctowerid_);
+ int cluster_phiOfTower_fullDetector = get_towerPhi_fromCardTowerInCard(ii, cluster_list_L2[ii][jj].ctowerid_);
+ float hcal_nrj = 0.0;
+ float isolation = 0.0;
+ int ntowers = 0;
+ for (int iii = 0; iii < n_towers_halfPhi; ++iii) { // The clusters have to be added to the isolation
+ for (unsigned int jjj = 0; jjj < 8 && jjj < cluster_list_L2[iii].size(); ++jjj) {
+ if (!(iii == ii && jjj == jj)) {
+ int cluster2_eta = get_towerEta_fromCardTowerInCard(iii, cluster_list_L2[iii][jjj].ctowerid_);
+ int cluster2_phi = get_towerPhi_fromCardTowerInCard(iii, cluster_list_L2[iii][jjj].ctowerid_);
+ if (abs(cluster2_eta - cluster_etaOfTower_fullDetector) <= 2 &&
+ (abs(cluster2_phi - cluster_phiOfTower_fullDetector) <= 2 or
+ abs(cluster2_phi - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) {
+ isolation += cluster_list_L2[iii][jjj].cpt;
+ }
+ }
+ }
+ }
+ for (int kk = 0; kk < n_towers_halfPhi; ++kk) { // 36 cards
+ for (int ll = 0; ll < n_links_card; ++ll) { // 4 links per card
+ for (int mm = 0; mm < n_towers_per_link; ++mm) { // 17 towers per link
+ int etaOftower_fullDetector = get_towerEta_fromCardLinkTower(kk, ll, mm);
+ int phiOftower_fullDetector = get_towerPhi_fromCardLinkTower(kk, ll, mm);
+ // First do ECAL
+ // The towers are within 3. Needs to stitch the two phi sides together
+ if (abs(etaOftower_fullDetector - cluster_etaOfTower_fullDetector) <= 2 &&
+ (abs(phiOftower_fullDetector - cluster_phiOfTower_fullDetector) <= 2 or
+ abs(phiOftower_fullDetector - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) {
+ // Remove the column outside of the L2 card: values (0,71), (23,26), (24,21), (47,50), (48,50), (71,2)
+ if (!((cluster_phiOfTower_fullDetector == 0 && phiOftower_fullDetector == 71) or
+ (cluster_phiOfTower_fullDetector == 23 && phiOftower_fullDetector == 26) or
+ (cluster_phiOfTower_fullDetector == 24 && phiOftower_fullDetector == 21) or
+ (cluster_phiOfTower_fullDetector == 47 && phiOftower_fullDetector == 50) or
+ (cluster_phiOfTower_fullDetector == 48 && phiOftower_fullDetector == 45) or
+ (cluster_phiOfTower_fullDetector == 71 && phiOftower_fullDetector == 2))) {
+ isolation += ECAL_tower_L1Card[ll][mm][kk];
+ ntowers++;
+ }
+ }
+ // Now do HCAL
+ // The towers are within 2. Needs to stitch the two phi sides together
+ if (abs(etaOftower_fullDetector - cluster_etaOfTower_fullDetector) <= 2 &&
+ (abs(phiOftower_fullDetector - cluster_phiOfTower_fullDetector) <= 2 or
+ abs(phiOftower_fullDetector - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) {
+ hcal_nrj += HCAL_tower_L1Card[ll][mm][kk];
+ }
+ }
+ }
+ }
+ cluster_list_L2[ii][jj].ciso_ = ((isolation) * (25.0 / ntowers)) / cluster_list_L2[ii][jj].cpt;
+ cluster_list_L2[ii][jj].chovere_ = hcal_nrj / cluster_list_L2[ii][jj].cpt;
+ }
+ }
+
+ //Reformat the information inside the 3 L2 cards
+ //First let's fill the towers
+ for (int ii = 0; ii < n_links_GCTcard; ++ii) {
+ for (int jj = 0; jj < n_towers_per_link; ++jj) {
+ for (int ll = 0; ll < 3; ++ll) {
+ ECAL_tower_L2Card[ii][jj][ll] =
+ ECAL_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
+ HCAL_tower_L2Card[ii][jj][ll] =
+ HCAL_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
+ iEta_tower_L2Card[ii][jj][ll] =
+ iEta_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
+ iPhi_tower_L2Card[ii][jj][ll] =
+ iPhi_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
+ }
+ }
+ }
+
+ //Second let's fill the clusters
+ for (int ii = 0; ii < n_towers_halfPhi; ++ii) { // The cluster list is still in the L1 like geometry
+ for (unsigned int jj = 0; jj < unsigned(cluster_list_L2[ii].size()) && jj < n_clusters_4link; ++jj) {
+ crystalID_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
+ [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ccrystalid_;
+ towerID_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
+ [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ctowerid_;
+ energy_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
+ [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cpt;
+ brem_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
+ [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cbrem_;
+ isolation_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
+ [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ciso_;
+ HE_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
+ [ii / n_clusters_4link] = cluster_list_L2[ii][jj].chovere_;
+ showerShape_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
+ [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cshowershape_;
+ showerShapeLooseTk_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
+ [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cshowershapeloosetk_;
+ photonShowerShape_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
+ [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cphotonshowershape_;
+ }
+ }
+
+ auto L1EGXtalClusters = std::make_unique();
+ auto L1EGammas = std::make_unique();
+ auto L1CaloTowers = std::make_unique();
+
+ // Fill the cluster collection and towers as well
+ for (int ii = 0; ii < n_links_GCTcard; ++ii) { // 48 links
+ for (int ll = 0; ll < n_GCTcards; ++ll) { // 3 cards
+ // For looping over the Towers a few lines below
+ for (int jj = 0; jj < 2; ++jj) { // 2 L1EGs
+ if (energy_cluster_L2Card[ii][jj][ll] > 0.45) {
+ reco::Candidate::PolarLorentzVector p4calibrated(
+ energy_cluster_L2Card[ii][jj][ll],
+ getEta_fromL2LinkCardTowerCrystal(
+ ii, ll, towerID_cluster_L2Card[ii][jj][ll], crystalID_cluster_L2Card[ii][jj][ll]),
+ getPhi_fromL2LinkCardTowerCrystal(
+ ii, ll, towerID_cluster_L2Card[ii][jj][ll], crystalID_cluster_L2Card[ii][jj][ll]),
+ 0.);
+ SimpleCaloHit centerhit;
+ bool is_iso = passes_iso(energy_cluster_L2Card[ii][jj][ll], isolation_cluster_L2Card[ii][jj][ll]);
+ bool is_looseTkiso =
+ passes_looseTkiso(energy_cluster_L2Card[ii][jj][ll], isolation_cluster_L2Card[ii][jj][ll]);
+ bool is_ss = (showerShape_cluster_L2Card[ii][jj][ll] == 1);
+ bool is_photon = (photonShowerShape_cluster_L2Card[ii][jj][ll] == 1) && is_ss && is_iso;
+ bool is_looseTkss = (showerShapeLooseTk_cluster_L2Card[ii][jj][ll] == 1);
+ // All the ID set to Standalone WP! Some dummy values for non calculated variables
+ l1tp2::CaloCrystalCluster cluster(p4calibrated,
+ energy_cluster_L2Card[ii][jj][ll],
+ HE_cluster_L2Card[ii][jj][ll],
+ isolation_cluster_L2Card[ii][jj][ll],
+ centerhit.id(),
+ -1000,
+ float(brem_cluster_L2Card[ii][jj][ll]),
+ -1000,
+ -1000,
+ energy_cluster_L2Card[ii][jj][ll],
+ -1,
+ is_iso&& is_ss,
+ is_iso&& is_ss,
+ is_photon,
+ is_iso&& is_ss,
+ is_looseTkiso&& is_looseTkss,
+ is_iso&& is_ss);
+ // Experimental parameters, don't want to bother with hardcoding them in data format
+ std::map params;
+ params["standaloneWP_showerShape"] = is_ss;
+ params["standaloneWP_isolation"] = is_iso;
+ params["trkMatchWP_showerShape"] = is_looseTkss;
+ params["trkMatchWP_isolation"] = is_looseTkiso;
+ params["TTiEta"] = getToweriEta_fromAbsoluteID(getTower_absoluteEtaID(p4calibrated.eta()));
+ params["TTiPhi"] = getToweriPhi_fromAbsoluteID(getTower_absolutePhiID(p4calibrated.phi()));
+ cluster.setExperimentalParams(params);
+ L1EGXtalClusters->push_back(cluster);
+
+ int standaloneWP = (int)(is_iso && is_ss);
+ int looseL1TkMatchWP = (int)(is_looseTkiso && is_looseTkss);
+ int photonWP = (int)(is_photon);
+ int quality =
+ (standaloneWP * std::pow(2, 0)) + (looseL1TkMatchWP * std::pow(2, 1)) + (photonWP * std::pow(2, 2));
+ L1EGammas->push_back(
+ 0, l1t::EGamma(p4calibrated, p4calibrated.pt(), p4calibrated.eta(), p4calibrated.phi(), quality, 1));
+ }
+ }
+ }
+ }
+
+ for (int ii = 0; ii < n_links_GCTcard; ++ii) { // 48 links
+ for (int ll = 0; ll < n_GCTcards; ++ll) { // 3 cards
+ // Fill the tower collection
+ for (int jj = 0; jj < n_towers_per_link; ++jj) { // 17 TTs
+ l1tp2::CaloTower l1CaloTower;
+ l1CaloTower.setEcalTowerEt(ECAL_tower_L2Card[ii][jj][ll]);
+ l1CaloTower.setHcalTowerEt(HCAL_tower_L2Card[ii][jj][ll]);
+ l1CaloTower.setTowerIEta(getToweriEta_fromAbsoluteID(iEta_tower_L2Card[ii][jj][ll]));
+ l1CaloTower.setTowerIPhi(getToweriPhi_fromAbsoluteID(iPhi_tower_L2Card[ii][jj][ll]));
+ l1CaloTower.setTowerEta(getTowerEta_fromAbsoluteID(iEta_tower_L2Card[ii][jj][ll]));
+ l1CaloTower.setTowerPhi(getTowerPhi_fromAbsoluteID(iPhi_tower_L2Card[ii][jj][ll]));
+ // Some towers have incorrect eta/phi because that wasn't initialized in certain cases, fix these
+ static float constexpr towerEtaUpperUnitialized = -80;
+ static float constexpr towerPhiUpperUnitialized = -90;
+ if (l1CaloTower.towerEta() < towerEtaUpperUnitialized && l1CaloTower.towerPhi() < towerPhiUpperUnitialized) {
+ l1CaloTower.setTowerEta(l1t::CaloTools::towerEta(l1CaloTower.towerIEta()));
+ l1CaloTower.setTowerPhi(l1t::CaloTools::towerPhi(l1CaloTower.towerIEta(), l1CaloTower.towerIPhi()));
+ }
+ l1CaloTower.setIsBarrel(true);
+
+ // Add L1EGs if they match in iEta / iPhi
+ // L1EGs are already pT ordered, we will take the ID info for the leading one, but pT as the sum
+ for (const auto& l1eg : *L1EGXtalClusters) {
+ if (l1eg.experimentalParam("TTiEta") != l1CaloTower.towerIEta())
+ continue;
+ if (l1eg.experimentalParam("TTiPhi") != l1CaloTower.towerIPhi())
+ continue;
+
+ int n_L1eg = l1CaloTower.nL1eg();
+ l1CaloTower.setNL1eg(n_L1eg++);
+ l1CaloTower.setL1egTowerEt(l1CaloTower.l1egTowerEt() + l1eg.pt());
+ // Don't record L1EG quality info for subleading L1EG
+ if (l1CaloTower.nL1eg() > 1)
+ continue;
+ l1CaloTower.setL1egTrkSS(l1eg.experimentalParam("trkMatchWP_showerShape"));
+ l1CaloTower.setL1egTrkIso(l1eg.experimentalParam("trkMatchWP_isolation"));
+ l1CaloTower.setL1egStandaloneSS(l1eg.experimentalParam("standaloneWP_showerShape"));
+ l1CaloTower.setL1egStandaloneIso(l1eg.experimentalParam("standaloneWP_isolation"));
+ }
+
+ L1CaloTowers->push_back(l1CaloTower);
+ }
+ }
+ }
+
+ iEvent.put(std::move(L1EGXtalClusters));
+ iEvent.put(std::move(L1EGammas));
+ iEvent.put(std::move(L1CaloTowers));
+}
+
+bool L1EGCrystalClusterEmulatorProducer::passes_iso(float pt, float iso) {
+ if (pt < slideIsoPtThreshold) {
+ if (!((a0_80 - a1_80 * pt) > iso))
+ return false;
+ } else {
+ if (iso > a0)
+ return false;
+ }
+ return true;
+}
+
+bool L1EGCrystalClusterEmulatorProducer::passes_looseTkiso(float pt, float iso) {
+ return (b0 + b1 * std::exp(-b2 * pt) > iso);
+}
+
+bool L1EGCrystalClusterEmulatorProducer::passes_ss(float pt, float ss) {
+ return ((c0_ss + c1_ss * std::exp(-c2_ss * pt)) <= ss);
+}
+
+bool L1EGCrystalClusterEmulatorProducer::passes_photon(float pt, float pss) { return (pss > d0 - d1 * pt); }
+
+bool L1EGCrystalClusterEmulatorProducer::passes_looseTkss(float pt, float ss) {
+ return ((e0_looseTkss - e1_looseTkss * std::exp(-e2_looseTkss * pt)) <= ss);
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(L1EGCrystalClusterEmulatorProducer);
diff --git a/L1Trigger/L1CaloTrigger/plugins/L1EGammaEEProducer.cc b/L1Trigger/L1CaloTrigger/plugins/L1EGammaEEProducer.cc
new file mode 100644
index 0000000000000..c2fbca3361aeb
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/plugins/L1EGammaEEProducer.cc
@@ -0,0 +1,169 @@
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "DataFormats/L1Trigger/interface/EGamma.h"
+#include "L1Trigger/L1CaloTrigger/interface/L1EGammaEECalibrator.h"
+#include "DataFormats/Math/interface/deltaPhi.h"
+
+namespace l1tp2 {
+ // we sort the clusters in pt
+ bool compare_cluster_pt(const l1t::HGCalMulticluster *cl1, const l1t::HGCalMulticluster *cl2) {
+ return cl1->pt() > cl2->pt();
+ }
+}; // namespace l1tp2
+
+int etaBin(const l1t::HGCalMulticluster *cl) {
+ static float constexpr eta_min = 1.;
+ static float constexpr eta_max = 4.;
+ static unsigned constexpr n_eta_bins = 150;
+ int eta_bin = floor((std::abs(cl->eta()) - eta_min) / ((eta_max - eta_min) / n_eta_bins));
+ if (cl->eta() < 0)
+ return -1 * eta_bin; // bin 0 doesn't exist
+ return eta_bin;
+}
+
+int get_phi_bin(const l1t::HGCalMulticluster *cl) {
+ static constexpr float phi_min = -M_PI;
+ static constexpr float phi_max = M_PI;
+ static constexpr unsigned n_phi_bins = 63;
+ return floor(std::abs(reco::deltaPhi(cl->phi(), phi_min)) / ((phi_max - phi_min) / n_phi_bins));
+}
+
+pair get_eta_phi_bin(const l1t::HGCalMulticluster *cl) { return std::make_pair(etaBin(cl), get_phi_bin(cl)); }
+
+class L1EGammaEEProducer : public edm::stream::EDProducer<> {
+public:
+ explicit L1EGammaEEProducer(const edm::ParameterSet &);
+
+private:
+ void produce(edm::Event &, const edm::EventSetup &) override;
+
+ edm::EDGetToken multiclusters_token_;
+ L1EGammaEECalibrator calibrator_;
+};
+
+L1EGammaEEProducer::L1EGammaEEProducer(const edm::ParameterSet &iConfig)
+ : multiclusters_token_(
+ consumes(iConfig.getParameter("Multiclusters"))),
+ calibrator_(iConfig.getParameter("calibrationConfig")) {
+ produces>("L1EGammaCollectionBXVWithCuts");
+}
+
+void L1EGammaEEProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
+ float minEt_ = 0;
+
+ std::unique_ptr> l1EgammaBxCollection(new l1t::EGammaBxCollection);
+
+ // retrieve clusters 3D
+ edm::Handle multiclusters_h;
+ iEvent.getByToken(multiclusters_token_, multiclusters_h);
+ const l1t::HGCalMulticlusterBxCollection &multiclusters = *multiclusters_h;
+
+ std::vector selected_multiclusters;
+ std::map, std::vector> etaphi_bins;
+
+ // here we loop on the TPGs
+ for (auto cl3d = multiclusters.begin(0); cl3d != multiclusters.end(0); cl3d++) {
+ if (cl3d->hwQual()) {
+ if (cl3d->et() > minEt_) {
+ int hw_quality = 1; // baseline EG ID passed
+ if (std::abs(cl3d->eta()) >= 1.52) {
+ hw_quality = 2; // baseline EG ID passed + cleanup of transition region
+ }
+
+ float calib_factor = calibrator_.calibrationFactor(cl3d->pt(), cl3d->eta());
+ l1t::EGamma eg =
+ l1t::EGamma(reco::Candidate::PolarLorentzVector(cl3d->pt() / calib_factor, cl3d->eta(), cl3d->phi(), 0.));
+ eg.setHwQual(hw_quality);
+ eg.setHwIso(1);
+ eg.setIsoEt(-1); // just temporarily as a dummy value
+ l1EgammaBxCollection->push_back(0, eg);
+ if (hw_quality == 2) {
+ // we build the EM interpreted EG object
+ l1t::EGamma eg_emint = l1t::EGamma(reco::Candidate::PolarLorentzVector(
+ cl3d->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM), cl3d->eta(), cl3d->phi(), 0.));
+ eg_emint.setHwQual(4);
+ eg_emint.setHwIso(1);
+ eg_emint.setIsoEt(-1); // just temporarily as a dummy value
+ l1EgammaBxCollection->push_back(0, eg_emint);
+ // we also prepare for the brem recovery procedure
+ selected_multiclusters.push_back(&(*cl3d));
+ auto eta_phi_bin = get_eta_phi_bin(&(*cl3d));
+ auto bucket = etaphi_bins.find(eta_phi_bin);
+ if (bucket == etaphi_bins.end()) {
+ std::vector vec;
+ vec.push_back(&(*cl3d));
+ etaphi_bins[eta_phi_bin] = vec;
+ } else {
+ bucket->second.push_back(&(*cl3d));
+ }
+ }
+ }
+ }
+ }
+
+ std::sort(selected_multiclusters.begin(), selected_multiclusters.end(), l1tp2::compare_cluster_pt);
+ std::set used_clusters;
+ for (const auto &cl3d : selected_multiclusters) {
+ if (used_clusters.find(cl3d) == used_clusters.end()) {
+ float pt = cl3d->pt();
+ // we drop the Had component of the energy
+ if (cl3d->hOverE() != -1)
+ pt = cl3d->pt() / (1 + cl3d->hOverE());
+ reco::Candidate::PolarLorentzVector mom(pt, cl3d->eta(), cl3d->phi(), 0.);
+ reco::Candidate::PolarLorentzVector mom_eint(
+ cl3d->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM), cl3d->eta(), cl3d->phi(), 0.);
+
+ // this is not yet used
+ used_clusters.insert(cl3d);
+ auto eta_phi_bin = get_eta_phi_bin(cl3d);
+
+ for (int eta_bin : {eta_phi_bin.first - 1, eta_phi_bin.first, eta_phi_bin.first + 1}) {
+ for (int phi_bin : {eta_phi_bin.second - 1, eta_phi_bin.second, eta_phi_bin.second + 1}) {
+ auto bucket = etaphi_bins.find(std::make_pair(eta_bin, phi_bin));
+ if (bucket != etaphi_bins.end()) {
+ // this bucket is not empty
+ for (const auto &other_cl_ptr : bucket->second) {
+ if (used_clusters.find(other_cl_ptr) == used_clusters.end()) {
+ if (std::abs(other_cl_ptr->eta() - cl3d->eta()) < 0.02) {
+ if (std::abs(reco::deltaPhi(other_cl_ptr->phi(), cl3d->phi())) < 0.1) {
+ float pt_other = other_cl_ptr->pt();
+ if (other_cl_ptr->hOverE() != -1)
+ pt_other = other_cl_ptr->pt() / (1 + other_cl_ptr->hOverE());
+ mom += reco::Candidate::PolarLorentzVector(pt_other, other_cl_ptr->eta(), other_cl_ptr->phi(), 0.);
+ mom_eint += reco::Candidate::PolarLorentzVector(
+ other_cl_ptr->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM),
+ other_cl_ptr->eta(),
+ other_cl_ptr->phi(),
+ 0.);
+ used_clusters.insert(other_cl_ptr);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ float calib_factor = calibrator_.calibrationFactor(mom.pt(), mom.eta());
+ l1t::EGamma eg =
+ l1t::EGamma(reco::Candidate::PolarLorentzVector(mom.pt() / calib_factor, mom.eta(), mom.phi(), 0.));
+ eg.setHwQual(3);
+ eg.setHwIso(1);
+ l1EgammaBxCollection->push_back(0, eg);
+
+ l1t::EGamma eg_emint_brec =
+ l1t::EGamma(reco::Candidate::PolarLorentzVector(mom_eint.pt(), mom_eint.eta(), mom_eint.phi(), 0.));
+ eg_emint_brec.setHwQual(5);
+ eg_emint_brec.setHwIso(1);
+ l1EgammaBxCollection->push_back(0, eg_emint_brec);
+ }
+ }
+
+ iEvent.put(std::move(l1EgammaBxCollection));
+}
+
+DEFINE_FWK_MODULE(L1EGammaEEProducer);
diff --git a/L1Trigger/L1CaloTrigger/python/L1EGammaCrystalsEmulatorProducer_cfi.py b/L1Trigger/L1CaloTrigger/python/L1EGammaCrystalsEmulatorProducer_cfi.py
new file mode 100644
index 0000000000000..750d12e168260
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/python/L1EGammaCrystalsEmulatorProducer_cfi.py
@@ -0,0 +1,30 @@
+import FWCore.ParameterSet.Config as cms
+
+L1EGammaClusterEmuProducer = cms.EDProducer("L1EGCrystalClusterEmulatorProducer",
+ ecalTPEB = cms.InputTag("simEcalEBTriggerPrimitiveDigis","","HLT"),
+ hcalTP = cms.InputTag("simHcalTriggerPrimitiveDigis","","HLT"),
+ calib = cms.PSet(
+
+ etaBins = cms.vdouble( 0.087 , 0.174 , 0.261 , 0.348 , 0.435 , 0.522 , 0.609 , 0.696 , 0.783 , 0.870 , 0.957 , 1.044 , 1.131 , 1.218 , 1.305 , 1.392 , 1.479),
+ ptBins = cms.vdouble( 12, 20, 30, 40, 55, 90, 1e6),
+ scale = cms.vdouble(
+ # pt < 12
+ 1.18*1.1 ,1.17*1.1 ,1.19*1.1 ,1.18*1.1 ,1.19*1.1 ,1.19*1.1 ,1.19*1.1 ,1.18*1.1 ,1.19*1.1 ,1.18*1.1 ,1.19*1.1 ,1.19*1.1 ,1.19*1.1 ,1.20*1.1 ,1.19*1.1 ,1.20*1.1 ,1.19*1.1,
+ # pt < 20
+ 1.14*1.03 ,1.13*1.03 ,1.13*1.03 ,1.13*1.03 ,1.13*1.03 ,1.13*1.03 ,1.13*1.03 ,1.14*1.03 ,1.14*1.03 ,1.13*1.03 ,1.13*1.03 ,1.14*1.03 ,1.13*1.03 ,1.13*1.03 ,1.14*1.03 ,1.14*1.03 ,1.12*1.03,
+ # pt < 30
+ 1.11 ,1.11 ,1.11 ,1.11 ,1.11 ,1.11 ,1.11 ,1.11 ,1.11 ,1.11 ,1.11 ,1.11 ,1.11 ,1.11 ,1.11 ,1.11 ,1.10,
+ # pt < 40
+ 1.09 ,1.09 ,1.09 ,1.09 ,1.09 ,1.09 ,1.09 ,1.09 ,1.09 ,1.09 ,1.09 ,1.09 ,1.09 ,1.09 ,1.09 ,1.09 ,1.09,
+ # pt < 55
+ 1.07 ,1.07 ,1.07 ,1.07 ,1.07 ,1.07 ,1.07 ,1.08 ,1.07 ,1.07 ,1.08 ,1.08 ,1.07 ,1.08 ,1.08 ,1.08 ,1.08,
+ # pt < 90
+ 1.06 ,1.06 ,1.06 ,1.06 ,1.05 ,1.05 ,1.06 ,1.06 ,1.06 ,1.06 ,1.06 ,1.06 ,1.06 ,1.06 ,1.06 ,1.06 ,1.06,
+ # pt < 1e6
+ 1.04 ,1.04 ,1.04 ,1.04 ,1.05 ,1.04 ,1.05 ,1.05 ,1.05 ,1.05 ,1.05 ,1.05 ,1.05 ,1.05 ,1.05 ,1.05 ,1.05,
+
+ ),
+
+ ),
+)
+
diff --git a/L1Trigger/L1CaloTrigger/python/L1TowerCalibrationProducer_cfi.py b/L1Trigger/L1CaloTrigger/python/L1TowerCalibrationProducer_cfi.py
new file mode 100644
index 0000000000000..dcadde6ad3a59
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/python/L1TowerCalibrationProducer_cfi.py
@@ -0,0 +1,184 @@
+import FWCore.ParameterSet.Config as cms
+
+L1TowerCalibrationProducer = cms.EDProducer("L1TowerCalibrator",
+ # Choosen settings 6 March 2019, 10_3_X MTD samples
+ HcalTpEtMin = cms.double(0.5),
+ EcalTpEtMin = cms.double(0.5),
+ HGCalHadTpEtMin = cms.double(0.25),
+ HGCalEmTpEtMin = cms.double(0.25),
+ HFTpEtMin = cms.double(0.5),
+ puThreshold = cms.double(5.0),
+ puThresholdL1eg = cms.double(2.0),
+ puThresholdHcalMin = cms.double(1.0),
+ puThresholdHcalMax = cms.double(2.0),
+ puThresholdEcalMin = cms.double(0.75),
+ puThresholdEcalMax = cms.double(1.5),
+ puThresholdHGCalEMMin = cms.double(1.25),
+ puThresholdHGCalEMMax = cms.double(1.75),
+ puThresholdHGCalHadMin = cms.double(0.75),
+ puThresholdHGCalHadMax = cms.double(1.25),
+ puThresholdHFMin = cms.double(10.0),
+ puThresholdHFMax = cms.double(15.0),
+
+ puThresholdEcal = cms.double(2.0),
+ puThresholdHcal = cms.double(3.0),
+
+ barrelSF = cms.double(1.0),
+ hgcalSF = cms.double(1.0),
+ hfSF = cms.double(1.0),
+ debug = cms.bool(False),
+ skipCalibrations = cms.bool(False),
+ l1CaloTowers = cms.InputTag("L1EGammaClusterEmuProducer",),
+ L1HgcalTowersInputTag = cms.InputTag("hgcalTowerProducer:HGCalTowerProcessor"),
+ hcalDigis = cms.InputTag("simHcalTriggerPrimitiveDigis"),
+ nHits_to_nvtx_params = cms.VPSet( # Parameters derived on 6 March 2019 on 10_3_X MTD samples
+ cms.PSet(
+ fit = cms.string( "hf" ),
+ params = cms.vdouble( 165.706, 0.153 )
+ ),
+ cms.PSet(
+ fit = cms.string( "ecal" ),
+ params = cms.vdouble( 168.055, 0.377 )
+ ),
+ cms.PSet(
+ fit = cms.string( "hgcalEM" ),
+ params = cms.vdouble( 157.522, 0.090 )
+ ),
+ cms.PSet(
+ fit = cms.string( "hgcalHad" ),
+ params = cms.vdouble( 159.295, 0.178 )
+ ),
+ cms.PSet(
+ fit = cms.string( "hcal" ),
+ params = cms.vdouble( 168.548, 0.362 )
+ ),
+ ),
+
+ nvtx_to_PU_sub_params = cms.VPSet(
+ cms.PSet(
+ calo = cms.string( "ecal" ),
+ iEta = cms.string( "er1to3" ),
+ params = cms.vdouble( 0.008955, 0.000298 )
+ ),
+ cms.PSet(
+ calo = cms.string( "ecal" ),
+ iEta = cms.string( "er4to6" ),
+ params = cms.vdouble( 0.009463, 0.000256 )
+ ),
+ cms.PSet(
+ calo = cms.string( "ecal" ),
+ iEta = cms.string( "er7to9" ),
+ params = cms.vdouble( 0.008783, 0.000293 )
+ ),
+ cms.PSet(
+ calo = cms.string( "ecal" ),
+ iEta = cms.string( "er10to12" ),
+ params = cms.vdouble( 0.009308, 0.000255 )
+ ),
+ cms.PSet(
+ calo = cms.string( "ecal" ),
+ iEta = cms.string( "er13to15" ),
+ params = cms.vdouble( 0.009290, 0.000221 )
+ ),
+ cms.PSet(
+ calo = cms.string( "ecal" ),
+ iEta = cms.string( "er16to18" ),
+ params = cms.vdouble( 0.009282, 0.000135 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hcal" ),
+ iEta = cms.string( "er1to3" ),
+ params = cms.vdouble( 0.009976, 0.000377 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hcal" ),
+ iEta = cms.string( "er4to6" ),
+ params = cms.vdouble( 0.009803, 0.000394 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hcal" ),
+ iEta = cms.string( "er7to9" ),
+ params = cms.vdouble( 0.009654, 0.000429 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hcal" ),
+ iEta = cms.string( "er10to12" ),
+ params = cms.vdouble( 0.009107, 0.000528 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hcal" ),
+ iEta = cms.string( "er13to15" ),
+ params = cms.vdouble( 0.008367, 0.000652 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hcal" ),
+ iEta = cms.string( "er16to18" ),
+ params = cms.vdouble( 0.009681, 0.000096 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hgcalEM" ),
+ iEta = cms.string( "er1p4to1p8" ),
+ params = cms.vdouble( -0.011772, 0.004142 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hgcalEM" ),
+ iEta = cms.string( "er1p8to2p1" ),
+ params = cms.vdouble( -0.015488, 0.005410 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hgcalEM" ),
+ iEta = cms.string( "er2p1to2p4" ),
+ params = cms.vdouble( -0.021150, 0.006078 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hgcalEM" ),
+ iEta = cms.string( "er2p4to2p7" ),
+ params = cms.vdouble( -0.015705, 0.005339 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hgcalEM" ),
+ iEta = cms.string( "er2p7to3p1" ),
+ params = cms.vdouble( -0.018492, 0.005620 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hgcalHad" ),
+ iEta = cms.string( "er1p4to1p8" ),
+ params = cms.vdouble( 0.005675, 0.000615 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hgcalHad" ),
+ iEta = cms.string( "er1p8to2p1" ),
+ params = cms.vdouble( 0.004560, 0.001099 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hgcalHad" ),
+ iEta = cms.string( "er2p1to2p4" ),
+ params = cms.vdouble( 0.000036, 0.001608 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hgcalHad" ),
+ iEta = cms.string( "er2p4to2p7" ),
+ params = cms.vdouble( 0.000869, 0.001754 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hgcalHad" ),
+ iEta = cms.string( "er2p7to3p1" ),
+ params = cms.vdouble( -0.006574, 0.003134 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hf" ),
+ iEta = cms.string( "er29to33" ),
+ params = cms.vdouble( -0.203291, 0.044096 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hf" ),
+ iEta = cms.string( "er34to37" ),
+ params = cms.vdouble( -0.210922, 0.045628 )
+ ),
+ cms.PSet(
+ calo = cms.string( "hf" ),
+ iEta = cms.string( "er38to41" ),
+ params = cms.vdouble( -0.229562, 0.050560 )
+ ),
+ )
+)
diff --git a/L1Trigger/L1CaloTrigger/python/l1EGammaCrystalsProducer_cfi.py b/L1Trigger/L1CaloTrigger/python/l1EGammaCrystalsProducer_cfi.py
new file mode 100644
index 0000000000000..3e90790d23df4
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/python/l1EGammaCrystalsProducer_cfi.py
@@ -0,0 +1,17 @@
+import FWCore.ParameterSet.Config as cms
+
+l1EGammaCrystalsProducer = cms.EDProducer("L1EGCrystalClusterProducer",
+ EtminForStore = cms.double(0.),
+ EcalTpEtMin = cms.untracked.double(0.5), # 500 MeV default per each Ecal TP
+ EtMinForSeedHit = cms.untracked.double(1.0), # 1 GeV decault for seed hit
+ debug = cms.untracked.bool(False),
+ useRecHits = cms.untracked.bool(False),
+ doBremClustering = cms.untracked.bool(True), # Should always be True when using for E/Gamma
+ ecalTPEB = cms.InputTag("simEcalEBTriggerPrimitiveDigis","","HLT"),
+ ecalRecHitEB = cms.InputTag("ecalRecHit","EcalRecHitsEB","RECO"),
+ hcalRecHit = cms.InputTag("hbhereco"),
+ hcalTP = cms.InputTag("simHcalTriggerPrimitiveDigis","","HLT"),
+ useTowerMap = cms.untracked.bool(False)
+)
+
+
diff --git a/L1Trigger/L1CaloTrigger/python/l1EGammaEEProducer_cfi.py b/L1Trigger/L1CaloTrigger/python/l1EGammaEEProducer_cfi.py
new file mode 100644
index 0000000000000..161e8e8faa27e
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/python/l1EGammaEEProducer_cfi.py
@@ -0,0 +1,5 @@
+import FWCore.ParameterSet.Config as cms
+
+l1EGammaEEProducer = cms.EDProducer("L1EGammaEEProducer",
+ calibrationConfig = cms.PSet(calibrationFile = cms.FileInPath('L1Trigger/L1TCalorimeter/data/calib_ee_v1.json')),
+ Multiclusters=cms.InputTag('hgcalBackEndLayer2Producer:HGCalBackendLayer2Processor3DClustering'))
diff --git a/L1Trigger/L1CaloTrigger/python/ntuple_cfi.py b/L1Trigger/L1CaloTrigger/python/ntuple_cfi.py
new file mode 100644
index 0000000000000..0abf49530f41d
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/python/ntuple_cfi.py
@@ -0,0 +1,16 @@
+import FWCore.ParameterSet.Config as cms
+
+ntuple_egammaEE = cms.PSet(
+ NtupleName = cms.string('L1TriggerNtupleEgammaEE'),
+ EgammaEE = cms.InputTag('l1EGammaEEProducer:L1EGammaCollectionBXVWithCuts')
+)
+
+ntuple_TTTracks = cms.PSet(
+ NtupleName = cms.string('L1TriggerNtupleTrackTrigger'),
+ TTTracks = cms.InputTag("TTTracksFromTrackletEmulation", "Level1TTTracks")
+)
+
+ntuple_tkEle = cms.PSet(
+ NtupleName = cms.string('L1TriggerNtupleTkElectrons'),
+ TkElectrons = cms.InputTag("L1TkElectrons","EG")
+)
diff --git a/L1Trigger/L1CaloTrigger/python/test_L1EGCrystalClusterEmulator_cfg.py b/L1Trigger/L1CaloTrigger/python/test_L1EGCrystalClusterEmulator_cfg.py
new file mode 100644
index 0000000000000..12d320d1ee80b
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/python/test_L1EGCrystalClusterEmulator_cfg.py
@@ -0,0 +1,76 @@
+import FWCore.ParameterSet.Config as cms
+
+from Configuration.StandardSequences.Eras import eras
+
+process = cms.Process("L1AlgoTest",eras.Phase2C4)
+
+process.load('Configuration.StandardSequences.Services_cff')
+process.load("FWCore.MessageService.MessageLogger_cfi")
+process.load('Configuration.EventContent.EventContent_cff')
+process.MessageLogger.categories = cms.untracked.vstring('L1EGRateStudies', 'FwkReport')
+process.MessageLogger.cerr.FwkReport = cms.untracked.PSet(
+ reportEvery = cms.untracked.int32(1)
+)
+
+process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(10) )
+
+process.source = cms.Source("PoolSource",
+ # Set to do test run on official Phase-2 L1T Ntuples
+ fileNames = cms.untracked.vstring('file:root://cmsxrootd.fnal.gov///store/mc/PhaseIISpring17D/SingleE_FlatPt-8to100/GEN-SIM-DIGI-RAW/PU200_90X_upgrade2023_realistic_v9-v1/120000/04B4BF1D-1E2C-E711-BE1C-7845C4FC39D1.root'),
+ inputCommands = cms.untracked.vstring(
+ "keep *",
+ "drop l1tEMTFHitExtras_simEmtfDigis_CSC_HLT",
+ "drop l1tEMTFHitExtras_simEmtfDigis_RPC_HLT",
+ "drop l1tEMTFTrackExtras_simEmtfDigis__HLT",
+ )
+)
+
+# All this stuff just runs the various EG algorithms that we are studying
+
+# ---- Global Tag :
+process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff')
+from Configuration.AlCa.GlobalTag import GlobalTag
+process.GlobalTag = GlobalTag(process.GlobalTag, '100X_upgrade2023_realistic_v1', '')
+
+# Choose a 2023 geometry!
+process.load('Configuration.Geometry.GeometryExtended2023D17Reco_cff')
+process.load('Configuration.StandardSequences.MagneticField_cff')
+
+# Add HCAL Transcoder
+process.load('SimCalorimetry.HcalTrigPrimProducers.hcaltpdigi_cff')
+process.load('CalibCalorimetry.CaloTPG.CaloTPGTranscoder_cfi')
+
+
+
+
+# --------------------------------------------------------------------------------------------
+#
+# ---- Produce the L1EGCrystal clusters using Emulator
+
+
+process.load('L1Trigger.L1CaloTrigger.L1EGammaCrystalsEmulatorProducer_cfi')
+
+process.pL1EG = cms.Path( process.L1EGammaClusterEmuProducer )
+
+
+
+
+process.Out = cms.OutputModule( "PoolOutputModule",
+ fileName = cms.untracked.string( "l1egCrystalTest.root" ),
+ fastCloning = cms.untracked.bool( False ),
+ outputCommands = cms.untracked.vstring(
+ "keep *_L1EGammaClusterEmuProducer_*_*",
+ "keep *_TriggerResults_*_*",
+ "keep *_simHcalTriggerPrimitiveDigis_*_*",
+ "keep *_EcalEBTrigPrimProducer_*_*"
+ )
+)
+
+process.end = cms.EndPath( process.Out )
+
+
+
+dump_file = open("dump_file.py", "w")
+dump_file.write(process.dumpPython())
+
+
diff --git a/L1Trigger/L1CaloTrigger/src/L1EGammaEECalibrator.cc b/L1Trigger/L1CaloTrigger/src/L1EGammaEECalibrator.cc
new file mode 100644
index 0000000000000..ddd9ee590a12c
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/src/L1EGammaEECalibrator.cc
@@ -0,0 +1,61 @@
+#include "L1Trigger/L1CaloTrigger/interface/L1EGammaEECalibrator.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/FileInPath.h"
+#include "boost/property_tree/ptree.hpp"
+#include "boost/property_tree/json_parser.hpp"
+#include
+
+namespace l1tp2 {
+ std::vector as_vector(boost::property_tree::ptree const& pt,
+ boost::property_tree::ptree::key_type const& key) {
+ std::vector ret;
+ for (const auto& item : pt.get_child(key))
+ ret.push_back(item.second.get_value());
+ return ret;
+ }
+}; // namespace l1tp2
+
+L1EGammaEECalibrator::L1EGammaEECalibrator(const edm::ParameterSet& pset) {
+ //read the JSON file and populate the eta - pt bins and the value container
+ boost::property_tree::ptree calibration_map;
+ read_json(pset.getParameter("calibrationFile").fullPath(), calibration_map);
+
+ auto eta_l = l1tp2::as_vector(calibration_map, "eta_l");
+ std::copy(eta_l.begin(), eta_l.end(), std::inserter(eta_bins, eta_bins.end()));
+ auto eta_h = l1tp2::as_vector(calibration_map, "eta_h");
+ eta_bins.insert(eta_h.back());
+
+ auto pt_l = l1tp2::as_vector(calibration_map, "pt_l");
+ std::copy(pt_l.begin(), pt_l.end(), std::inserter(pt_bins, pt_bins.end()));
+ auto pt_h = l1tp2::as_vector(calibration_map, "pt_h");
+ pt_bins.insert(pt_h.back());
+
+ auto calib_data = l1tp2::as_vector(calibration_map, "calib");
+ auto n_bins_eta = eta_bins.size();
+ auto n_bins_pt = pt_bins.size();
+ calib_factors.reserve(n_bins_eta * n_bins_pt);
+ for (auto calib_f = calib_data.begin(); calib_f != calib_data.end(); ++calib_f) {
+ auto index = calib_f - calib_data.begin();
+ int eta_bin = etaBin(eta_l[index]);
+ int pt_bin = ptBin(pt_l[index]);
+ calib_factors[(eta_bin * n_bins_pt) + pt_bin] = *calib_f;
+ }
+}
+
+int L1EGammaEECalibrator::bin(const std::set& container, float value) const {
+ auto bin_l = container.upper_bound(value);
+ if (bin_l == container.end()) {
+ // value not mapped to any bin
+ return -1;
+ }
+ return std::distance(container.begin(), bin_l) - 1;
+}
+
+float L1EGammaEECalibrator::calibrationFactor(const float& pt, const float& eta) const {
+ int bin_eta = etaBin(eta);
+ int bin_pt = ptBin(pt);
+ if (bin_eta == -1 || bin_pt == -1)
+ return 1.;
+ auto n_bins_pt = pt_bins.size();
+ return calib_factors[(bin_eta * n_bins_pt) + bin_pt];
+}
diff --git a/L1Trigger/L1CaloTrigger/src/ParametricCalibration.cc b/L1Trigger/L1CaloTrigger/src/ParametricCalibration.cc
new file mode 100644
index 0000000000000..60560f431fc01
--- /dev/null
+++ b/L1Trigger/L1CaloTrigger/src/ParametricCalibration.cc
@@ -0,0 +1,49 @@
+#include "L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h"
+
+namespace l1tp2 {
+ ParametricCalibration::ParametricCalibration(const edm::ParameterSet& cpset) {
+ std::vector etaBins = cpset.getParameter>("etaBins");
+ std::vector ptBins = cpset.getParameter>("ptBins");
+ std::vector scale = cpset.getParameter>("scale");
+ etas.insert(etas.end(), etaBins.begin(), etaBins.end());
+ pts.insert(pts.end(), ptBins.begin(), ptBins.end());
+ scales.insert(scales.end(), scale.begin(), scale.end());
+
+ if (pts.size() * etas.size() != scales.size())
+ throw cms::Exception("Configuration",
+ "Bad number of calibration scales, pts.size() * etas.size() != scales.size()");
+ }
+
+ // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
+ void ParametricCalibration::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+ edm::ParameterSetDescription desc;
+ desc.setComment("");
+ desc.addUntracked>("etaBins", std::vector{});
+ desc.addUntracked>("ptBins", std::vector{});
+ desc.addUntracked>("scale", std::vector{});
+ descriptions.add("createIdealTkAlRecords", desc);
+ }
+
+ float ParametricCalibration::operator()(const float pt, const float abseta) const {
+ int ptBin = -1;
+ for (unsigned int i = 0, n = pts.size(); i < n; ++i) {
+ if (pt < pts[i]) {
+ ptBin = i;
+ break;
+ }
+ }
+ int etaBin = -1;
+ for (unsigned int i = 0, n = etas.size(); i < n; ++i) {
+ if (abseta < etas[i]) {
+ etaBin = i;
+ break;
+ }
+ }
+
+ if (ptBin == -1 || etaBin == -1)
+ return 1;
+ else
+ return scales[ptBin * etas.size() + etaBin];
+ }
+
+}; // namespace l1tp2