Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:36:49 

0001 // -*- C++ -*-
0002 //
0003 // Package: L1CaloTrigger
0004 // Class: L1EGammaCrystalsEmulatorProducer
0005 //
0006 /**\class L1EGammaCrystalsEmulatorProducer L1EGammaCrystalsEmulatorProducer.cc L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc
0007 Description: Produces crystal clusters using crystal-level information and hardware constraints
0008 Implementation:
0009 [Notes on implementation]
0010 */
0011 //
0012 // Original Author: Cecile Caillol
0013 // Created: Tue Aug 10 2018
0014 //
0015 // $Id$
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 #include <array>
0022 #include <iostream>
0023 #include <cmath>
0024 
0025 // user include files
0026 #include "FWCore/Framework/interface/stream/EDProducer.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/EventSetup.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 
0031 #include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h"
0032 #include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h"
0033 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0034 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0035 #include "Geometry/HcalTowerAlgo/interface/HcalTrigTowerGeometry.h"
0036 #include "FWCore/Framework/interface/MakerMacros.h"
0037 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0038 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0039 
0040 // ECAL TPs
0041 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0042 
0043 // HCAL TPs
0044 #include "DataFormats/HcalDigi/interface/HcalTriggerPrimitiveDigi.h"
0045 
0046 // Output tower collection
0047 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloCrystalCluster.h"
0048 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h"
0049 #include "DataFormats/L1Trigger/interface/EGamma.h"
0050 
0051 #include "L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h"
0052 #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
0053 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0054 
0055 static constexpr bool do_brem = true;
0056 
0057 static constexpr int n_eta_bins = 2;
0058 static constexpr int n_borders_phi = 18;
0059 static constexpr int n_borders_eta = 18;
0060 static constexpr int n_clusters_max = 5;
0061 static constexpr int n_clusters_link = 3;
0062 static constexpr int n_clusters_4link = 4 * 3;
0063 static constexpr int n_crystals_towerEta = 5;
0064 static constexpr int n_crystals_towerPhi = 5;
0065 static constexpr int n_crystals_3towers = 3 * 5;
0066 static constexpr int n_crystals_2towers = 2 * 5;
0067 static constexpr int n_towers_per_link = 17;
0068 static constexpr int n_clusters_per_link = 2;
0069 static constexpr int n_clusters_per_L1card = 8;
0070 static constexpr int n_towers_Eta = 34;
0071 static constexpr int n_towers_Phi = 72;
0072 static constexpr int n_towers_halfPhi = 36;
0073 static constexpr int n_links_card = 4;
0074 static constexpr int n_links_GCTcard = 48;
0075 static constexpr int n_GCTcards = 3;
0076 static constexpr float ECAL_eta_range = 1.4841;
0077 static constexpr float half_crystal_size = 0.00873;
0078 static constexpr float slideIsoPtThreshold = 80;
0079 static constexpr float plateau_ss = 130.0;
0080 static constexpr float a0_80 = 0.85, a1_80 = 0.0080, a0 = 0.21;                        // passes_iso
0081 static constexpr float b0 = 0.38, b1 = 1.9, b2 = 0.05;                                 //passes_looseTkiso
0082 static constexpr float c0_ss = 0.94, c1_ss = 0.052, c2_ss = 0.044;                     //passes_ss
0083 static constexpr float d0 = 0.96, d1 = 0.0003;                                         //passes_photon
0084 static constexpr float e0_looseTkss = 0.944, e1_looseTkss = 0.65, e2_looseTkss = 0.4;  //passes_looseTkss
0085 static constexpr float cut_500_MeV = 0.5;
0086 
0087 // absolue IDs range from 0-33
0088 // 0-16 are iEta -17 to -1
0089 // 17 to 33 are iEta 1 to 17
0090 static constexpr int toweriEta_fromAbsoluteID_shift = 16;
0091 
0092 // absolue IDs range from 0-71.
0093 // To align with detector tower IDs (1 - n_towers_Phi)
0094 // shift all indices by 37 and loop over after 72
0095 static constexpr int toweriPhi_fromAbsoluteID_shift_lowerHalf = 37;
0096 static constexpr int toweriPhi_fromAbsoluteID_shift_upperHalf = 35;
0097 
0098 float getEta_fromL2LinkCardTowerCrystal(int link, int card, int tower, int crystal) {
0099   int etaID = n_crystals_towerEta * (n_towers_per_link * ((link / n_links_card) % 2) + (tower % n_towers_per_link)) +
0100               crystal % n_crystals_towerEta;
0101   float size_cell = 2 * ECAL_eta_range / (n_crystals_towerEta * n_towers_Eta);
0102   return etaID * size_cell - ECAL_eta_range + half_crystal_size;
0103 }
0104 
0105 float getPhi_fromL2LinkCardTowerCrystal(int link, int card, int tower, int crystal) {
0106   int phiID = n_crystals_towerPhi * ((card * 24) + (n_links_card * (link / 8)) + (tower / n_towers_per_link)) +
0107               crystal / n_crystals_towerPhi;
0108   float size_cell = 2 * M_PI / (n_crystals_towerPhi * n_towers_Phi);
0109   return phiID * size_cell - M_PI + half_crystal_size;
0110 }
0111 
0112 int getCrystal_etaID(float eta) {
0113   float size_cell = 2 * ECAL_eta_range / (n_crystals_towerEta * n_towers_Eta);
0114   int etaID = int((eta + ECAL_eta_range) / size_cell);
0115   return etaID;
0116 }
0117 
0118 int convert_L2toL1_link(int link) { return link % n_links_card; }
0119 
0120 int convert_L2toL1_tower(int tower) { return tower; }
0121 
0122 int convert_L2toL1_card(int card, int link) { return card * n_clusters_4link + link / n_links_card; }
0123 
0124 int getCrystal_phiID(float phi) {
0125   float size_cell = 2 * M_PI / (n_crystals_towerPhi * n_towers_Phi);
0126   int phiID = int((phi + M_PI) / size_cell);
0127   return phiID;
0128 }
0129 
0130 int getTower_absoluteEtaID(float eta) {
0131   float size_cell = 2 * ECAL_eta_range / n_towers_Eta;
0132   int etaID = int((eta + ECAL_eta_range) / size_cell);
0133   return etaID;
0134 }
0135 
0136 int getTower_absolutePhiID(float phi) {
0137   float size_cell = 2 * M_PI / n_towers_Phi;
0138   int phiID = int((phi + M_PI) / size_cell);
0139   return phiID;
0140 }
0141 
0142 int getToweriEta_fromAbsoluteID(int id) {
0143   if (id < n_towers_per_link)
0144     return id - n_towers_per_link;
0145   else
0146     return id - toweriEta_fromAbsoluteID_shift;
0147 }
0148 
0149 int getToweriPhi_fromAbsoluteID(int id) {
0150   if (id < n_towers_Phi / 2)
0151     return id + toweriPhi_fromAbsoluteID_shift_lowerHalf;
0152   else
0153     return id - toweriPhi_fromAbsoluteID_shift_upperHalf;
0154 }
0155 
0156 float getTowerEta_fromAbsoluteID(int id) {
0157   float size_cell = 2 * ECAL_eta_range / n_towers_Eta;
0158   float eta = (id * size_cell) - ECAL_eta_range + 0.5 * size_cell;
0159   return eta;
0160 }
0161 
0162 float getTowerPhi_fromAbsoluteID(int id) {
0163   float size_cell = 2 * M_PI / n_towers_Phi;
0164   float phi = (id * size_cell) - M_PI + 0.5 * size_cell;
0165   return phi;
0166 }
0167 
0168 int getCrystalIDInTower(int etaID, int phiID) {
0169   return int(n_crystals_towerPhi * (phiID % n_crystals_towerPhi) + (etaID % n_crystals_towerEta));
0170 }
0171 
0172 int get_towerEta_fromCardTowerInCard(int card, int towerincard) {
0173   return n_towers_per_link * (card % 2) + towerincard % n_towers_per_link;
0174 }
0175 
0176 int get_towerPhi_fromCardTowerInCard(int card, int towerincard) {
0177   return 4 * (card / 2) + towerincard / n_towers_per_link;
0178 }
0179 
0180 int get_towerEta_fromCardLinkTower(int card, int link, int tower) { return n_towers_per_link * (card % 2) + tower; }
0181 
0182 int get_towerPhi_fromCardLinkTower(int card, int link, int tower) { return 4 * (card / 2) + link; }
0183 
0184 int getTowerID(int etaID, int phiID) {
0185   return int(n_towers_per_link * ((phiID / n_crystals_towerPhi) % 4) +
0186              (etaID / n_crystals_towerEta) % n_towers_per_link);
0187 }
0188 
0189 int getTower_phiID(int cluster_phiID) {  // Tower ID in card given crystal ID in total detector
0190   return int((cluster_phiID / n_crystals_towerPhi) % 4);
0191 }
0192 
0193 int getTower_etaID(int cluster_etaID) {  // Tower ID in card given crystal ID in total detector
0194   return int((cluster_etaID / n_crystals_towerEta) % n_towers_per_link);
0195 }
0196 
0197 int getEtaMax_card(int card) {
0198   int etamax = 0;
0199   if (card % 2 == 0)
0200     etamax = n_towers_per_link * n_crystals_towerEta - 1;  // First eta half. 5 crystals in eta in 1 tower.
0201   else
0202     etamax = n_towers_Eta * n_crystals_towerEta - 1;
0203   return etamax;
0204 }
0205 
0206 int getEtaMin_card(int card) {
0207   int etamin = 0;
0208   if (card % 2 == 0)
0209     etamin = 0 * n_crystals_towerEta;  // First eta half. 5 crystals in eta in 1 tower.
0210   else
0211     etamin = n_towers_per_link * n_crystals_towerEta;
0212   return etamin;
0213 }
0214 
0215 int getPhiMax_card(int card) {
0216   int phimax = ((card / 2) + 1) * 4 * n_crystals_towerPhi - 1;
0217   return phimax;
0218 }
0219 
0220 int getPhiMin_card(int card) {
0221   int phimin = (card / 2) * 4 * n_crystals_towerPhi;
0222   return phimin;
0223 }
0224 
0225 /* 
0226  * Replace in-line region boundary arithmetic with function that accounts for region 0 in negative eta cards
0227  * In the indexing convention of the old emulator,  region 0 is the region overlapping with the endcap, and is
0228  * only two towers wide in eta.
0229  */
0230 int getEtaMin_region(int card, int nregion) {
0231   // Special handling for negative-eta cards
0232   if (card % 2 == 0) {
0233     if (nregion == 0) {
0234       return getEtaMin_card(card);
0235     } else {
0236       return getEtaMin_card(card) + n_crystals_2towers + n_crystals_3towers * (nregion - 1);
0237     }
0238   }
0239   // Positive-eta cards: same as original in-line arithmetic
0240   else {
0241     return getEtaMin_card(card) + n_crystals_3towers * nregion;
0242   }
0243 }
0244 
0245 /* 
0246  * Replace in-line region boundary arithmetic that accounts for region 0 in negative eta cards.
0247  * Same as above but for max eta of the region. 
0248  */
0249 int getEtaMax_region(int card, int nregion) {
0250   // Special handling for negative-eta cards
0251   if (card % 2 == 0) {
0252     if (nregion == 0) {
0253       return getEtaMin_card(card) + n_crystals_2towers;
0254     } else {
0255       return getEtaMin_card(card) + n_crystals_2towers + (n_crystals_3towers * nregion);
0256     }
0257   }
0258   // Positive-eta cards: same as original in-line arithmetic
0259   else {
0260     return getEtaMin_card(card) + n_crystals_3towers * (nregion + 1);
0261   }
0262 }
0263 
0264 class L1EGCrystalClusterEmulatorProducer : public edm::stream::EDProducer<> {
0265 public:
0266   explicit L1EGCrystalClusterEmulatorProducer(const edm::ParameterSet&);
0267   ~L1EGCrystalClusterEmulatorProducer() override;
0268 
0269 private:
0270   void produce(edm::Event&, const edm::EventSetup&) override;
0271   bool passes_ss(float pt, float ss);
0272   bool passes_photon(float pt, float pss);
0273   bool passes_iso(float pt, float iso);
0274   bool passes_looseTkss(float pt, float ss);
0275   bool passes_looseTkiso(float pt, float iso);
0276   float get_calibrate(float uncorr);
0277 
0278   edm::EDGetTokenT<EcalEBTrigPrimDigiCollection> ecalTPEBToken_;
0279   edm::EDGetTokenT<edm::SortedCollection<HcalTriggerPrimitiveDigi> > hcalTPToken_;
0280   edm::ESGetToken<CaloTPGTranscoder, CaloTPGRecord> decoderTag_;
0281 
0282   l1tp2::ParametricCalibration calib_;
0283 
0284   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryTag_;
0285   const CaloSubdetectorGeometry* ebGeometry;
0286   const CaloSubdetectorGeometry* hbGeometry;
0287   edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> hbTopologyTag_;
0288   const HcalTopology* hcTopology_;
0289 
0290   struct mycluster {
0291     float c2x2_;
0292     float c2x5_;
0293     float c5x5_;
0294     int cshowershape_;
0295     int cshowershapeloosetk_;
0296     float cvalueshowershape_;
0297     int cphotonshowershape_;
0298     float cpt;   // ECAL pt
0299     int cbrem_;  // if brem corrections were applied
0300     float cWeightedEta_;
0301     float cWeightedPhi_;
0302     float ciso_;      // pt of cluster divided by 7x7 ECAL towers
0303     float crawIso_;   // raw isolation sum
0304     float chovere_;   // 5x5 HCAL towers divided by the ECAL cluster pt
0305     float craweta_;   // coordinates between -1.44 and 1.44
0306     float crawphi_;   // coordinates between -pi and pi
0307     float chcal_;     // 5x5 HCAL towers
0308     float ceta_;      // eta ID in the whole detector (between 0 and 5*34-1)
0309     float cphi_;      // phi ID in the whole detector (between 0 and 5*72-1)
0310     int ccrystalid_;  // crystal ID inside tower (between 0 and 24)
0311     int cinsidecrystalid_;
0312     int ctowerid_;  // tower ID inside card (between 0 and 4*n_towers_per_link-1)
0313   };
0314 
0315   class SimpleCaloHit {
0316   private:
0317     float pt_ = 0;
0318     float energy_ = 0.;
0319     bool isEndcapHit_ = false;  // If using endcap, we won't be using integer crystal indices
0320     bool stale_ = false;        // Hits become stale once used in clustering algorithm to prevent overlap in clusters
0321     bool used_ = false;
0322     GlobalVector position_;  // As opposed to GlobalPoint, so we can add them (for weighted average)
0323     HcalDetId id_hcal_;
0324     EBDetId id_;
0325 
0326   public:
0327     // tool functions
0328     inline void setPt() { pt_ = (position_.mag2() > 0) ? energy_ * sin(position_.theta()) : 0; };
0329     inline void setEnergy(float et) { energy_ = et / sin(position_.theta()); };
0330     inline void setIsEndcapHit(bool isEC) { isEndcapHit_ = isEC; };
0331     inline void setUsed(bool isUsed) { used_ = isUsed; };
0332     inline void setPosition(const GlobalVector& pos) { position_ = pos; };
0333     inline void setIdHcal(const HcalDetId& idhcal) { id_hcal_ = idhcal; };
0334     inline void setId(const EBDetId& id) { id_ = id; };
0335 
0336     inline float pt() const { return pt_; };
0337     inline float energy() const { return energy_; };
0338     inline bool isEndcapHit() const { return isEndcapHit_; };
0339     inline bool used() const { return used_; };
0340     inline const GlobalVector& position() const { return position_; };
0341     inline const EBDetId& id() const { return id_; };
0342 
0343     /*
0344      * Check if it falls within the boundary of a card.
0345      */
0346     bool isInCard(int cc) const {
0347       return (getCrystal_phiID(position().phi()) <= getPhiMax_card(cc) &&
0348               getCrystal_phiID(position().phi()) >= getPhiMin_card(cc) &&
0349               getCrystal_etaID(position().eta()) <= getEtaMax_card(cc) &&
0350               getCrystal_etaID(position().eta()) >= getEtaMin_card(cc));
0351     };
0352 
0353     /* 
0354      * Check if it falls within the boundary card AND a region in the card.
0355      */
0356     bool isInCardAndRegion(int cc, int nregion) const {
0357       bool isInRegionEta = (getCrystal_etaID(position().eta()) < getEtaMax_region(cc, nregion) &&
0358                             getCrystal_etaID(position().eta()) >= getEtaMin_region(cc, nregion));
0359       return (isInCard(cc) && isInRegionEta);
0360     }
0361 
0362     // Comparison functions with other SimpleCaloHit instances
0363     inline float deta(SimpleCaloHit& other) const { return position_.eta() - other.position().eta(); };
0364     int dieta(SimpleCaloHit& other) const {
0365       if (isEndcapHit_ || other.isEndcapHit())
0366         return 9999;  // We shouldn't compare integer indices in endcap, the map is not linear
0367       if (id_.ieta() * other.id().ieta() > 0)
0368         return id_.ieta() - other.id().ieta();
0369       return id_.ieta() - other.id().ieta() - 1;
0370     };
0371     inline float dphi(SimpleCaloHit& other) const {
0372       return reco::deltaPhi(static_cast<float>(position_.phi()), static_cast<float>(other.position().phi()));
0373     };
0374     int diphi(SimpleCaloHit& other) const {
0375       if (isEndcapHit_ || other.isEndcapHit())
0376         return 9999;  // We shouldn't compare integer indices in endcap, the map is not linear
0377       // Logic from EBDetId::distancePhi() without the abs()
0378       static constexpr int PI = 180;
0379       int result = id().iphi() - other.id().iphi();
0380       while (result > PI)
0381         result -= 2 * PI;
0382       while (result <= -PI)
0383         result += 2 * PI;
0384       return result;
0385     };
0386     int dieta_byCrystalID(SimpleCaloHit& other) const {
0387       return (getCrystal_etaID(position_.eta()) - getCrystal_etaID(other.position().eta()));
0388     };
0389     int diphi_byCrystalID(SimpleCaloHit& other) const {
0390       return (getCrystal_phiID(position_.phi()) - getCrystal_phiID(other.position().phi()));
0391     };
0392     int id_iEta() { return id_.ieta(); }
0393     int id_iPhi() { return id_.iphi(); }
0394     float distanceTo(SimpleCaloHit& other) const {
0395       // Treat position as a point, measure 3D distance
0396       // This is used for endcap hits, where we don't have a rectangular mapping
0397       return (position() - other.position()).mag();
0398     };
0399     bool operator==(SimpleCaloHit& other) const {
0400       return (id_ == other.id() && position() == other.position() && energy_ == other.energy() &&
0401               isEndcapHit_ == other.isEndcapHit());
0402     };
0403   };
0404 };
0405 
0406 L1EGCrystalClusterEmulatorProducer::L1EGCrystalClusterEmulatorProducer(const edm::ParameterSet& iConfig)
0407     : ecalTPEBToken_(consumes<EcalEBTrigPrimDigiCollection>(iConfig.getParameter<edm::InputTag>("ecalTPEB"))),
0408       hcalTPToken_(
0409           consumes<edm::SortedCollection<HcalTriggerPrimitiveDigi> >(iConfig.getParameter<edm::InputTag>("hcalTP"))),
0410       decoderTag_(esConsumes<CaloTPGTranscoder, CaloTPGRecord>(edm::ESInputTag("", ""))),
0411       calib_(iConfig.getParameter<edm::ParameterSet>("calib")),
0412       caloGeometryTag_(esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag("", ""))),
0413       hbTopologyTag_(esConsumes<HcalTopology, HcalRecNumberingRecord>(edm::ESInputTag("", ""))) {
0414   produces<l1tp2::CaloCrystalClusterCollection>();
0415   produces<BXVector<l1t::EGamma> >();
0416   produces<l1tp2::CaloTowerCollection>("L1CaloTowerCollection");
0417 }
0418 
0419 L1EGCrystalClusterEmulatorProducer::~L1EGCrystalClusterEmulatorProducer() {}
0420 
0421 void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0422   using namespace edm;
0423 
0424   edm::Handle<EcalEBTrigPrimDigiCollection> pcalohits;
0425   iEvent.getByToken(ecalTPEBToken_, pcalohits);
0426 
0427   const auto& caloGeometry = iSetup.getData(caloGeometryTag_);
0428   ebGeometry = caloGeometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0429   hbGeometry = caloGeometry.getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0430   const auto& hbTopology = iSetup.getData(hbTopologyTag_);
0431   hcTopology_ = &hbTopology;
0432   HcalTrigTowerGeometry theTrigTowerGeometry(hcTopology_);
0433 
0434   const auto& decoder = iSetup.getData(decoderTag_);
0435 
0436   //****************************************************************
0437   //******************* Get all the hits ***************************
0438   //****************************************************************
0439 
0440   // Get all the ECAL hits
0441   iEvent.getByToken(ecalTPEBToken_, pcalohits);
0442   std::vector<SimpleCaloHit> ecalhits;
0443 
0444   for (const auto& hit : *pcalohits.product()) {
0445     if (hit.encodedEt() > 0)  // hit.encodedEt() returns an int corresponding to 2x the crystal Et
0446     {
0447       // Et is 10 bit, by keeping the ADC saturation Et at 120 GeV it means that you have to divide by 8
0448       float et = hit.encodedEt() / 8.;
0449       if (et < cut_500_MeV)
0450         continue;  // keep the 500 MeV ET Cut
0451 
0452       auto cell = ebGeometry->getGeometry(hit.id());
0453 
0454       SimpleCaloHit ehit;
0455       ehit.setId(hit.id());
0456       ehit.setPosition(GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z()));
0457       ehit.setEnergy(et);
0458       ehit.setPt();
0459       ecalhits.push_back(ehit);
0460     }
0461   }
0462 
0463   // Get all the HCAL hits
0464   std::vector<SimpleCaloHit> hcalhits;
0465   edm::Handle<edm::SortedCollection<HcalTriggerPrimitiveDigi> > hbhecoll;
0466   iEvent.getByToken(hcalTPToken_, hbhecoll);
0467   for (const auto& hit : *hbhecoll.product()) {
0468     float et = decoder.hcaletValue(hit.id(), hit.t0());
0469     if (et <= 0)
0470       continue;
0471     if (!(hcTopology_->validHT(hit.id()))) {
0472       LogError("L1EGCrystalClusterEmulatorProducer")
0473           << " -- Hcal hit DetID not present in HCAL Geom: " << hit.id() << std::endl;
0474       throw cms::Exception("L1EGCrystalClusterEmulatorProducer");
0475       continue;
0476     }
0477     const std::vector<HcalDetId>& hcId = theTrigTowerGeometry.detIds(hit.id());
0478     if (hcId.empty()) {
0479       LogError("L1EGCrystalClusterEmulatorProducer")
0480           << "Cannot find any HCalDetId corresponding to " << hit.id() << std::endl;
0481       throw cms::Exception("L1EGCrystalClusterEmulatorProducer");
0482       continue;
0483     }
0484     if (hcId[0].subdetId() > 1)
0485       continue;
0486     GlobalVector hcal_tp_position = GlobalVector(0., 0., 0.);
0487     for (const auto& hcId_i : hcId) {
0488       if (hcId_i.subdetId() > 1)
0489         continue;
0490       auto cell = hbGeometry->getGeometry(hcId_i);
0491       if (cell == nullptr)
0492         continue;
0493       GlobalVector tmpVector = GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z());
0494       hcal_tp_position = tmpVector;
0495       break;
0496     }
0497     SimpleCaloHit hhit;
0498     hhit.setId(hit.id());
0499     hhit.setIdHcal(hit.id());
0500     hhit.setPosition(hcal_tp_position);
0501     hhit.setEnergy(et);
0502     hhit.setPt();
0503     hcalhits.push_back(hhit);
0504   }
0505 
0506   //*******************************************************************
0507   //********************** Do layer 1 *********************************
0508   //*******************************************************************
0509 
0510   // Definition of L1 outputs
0511   // 36 L1 cards send each 4 links with 17 towers
0512   float ECAL_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
0513   float HCAL_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
0514   int iEta_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
0515   int iPhi_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
0516   // 36 L1 cards send each 4 links with 3 clusters
0517   float energy_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0518   // 36 L1 cards send each 4 links with 3 clusters
0519   int brem_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0520   int towerID_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0521   int crystalID_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0522   int showerShape_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0523   int showerShapeLooseTk_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0524   int photonShowerShape_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0525 
0526   for (int ii = 0; ii < n_links_card; ++ii) {
0527     for (int jj = 0; jj < n_towers_per_link; ++jj) {
0528       for (int ll = 0; ll < n_towers_halfPhi; ++ll) {
0529         ECAL_tower_L1Card[ii][jj][ll] = 0;
0530         HCAL_tower_L1Card[ii][jj][ll] = 0;
0531         iPhi_tower_L1Card[ii][jj][ll] = -999;
0532         iEta_tower_L1Card[ii][jj][ll] = -999;
0533       }
0534     }
0535   }
0536   for (int ii = 0; ii < n_links_card; ++ii) {
0537     for (int jj = 0; jj < n_clusters_link; ++jj) {
0538       for (int ll = 0; ll < n_towers_halfPhi; ++ll) {
0539         energy_cluster_L1Card[ii][jj][ll] = 0;
0540         brem_cluster_L1Card[ii][jj][ll] = 0;
0541         towerID_cluster_L1Card[ii][jj][ll] = 0;
0542         crystalID_cluster_L1Card[ii][jj][ll] = 0;
0543       }
0544     }
0545   }
0546 
0547   // There is one list of clusters per card. We take the 12 highest pt per card
0548   vector<mycluster> cluster_list[n_towers_halfPhi];
0549   // After merging the clusters in different regions of a single L1 card
0550   vector<mycluster> cluster_list_merged[n_towers_halfPhi];
0551 
0552   for (int cc = 0; cc < n_towers_halfPhi; ++cc) {  // Loop over 36 L1 cards
0553     // Loop over 3x4 etaxphi regions to search for max 5 clusters
0554     for (int nregion = 0; nregion <= n_clusters_max; ++nregion) {
0555       int nclusters = 0;
0556       bool build_cluster = true;
0557 
0558       // Continue until 5 clusters have been built or there is no cluster left
0559       while (nclusters < n_clusters_max && build_cluster) {
0560         build_cluster = false;
0561         SimpleCaloHit centerhit;
0562 
0563         for (const auto& hit : ecalhits) {
0564           // Highest hit in good region with pt>1 and not used in any other cluster
0565           if (hit.isInCardAndRegion(cc, nregion) && !hit.used() && hit.pt() >= 1.0 && hit.pt() > centerhit.pt()) {
0566             centerhit = hit;
0567             build_cluster = true;
0568           }
0569         }
0570         if (build_cluster)
0571           nclusters++;
0572 
0573         // Use only the 5 most energetic clusters
0574         if (build_cluster && nclusters > 0 && nclusters <= n_clusters_max) {
0575           mycluster mc1;
0576           mc1.cpt = 0.0;
0577           mc1.cWeightedEta_ = 0.0;
0578           mc1.cWeightedPhi_ = 0.0;
0579           float leftlobe = 0;
0580           float rightlobe = 0;
0581           float e5x5 = 0;
0582           float e2x5_1 = 0;
0583           float e2x5_2 = 0;
0584           float e2x2_1 = 0;
0585           float e2x2_2 = 0;
0586           float e2x2_3 = 0;
0587           float e2x2_4 = 0;
0588           for (auto& hit : ecalhits) {
0589             if (hit.isInCardAndRegion(cc, nregion) && (hit.pt() > 0)) {
0590               if (abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) > 2 && hit.diphi(centerhit) <= 7) {
0591                 rightlobe += hit.pt();
0592               }
0593               if (abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) < -2 && hit.diphi(centerhit) >= -7) {
0594                 leftlobe += hit.pt();
0595               }
0596               if (abs(hit.dieta(centerhit)) <= 2 && abs(hit.diphi(centerhit)) <= 2) {
0597                 e5x5 += hit.energy();
0598               }
0599               if ((hit.dieta(centerhit) == 1 or hit.dieta(centerhit) == 0) &&
0600                   (hit.diphi(centerhit) == 1 or hit.diphi(centerhit) == 0)) {
0601                 e2x2_1 += hit.energy();
0602               }
0603               if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) &&
0604                   (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == 1)) {
0605                 e2x2_2 += hit.energy();
0606               }
0607               if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == 1) &&
0608                   (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == -1)) {
0609                 e2x2_3 += hit.energy();
0610               }
0611               if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) &&
0612                   (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == -1)) {
0613                 e2x2_4 += hit.energy();
0614               }
0615               if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == 1) && abs(hit.diphi(centerhit)) <= 2) {
0616                 e2x5_1 += hit.energy();
0617               }
0618               if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) && abs(hit.diphi(centerhit)) <= 2) {
0619                 e2x5_2 += hit.energy();
0620               }
0621             }
0622             if (hit.isInCardAndRegion(cc, nregion) && !hit.used() && hit.pt() > 0 && abs(hit.dieta(centerhit)) <= 1 &&
0623                 abs(hit.diphi(centerhit)) <= 2) {
0624               // clusters 3x5 in etaxphi using only the hits in the corresponding card and in the corresponding 3x4 region
0625               hit.setUsed(true);
0626               mc1.cpt += hit.pt();
0627               mc1.cWeightedEta_ += float(hit.pt()) * float(hit.position().eta());
0628               mc1.cWeightedPhi_ = mc1.cWeightedPhi_ + (float(hit.pt()) * float(hit.position().phi()));
0629             }
0630           }
0631           if (do_brem && (rightlobe > 0.10 * mc1.cpt or leftlobe > 0.10 * mc1.cpt)) {
0632             for (auto& hit : ecalhits) {
0633               if (hit.isInCardAndRegion(cc, nregion) && hit.pt() > 0 && !hit.used()) {
0634                 if (rightlobe > 0.10 * mc1.cpt && (leftlobe < 0.10 * mc1.cpt or rightlobe > leftlobe) &&
0635                     abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) > 2 && hit.diphi(centerhit) <= 7) {
0636                   mc1.cpt += hit.pt();
0637                   hit.setUsed(true);
0638                   mc1.cbrem_ = 1;
0639                 }
0640                 if (leftlobe > 0.10 * mc1.cpt && (rightlobe < 0.10 * mc1.cpt or leftlobe >= rightlobe) &&
0641                     abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) < -2 && hit.diphi(centerhit) >= -7) {
0642                   mc1.cpt += hit.pt();
0643                   hit.setUsed(true);
0644                   mc1.cbrem_ = 1;
0645                 }
0646               }
0647             }
0648           }
0649           mc1.c5x5_ = e5x5;
0650           mc1.c2x5_ = max(e2x5_1, e2x5_2);
0651           mc1.c2x2_ = e2x2_1;
0652           if (e2x2_2 > mc1.c2x2_)
0653             mc1.c2x2_ = e2x2_2;
0654           if (e2x2_3 > mc1.c2x2_)
0655             mc1.c2x2_ = e2x2_3;
0656           if (e2x2_4 > mc1.c2x2_)
0657             mc1.c2x2_ = e2x2_4;
0658           mc1.cWeightedEta_ = mc1.cWeightedEta_ / mc1.cpt;
0659           mc1.cWeightedPhi_ = mc1.cWeightedPhi_ / mc1.cpt;
0660           mc1.ceta_ = getCrystal_etaID(centerhit.position().eta());
0661           mc1.cphi_ = getCrystal_phiID(centerhit.position().phi());
0662           mc1.crawphi_ = centerhit.position().phi();
0663           mc1.craweta_ = centerhit.position().eta();
0664           cluster_list[cc].push_back(mc1);
0665         }  // End if 5 clusters per region
0666       }  // End while to find the 5 clusters
0667     }  // End loop over regions to search for clusters
0668     std::sort(begin(cluster_list[cc]), end(cluster_list[cc]), [](mycluster a, mycluster b) { return a.cpt > b.cpt; });
0669 
0670     // Merge clusters from different regions
0671     for (unsigned int jj = 0; jj < unsigned(cluster_list[cc].size()); ++jj) {
0672       for (unsigned int kk = jj + 1; kk < unsigned(cluster_list[cc].size()); ++kk) {
0673         if (std::abs(cluster_list[cc][jj].ceta_ - cluster_list[cc][kk].ceta_) < 2 &&
0674             std::abs(cluster_list[cc][jj].cphi_ - cluster_list[cc][kk].cphi_) < 2) {  // Diagonal + exact neighbors
0675           if (cluster_list[cc][kk].cpt > cluster_list[cc][jj].cpt) {
0676             cluster_list[cc][kk].cpt += cluster_list[cc][jj].cpt;
0677             cluster_list[cc][kk].c5x5_ += cluster_list[cc][jj].c5x5_;
0678             cluster_list[cc][kk].c2x5_ += cluster_list[cc][jj].c2x5_;
0679             cluster_list[cc][jj].cpt = 0;
0680             cluster_list[cc][jj].c5x5_ = 0;
0681             cluster_list[cc][jj].c2x5_ = 0;
0682             cluster_list[cc][jj].c2x2_ = 0;
0683           } else {
0684             cluster_list[cc][jj].cpt += cluster_list[cc][kk].cpt;
0685             cluster_list[cc][jj].c5x5_ += cluster_list[cc][kk].c5x5_;
0686             cluster_list[cc][jj].c2x5_ += cluster_list[cc][kk].c2x5_;
0687             cluster_list[cc][kk].cpt = 0;
0688             cluster_list[cc][kk].c2x2_ = 0;
0689             cluster_list[cc][kk].c2x5_ = 0;
0690             cluster_list[cc][kk].c5x5_ = 0;
0691           }
0692         }
0693       }
0694       if (cluster_list[cc][jj].cpt > 0) {
0695         cluster_list[cc][jj].cpt =
0696             cluster_list[cc][jj].cpt *
0697             calib_(cluster_list[cc][jj].cpt,
0698                    std::abs(cluster_list[cc][jj].craweta_));  //Mark's calibration as a function of eta and pt
0699         cluster_list_merged[cc].push_back(cluster_list[cc][jj]);
0700       }
0701     }
0702     std::sort(begin(cluster_list_merged[cc]), end(cluster_list_merged[cc]), [](mycluster a, mycluster b) {
0703       return a.cpt > b.cpt;
0704     });
0705 
0706     // Fill cluster information in the arrays. We keep max 12 clusters (distributed between 4 links)
0707     for (unsigned int jj = 0; jj < unsigned(cluster_list_merged[cc].size()) && jj < n_clusters_4link; ++jj) {
0708       crystalID_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] =
0709           getCrystalIDInTower(cluster_list_merged[cc][jj].ceta_, cluster_list_merged[cc][jj].cphi_);
0710       towerID_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] =
0711           getTowerID(cluster_list_merged[cc][jj].ceta_, cluster_list_merged[cc][jj].cphi_);
0712       energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = cluster_list_merged[cc][jj].cpt;
0713       brem_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = cluster_list_merged[cc][jj].cbrem_;
0714       if (passes_ss(cluster_list_merged[cc][jj].cpt,
0715                     cluster_list_merged[cc][jj].c2x5_ / cluster_list_merged[cc][jj].c5x5_))
0716         showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1;
0717       else
0718         showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0;
0719       if (passes_looseTkss(cluster_list_merged[cc][jj].cpt,
0720                            cluster_list_merged[cc][jj].c2x5_ / cluster_list_merged[cc][jj].c5x5_))
0721         showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1;
0722       else
0723         showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0;
0724       if (passes_photon(cluster_list_merged[cc][jj].cpt,
0725                         cluster_list_merged[cc][jj].c2x2_ / cluster_list_merged[cc][jj].c2x5_))
0726         photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1;
0727       else
0728         photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0;
0729     }
0730 
0731     // Loop over calo ecal hits to get the ECAL towers. Take only hits that have not been used to make clusters
0732     for (const auto& hit : ecalhits) {
0733       if (hit.isInCard(cc) && !hit.used()) {
0734         for (int jj = 0; jj < n_links_card; ++jj) {                                        // loop over 4 links per card
0735           if ((getCrystal_phiID(hit.position().phi()) / n_crystals_towerPhi) % 4 == jj) {  // Go to ID tower modulo 4
0736             for (int ii = 0; ii < n_towers_per_link; ++ii) {
0737               // Apply Mark's calibration at the same time (row of the lowest pT, as a function of eta)
0738               if ((getCrystal_etaID(hit.position().eta()) / n_crystals_towerEta) % n_towers_per_link == ii) {
0739                 ECAL_tower_L1Card[jj][ii][cc] += hit.pt() * calib_(0, std::abs(hit.position().eta()));
0740                 iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(hit.position().eta());
0741                 iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(hit.position().phi());
0742               }
0743             }  // end of loop over eta towers
0744           }
0745         }  // end of loop over phi links
0746         // Make sure towers with 0 ET are initialized with proper iEta, iPhi coordinates
0747         static constexpr float tower_width = 0.0873;
0748         for (int jj = 0; jj < n_links_card; ++jj) {
0749           for (int ii = 0; ii < n_towers_per_link; ++ii) {
0750             float phi = getPhiMin_card(cc) * tower_width / n_crystals_towerPhi - M_PI + (jj + 0.5) * tower_width;
0751             float eta = getEtaMin_card(cc) * tower_width / n_crystals_towerEta - n_towers_per_link * tower_width +
0752                         (ii + 0.5) * tower_width;
0753             iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(eta);
0754             iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(phi);
0755           }
0756         }
0757       }  // end of check if inside card
0758     }  // end of loop over hits to build towers
0759 
0760     // Loop over hcal hits to get the HCAL towers.
0761     for (const auto& hit : hcalhits) {
0762       if (hit.isInCard(cc) && hit.pt() > 0) {
0763         for (int jj = 0; jj < n_links_card; ++jj) {
0764           if ((getCrystal_phiID(hit.position().phi()) / n_crystals_towerPhi) % n_links_card == jj) {
0765             for (int ii = 0; ii < n_towers_per_link; ++ii) {
0766               if ((getCrystal_etaID(hit.position().eta()) / n_crystals_towerEta) % n_towers_per_link == ii) {
0767                 HCAL_tower_L1Card[jj][ii][cc] += hit.pt();
0768                 iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(hit.position().eta());
0769                 iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(hit.position().phi());
0770               }
0771             }  // end of loop over eta towers
0772           }
0773         }  // end of loop over phi links
0774       }  // end of check if inside card
0775     }  // end of loop over hits to build towers
0776 
0777     // Give back energy of not used clusters to the towers (if there are more than 12 clusters)
0778     for (unsigned int kk = n_clusters_4link; kk < cluster_list_merged[cc].size(); ++kk) {
0779       if (cluster_list_merged[cc][kk].cpt > 0) {
0780         ECAL_tower_L1Card[getTower_phiID(cluster_list_merged[cc][kk].cphi_)]
0781                          [getTower_etaID(cluster_list_merged[cc][kk].ceta_)][cc] += cluster_list_merged[cc][kk].cpt;
0782       }
0783     }
0784   }  //End of loop over cards
0785 
0786   //*********************************************************
0787   //******************** Do layer 2 *************************
0788   //*********************************************************
0789 
0790   // Definition of L2 outputs
0791   float HCAL_tower_L2Card[n_links_GCTcard][n_towers_per_link]
0792                          [n_GCTcards];  // 3 L2 cards send each 48 links with 17 towers
0793   float ECAL_tower_L2Card[n_links_GCTcard][n_towers_per_link][n_GCTcards];
0794   int iEta_tower_L2Card[n_links_GCTcard][n_towers_per_link][n_GCTcards];
0795   int iPhi_tower_L2Card[n_links_GCTcard][n_towers_per_link][n_GCTcards];
0796   float energy_cluster_L2Card[n_links_GCTcard][n_clusters_per_link]
0797                              [n_GCTcards];  // 3 L2 cards send each 48 links with 2 clusters
0798   float brem_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0799   int towerID_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0800   int crystalID_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0801   float isolation_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0802   float HE_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0803   int showerShape_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0804   int showerShapeLooseTk_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0805   int photonShowerShape_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0806 
0807   for (int ii = 0; ii < n_links_GCTcard; ++ii) {
0808     for (int jj = 0; jj < n_towers_per_link; ++jj) {
0809       for (int ll = 0; ll < n_GCTcards; ++ll) {
0810         HCAL_tower_L2Card[ii][jj][ll] = 0;
0811         ECAL_tower_L2Card[ii][jj][ll] = 0;
0812         iEta_tower_L2Card[ii][jj][ll] = 0;
0813         iPhi_tower_L2Card[ii][jj][ll] = 0;
0814       }
0815     }
0816   }
0817   for (int ii = 0; ii < n_links_GCTcard; ++ii) {
0818     for (int jj = 0; jj < n_clusters_per_link; ++jj) {
0819       for (int ll = 0; ll < n_GCTcards; ++ll) {
0820         energy_cluster_L2Card[ii][jj][ll] = 0;
0821         brem_cluster_L2Card[ii][jj][ll] = 0;
0822         towerID_cluster_L2Card[ii][jj][ll] = 0;
0823         crystalID_cluster_L2Card[ii][jj][ll] = 0;
0824         isolation_cluster_L2Card[ii][jj][ll] = 0;
0825         HE_cluster_L2Card[ii][jj][ll] = 0;
0826         photonShowerShape_cluster_L2Card[ii][jj][ll] = 0;
0827         showerShape_cluster_L2Card[ii][jj][ll] = 0;
0828         showerShapeLooseTk_cluster_L2Card[ii][jj][ll] = 0;
0829       }
0830     }
0831   }
0832 
0833   // There is one list of clusters per equivalent of L1 card. We take the 8 highest pt.
0834   vector<mycluster> cluster_list_L2[n_towers_halfPhi];
0835 
0836   // Merge clusters on the phi edges
0837   for (int ii = 0; ii < n_borders_phi; ++ii) {  // 18 borders in phi
0838     for (int jj = 0; jj < n_eta_bins; ++jj) {   // 2 eta bins
0839       int card_left = 2 * ii + jj;
0840       int card_right = 2 * ii + jj + 2;
0841       if (card_right > 35)
0842         card_right = card_right - n_towers_halfPhi;
0843       for (int kk = 0; kk < n_clusters_4link; ++kk) {  // 12 clusters in the first card. We check the right side
0844         if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 50 &&
0845             crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 19 &&
0846             energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 0) {
0847           for (int ll = 0; ll < n_clusters_4link; ++ll) {  // We check the 12 clusters in the card on the right
0848             if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_per_link &&
0849                 crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < 5 &&
0850                 std::abs(
0851                     5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) % n_towers_per_link +
0852                     crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] % 5 -
0853                     5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) % n_towers_per_link -
0854                     crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] % 5) < 2) {
0855               if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] >
0856                   energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) {
0857                 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] +=
0858                     energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right];
0859                 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 0;
0860               }  // The least energetic cluster is merged to the most energetic one
0861               else {
0862                 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] +=
0863                     energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left];
0864                 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 0;
0865               }
0866             }
0867           }
0868         }
0869       }
0870     }
0871   }
0872 
0873   // Bremsstrahlung corrections. Merge clusters on the phi edges depending on pT (pt more than 10 percent, dphi leq 5, deta leq 1)
0874   if (do_brem) {
0875     for (int ii = 0; ii < n_borders_phi; ++ii) {  // 18 borders in phi
0876       for (int jj = 0; jj < n_eta_bins; ++jj) {   // 2 eta bins
0877         int card_left = 2 * ii + jj;
0878         int card_right = 2 * ii + jj + 2;
0879         if (card_right > 35)
0880           card_right = card_right - n_towers_halfPhi;
0881         // 12 clusters in the first card. We check the right side crystalID_cluster_L1Card[kk%4][kk/4][card_left]
0882         for (int kk = 0; kk < n_clusters_4link; ++kk) {
0883           // if the tower is at the edge there might be brem corrections, whatever the crystal ID
0884           if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 50 &&
0885               energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 0) {
0886             for (int ll = 0; ll < n_clusters_4link; ++ll) {  // We check the 12 clusters in the card on the right
0887               //Distance of 1 max in eta
0888               if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_per_link &&
0889                   std::abs(5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) %
0890                                n_towers_per_link +
0891                            crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] % 5 -
0892                            5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) %
0893                                n_towers_per_link -
0894                            crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] % 5) <= 1) {
0895                 //Distance of 5 max in phi
0896                 if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_per_link &&
0897                     std::abs(5 + crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] / 5 -
0898                              (crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] / 5)) <= 5) {
0899                   if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] >
0900                           energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] &&
0901                       energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] >
0902                           0.10 * energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) {
0903                     energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] +=
0904                         energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right];
0905                     energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 0;
0906                     brem_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 1;
0907                   }
0908                   // The least energetic cluster is merged to the most energetic one, if its pt is at least ten percent
0909                   else if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_right] >=
0910                                energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_left] &&
0911                            energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_left] >
0912                                0.10 * energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_right]) {
0913                     energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] +=
0914                         energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left];
0915                     energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 0;
0916                     brem_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 1;
0917                   }
0918                 }  //max distance eta
0919               }  //max distance phi
0920             }  //max distance phi
0921           }
0922         }
0923       }
0924     }
0925   }
0926 
0927   // Merge clusters on the eta edges
0928   for (int ii = 0; ii < n_borders_eta; ++ii) {  // 18 borders in eta
0929     int card_bottom = 2 * ii;
0930     int card_top = 2 * ii + 1;
0931     for (int kk = 0; kk < n_clusters_4link; ++kk) {  // 12 clusters in the first card. We check the top side
0932       if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] % n_towers_per_link == 16 &&
0933           crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] % 5 == n_links_card &&
0934           energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] >
0935               0) {                                       // If there is one cluster on the right side of the first card
0936         for (int ll = 0; ll < n_clusters_4link; ++ll) {  // We check the card on the right
0937           if (std::abs(
0938                   5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] / n_towers_per_link) +
0939                   crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] / 5 -
0940                   5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] / n_towers_per_link) -
0941                   crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] / 5) < 2) {
0942             if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] >
0943                 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_bottom]) {
0944               energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] +=
0945                   energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top];
0946               energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] = 0;
0947             } else {
0948               energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] +=
0949                   energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom];
0950               energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] = 0;
0951             }
0952           }
0953         }
0954       }
0955     }
0956   }
0957 
0958   // Regroup the new clusters per equivalent of L1 card geometry
0959   for (int ii = 0; ii < n_towers_halfPhi; ++ii) {
0960     for (int jj = 0; jj < n_clusters_4link; ++jj) {
0961       if (energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii] > 0) {
0962         mycluster mc1;
0963         mc1.cpt = energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0964         mc1.cbrem_ = brem_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0965         mc1.ctowerid_ = towerID_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0966         mc1.ccrystalid_ = crystalID_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0967         mc1.cshowershape_ = showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0968         mc1.cshowershapeloosetk_ = showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0969         mc1.cphotonshowershape_ = photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0970         cluster_list_L2[ii].push_back(mc1);
0971       }
0972     }
0973     std::sort(
0974         begin(cluster_list_L2[ii]), end(cluster_list_L2[ii]), [](mycluster a, mycluster b) { return a.cpt > b.cpt; });
0975   }
0976 
0977   // If there are more than 8 clusters per equivalent of L1 card we need to put them back in the towers
0978   for (int ii = 0; ii < n_towers_halfPhi; ++ii) {
0979     for (unsigned int jj = n_clusters_per_L1card; jj < n_clusters_4link && jj < cluster_list_L2[ii].size(); ++jj) {
0980       if (cluster_list_L2[ii][jj].cpt > 0) {
0981         ECAL_tower_L1Card[cluster_list_L2[ii][jj].ctowerid_ / n_towers_per_link]
0982                          [cluster_list_L2[ii][jj].ctowerid_ % n_towers_per_link][ii] += cluster_list_L2[ii][jj].cpt;
0983         cluster_list_L2[ii][jj].cpt = 0;
0984         cluster_list_L2[ii][jj].ctowerid_ = 0;
0985         cluster_list_L2[ii][jj].ccrystalid_ = 0;
0986       }
0987     }
0988   }
0989 
0990   // Compute isolation (7*7 ECAL towers) and HCAL energy (5x5 HCAL towers)
0991   for (int ii = 0; ii < n_towers_halfPhi; ++ii) {  // Loop over the new cluster list (stored in 36x8 format)
0992     for (unsigned int jj = 0; jj < n_clusters_per_L1card && jj < cluster_list_L2[ii].size(); ++jj) {
0993       int cluster_etaOfTower_fullDetector = get_towerEta_fromCardTowerInCard(ii, cluster_list_L2[ii][jj].ctowerid_);
0994       int cluster_phiOfTower_fullDetector = get_towerPhi_fromCardTowerInCard(ii, cluster_list_L2[ii][jj].ctowerid_);
0995       float hcal_nrj = 0.0;
0996       float isolation = 0.0;
0997       int ntowers = 0;
0998 
0999       for (int iii = 0; iii < n_towers_halfPhi; ++iii) {  // The clusters have to be added to the isolation
1000         for (unsigned int jjj = 0; jjj < n_clusters_per_L1card && jjj < cluster_list_L2[iii].size(); ++jjj) {
1001           if (!(iii == ii && jjj == jj)) {
1002             int cluster2_eta = get_towerEta_fromCardTowerInCard(iii, cluster_list_L2[iii][jjj].ctowerid_);
1003             int cluster2_phi = get_towerPhi_fromCardTowerInCard(iii, cluster_list_L2[iii][jjj].ctowerid_);
1004             if (abs(cluster2_eta - cluster_etaOfTower_fullDetector) <= 2 &&
1005                 (abs(cluster2_phi - cluster_phiOfTower_fullDetector) <= 2 or
1006                  abs(cluster2_phi - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) {
1007               isolation += cluster_list_L2[iii][jjj].cpt;
1008             }
1009           }
1010         }
1011       }
1012       for (int kk = 0; kk < n_towers_halfPhi; ++kk) {       // 36 cards
1013         for (int ll = 0; ll < n_links_card; ++ll) {         // 4 links per card
1014           for (int mm = 0; mm < n_towers_per_link; ++mm) {  // 17 towers per link
1015             int etaOftower_fullDetector = get_towerEta_fromCardLinkTower(kk, ll, mm);
1016             int phiOftower_fullDetector = get_towerPhi_fromCardLinkTower(kk, ll, mm);
1017             // First do ECAL
1018             // The towers are within 3. Needs to stitch the two phi sides together
1019             if (abs(etaOftower_fullDetector - cluster_etaOfTower_fullDetector) <= 2 &&
1020                 (abs(phiOftower_fullDetector - cluster_phiOfTower_fullDetector) <= 2 or
1021                  abs(phiOftower_fullDetector - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) {
1022               // Remove the column outside of the L2 card:  values (0,71), (23,26), (24,21), (47,50), (48,50), (71,2)
1023               if (!((cluster_phiOfTower_fullDetector == 0 && phiOftower_fullDetector == 71) or
1024                     (cluster_phiOfTower_fullDetector == 23 && phiOftower_fullDetector == 26) or
1025                     (cluster_phiOfTower_fullDetector == 24 && phiOftower_fullDetector == 21) or
1026                     (cluster_phiOfTower_fullDetector == 47 && phiOftower_fullDetector == 50) or
1027                     (cluster_phiOfTower_fullDetector == 48 && phiOftower_fullDetector == 45) or
1028                     (cluster_phiOfTower_fullDetector == 71 && phiOftower_fullDetector == 2))) {
1029                 isolation += ECAL_tower_L1Card[ll][mm][kk];
1030                 ntowers++;
1031               }
1032             }
1033             // Now do HCAL
1034             // The towers are within 2. Needs to stitch the two phi sides together
1035             if (abs(etaOftower_fullDetector - cluster_etaOfTower_fullDetector) <= 2 &&
1036                 (abs(phiOftower_fullDetector - cluster_phiOfTower_fullDetector) <= 2 or
1037                  abs(phiOftower_fullDetector - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) {
1038               hcal_nrj += HCAL_tower_L1Card[ll][mm][kk];
1039             }
1040           }
1041         }
1042       }
1043       // If we summed over fewer than 5*5 = 25 towers (because the cluster was near the edge), scale up the isolation sum
1044       int nTowersIn5x5Window = 5 * 5;
1045       cluster_list_L2[ii][jj].ciso_ = ((isolation) * (nTowersIn5x5Window / ntowers)) / cluster_list_L2[ii][jj].cpt;
1046       cluster_list_L2[ii][jj].crawIso_ = ((isolation) * (nTowersIn5x5Window / ntowers));
1047       cluster_list_L2[ii][jj].chovere_ = hcal_nrj / cluster_list_L2[ii][jj].cpt;
1048     }
1049   }
1050 
1051   //Reformat the information inside the 3 L2 cards
1052   //First let's fill the towers
1053   for (int ii = 0; ii < n_links_GCTcard; ++ii) {
1054     for (int jj = 0; jj < n_towers_per_link; ++jj) {
1055       for (int ll = 0; ll < 3; ++ll) {
1056         ECAL_tower_L2Card[ii][jj][ll] =
1057             ECAL_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
1058         HCAL_tower_L2Card[ii][jj][ll] =
1059             HCAL_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
1060         iEta_tower_L2Card[ii][jj][ll] =
1061             iEta_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
1062         iPhi_tower_L2Card[ii][jj][ll] =
1063             iPhi_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
1064       }
1065     }
1066   }
1067 
1068   //Second let's fill the clusters
1069   for (int ii = 0; ii < n_towers_halfPhi; ++ii) {  // The cluster list is still in the L1 like geometry
1070     for (unsigned int jj = 0; jj < unsigned(cluster_list_L2[ii].size()) && jj < n_clusters_per_L1card; ++jj) {
1071       crystalID_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1072                               [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ccrystalid_;
1073       towerID_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1074                             [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ctowerid_;
1075       energy_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1076                            [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cpt;
1077       brem_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1078                          [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cbrem_;
1079       isolation_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1080                               [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ciso_;
1081       HE_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1082                        [ii / n_clusters_4link] = cluster_list_L2[ii][jj].chovere_;
1083       showerShape_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1084                                 [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cshowershape_;
1085       showerShapeLooseTk_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1086                                        [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cshowershapeloosetk_;
1087       photonShowerShape_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1088                                       [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cphotonshowershape_;
1089     }
1090   }
1091 
1092   auto L1EGXtalClusters = std::make_unique<l1tp2::CaloCrystalClusterCollection>();
1093   auto L1EGammas = std::make_unique<l1t::EGammaBxCollection>();
1094   auto L1CaloTowers = std::make_unique<l1tp2::CaloTowerCollection>();
1095 
1096   // Fill the cluster collection and towers as well
1097   for (int ii = 0; ii < n_links_GCTcard; ++ii) {  // 48 links
1098     for (int ll = 0; ll < n_GCTcards; ++ll) {     // 3 cards
1099       // For looping over the Towers a few lines below
1100       for (int jj = 0; jj < 2; ++jj) {  // 2 L1EGs
1101         if (energy_cluster_L2Card[ii][jj][ll] > 0.45) {
1102           reco::Candidate::PolarLorentzVector p4calibrated(
1103               energy_cluster_L2Card[ii][jj][ll],
1104               getEta_fromL2LinkCardTowerCrystal(
1105                   ii, ll, towerID_cluster_L2Card[ii][jj][ll], crystalID_cluster_L2Card[ii][jj][ll]),
1106               getPhi_fromL2LinkCardTowerCrystal(
1107                   ii, ll, towerID_cluster_L2Card[ii][jj][ll], crystalID_cluster_L2Card[ii][jj][ll]),
1108               0.);
1109           SimpleCaloHit centerhit;
1110           bool is_iso = passes_iso(energy_cluster_L2Card[ii][jj][ll], isolation_cluster_L2Card[ii][jj][ll]);
1111           bool is_looseTkiso =
1112               passes_looseTkiso(energy_cluster_L2Card[ii][jj][ll], isolation_cluster_L2Card[ii][jj][ll]);
1113           bool is_ss = (showerShape_cluster_L2Card[ii][jj][ll] == 1);
1114           bool is_photon = (photonShowerShape_cluster_L2Card[ii][jj][ll] == 1) && is_ss && is_iso;
1115           bool is_looseTkss = (showerShapeLooseTk_cluster_L2Card[ii][jj][ll] == 1);
1116           // All the ID set to Standalone WP! Some dummy values for non calculated variables
1117           l1tp2::CaloCrystalCluster cluster(p4calibrated,
1118                                             energy_cluster_L2Card[ii][jj][ll],
1119                                             HE_cluster_L2Card[ii][jj][ll],
1120                                             isolation_cluster_L2Card[ii][jj][ll],
1121                                             centerhit.id(),
1122                                             -1000,
1123                                             float(brem_cluster_L2Card[ii][jj][ll]),
1124                                             -1000,
1125                                             -1000,
1126                                             energy_cluster_L2Card[ii][jj][ll],
1127                                             -1,
1128                                             is_iso && is_ss,
1129                                             is_iso && is_ss,
1130                                             is_photon,
1131                                             is_iso && is_ss,
1132                                             is_looseTkiso && is_looseTkss,
1133                                             is_iso && is_ss);
1134           // Experimental parameters, don't want to bother with hardcoding them in data format
1135           std::map<std::string, float> params;
1136           params["standaloneWP_showerShape"] = is_ss;
1137           params["standaloneWP_isolation"] = is_iso;
1138           params["trkMatchWP_showerShape"] = is_looseTkss;
1139           params["trkMatchWP_isolation"] = is_looseTkiso;
1140           params["TTiEta"] = getToweriEta_fromAbsoluteID(getTower_absoluteEtaID(p4calibrated.eta()));
1141           params["TTiPhi"] = getToweriPhi_fromAbsoluteID(getTower_absolutePhiID(p4calibrated.phi()));
1142           cluster.setExperimentalParams(params);
1143           L1EGXtalClusters->push_back(cluster);
1144 
1145           int standaloneWP = (int)(is_iso && is_ss);
1146           int looseL1TkMatchWP = (int)(is_looseTkiso && is_looseTkss);
1147           int photonWP = (int)(is_photon);
1148           int quality =
1149               (standaloneWP * std::pow(2, 0)) + (looseL1TkMatchWP * std::pow(2, 1)) + (photonWP * std::pow(2, 2));
1150           L1EGammas->push_back(
1151               0, l1t::EGamma(p4calibrated, p4calibrated.pt(), p4calibrated.eta(), p4calibrated.phi(), quality, 1));
1152         }
1153       }
1154     }
1155   }
1156 
1157   for (int ii = 0; ii < n_links_GCTcard; ++ii) {  // 48 links
1158     for (int ll = 0; ll < n_GCTcards; ++ll) {     // 3 cards
1159       // Fill the tower collection
1160       for (int jj = 0; jj < n_towers_per_link; ++jj) {  // 17 TTs
1161         l1tp2::CaloTower l1CaloTower;
1162         l1CaloTower.setEcalTowerEt(ECAL_tower_L2Card[ii][jj][ll]);
1163         l1CaloTower.setHcalTowerEt(HCAL_tower_L2Card[ii][jj][ll]);
1164         l1CaloTower.setTowerIEta(getToweriEta_fromAbsoluteID(iEta_tower_L2Card[ii][jj][ll]));
1165         l1CaloTower.setTowerIPhi(getToweriPhi_fromAbsoluteID(iPhi_tower_L2Card[ii][jj][ll]));
1166         l1CaloTower.setTowerEta(getTowerEta_fromAbsoluteID(iEta_tower_L2Card[ii][jj][ll]));
1167         l1CaloTower.setTowerPhi(getTowerPhi_fromAbsoluteID(iPhi_tower_L2Card[ii][jj][ll]));
1168         // Some towers have incorrect eta/phi because that wasn't initialized in certain cases, fix these
1169         static float constexpr towerEtaUpperUnitialized = -80;
1170         static float constexpr towerPhiUpperUnitialized = -90;
1171         if (l1CaloTower.towerEta() < towerEtaUpperUnitialized && l1CaloTower.towerPhi() < towerPhiUpperUnitialized) {
1172           l1CaloTower.setTowerEta(l1t::CaloTools::towerEta(l1CaloTower.towerIEta()));
1173           l1CaloTower.setTowerPhi(l1t::CaloTools::towerPhi(l1CaloTower.towerIEta(), l1CaloTower.towerIPhi()));
1174         }
1175         l1CaloTower.setIsBarrel(true);
1176 
1177         // Add L1EGs if they match in iEta / iPhi
1178         // L1EGs are already pT ordered, we will take the ID info for the leading one, but pT as the sum
1179         for (const auto& l1eg : *L1EGXtalClusters) {
1180           if (l1eg.experimentalParam("TTiEta") != l1CaloTower.towerIEta())
1181             continue;
1182           if (l1eg.experimentalParam("TTiPhi") != l1CaloTower.towerIPhi())
1183             continue;
1184 
1185           int n_L1eg = l1CaloTower.nL1eg();
1186           l1CaloTower.setNL1eg(n_L1eg++);
1187           l1CaloTower.setL1egTowerEt(l1CaloTower.l1egTowerEt() + l1eg.pt());
1188           // Don't record L1EG quality info for subleading L1EG
1189           if (l1CaloTower.nL1eg() > 1)
1190             continue;
1191           l1CaloTower.setL1egTrkSS(l1eg.experimentalParam("trkMatchWP_showerShape"));
1192           l1CaloTower.setL1egTrkIso(l1eg.experimentalParam("trkMatchWP_isolation"));
1193           l1CaloTower.setL1egStandaloneSS(l1eg.experimentalParam("standaloneWP_showerShape"));
1194           l1CaloTower.setL1egStandaloneIso(l1eg.experimentalParam("standaloneWP_isolation"));
1195         }
1196 
1197         L1CaloTowers->push_back(l1CaloTower);
1198       }
1199     }
1200   }
1201 
1202   iEvent.put(std::move(L1EGXtalClusters));
1203   iEvent.put(std::move(L1EGammas));
1204   iEvent.put(std::move(L1CaloTowers), "L1CaloTowerCollection");
1205 }
1206 
1207 bool L1EGCrystalClusterEmulatorProducer::passes_iso(float pt, float iso) {
1208   bool is_iso = true;
1209   if (pt < slideIsoPtThreshold) {
1210     if (!((a0_80 - a1_80 * pt) > iso))
1211       is_iso = false;
1212   } else {
1213     if (iso > a0)
1214       is_iso = false;
1215   }
1216   if (pt > plateau_ss)
1217     is_iso = true;
1218   return is_iso;
1219 }
1220 
1221 bool L1EGCrystalClusterEmulatorProducer::passes_looseTkiso(float pt, float iso) {
1222   bool is_iso = (b0 + b1 * std::exp(-b2 * pt) > iso);
1223   if (pt > plateau_ss)
1224     is_iso = true;
1225   return is_iso;
1226 }
1227 
1228 bool L1EGCrystalClusterEmulatorProducer::passes_ss(float pt, float ss) {
1229   bool is_ss = ((c0_ss + c1_ss * std::exp(-c2_ss * pt)) <= ss);
1230   if (pt > plateau_ss)
1231     is_ss = true;
1232   return is_ss;
1233 }
1234 
1235 bool L1EGCrystalClusterEmulatorProducer::passes_photon(float pt, float pss) {
1236   bool is_ss = (pss > d0 - d1 * pt);
1237   if (pt > plateau_ss)
1238     is_ss = true;
1239   return is_ss;
1240 }
1241 
1242 bool L1EGCrystalClusterEmulatorProducer::passes_looseTkss(float pt, float ss) {
1243   bool is_ss = ((e0_looseTkss - e1_looseTkss * std::exp(-e2_looseTkss * pt)) <= ss);
1244   if (pt > plateau_ss)
1245     is_ss = true;
1246   return is_ss;
1247 }
1248 
1249 //define this as a plug-in
1250 DEFINE_FWK_MODULE(L1EGCrystalClusterEmulatorProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!