Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-05-07 00:41:28 

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