Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-03-26 01:51:22

0001 //------------------------------------
0002 // Helper functions for Phase2L1CaloEGammaEmulator.cc
0003 //------------------------------------
0004 #ifndef L1Trigger_L1CaloTrigger_Phase2L1CaloEGammaUtils
0005 #define L1Trigger_L1CaloTrigger_Phase2L1CaloEGammaUtils
0006 
0007 #include <ap_int.h>
0008 #include <array>
0009 #include <cstdio>
0010 #include <fstream>
0011 #include <iomanip>
0012 #include <iostream>
0013 
0014 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0015 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0016 #include "L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h"
0017 
0018 // Output collections
0019 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloCrystalCluster.h"
0020 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h"
0021 #include "DataFormats/L1TCalorimeterPhase2/interface/DigitizedClusterCorrelator.h"
0022 #include "DataFormats/L1TCalorimeterPhase2/interface/DigitizedTowerCorrelator.h"
0023 #include "DataFormats/L1TCalorimeterPhase2/interface/DigitizedClusterGT.h"
0024 #include "DataFormats/L1TCalorimeterPhase2/interface/GCTEmDigiCluster.h"
0025 #include "DataFormats/L1TCalorimeterPhase2/interface/GCTHadDigiCluster.h"
0026 
0027 #include "DataFormats/L1Trigger/interface/BXVector.h"
0028 #include "DataFormats/L1Trigger/interface/EGamma.h"
0029 
0030 namespace p2eg {
0031 
0032   static constexpr int n_towers_Eta = 34;
0033   static constexpr int n_towers_Phi = 72;
0034   static constexpr int n_towers_halfPhi = 36;
0035   static constexpr int n_towers_cardEta = 17;  // new: equivalent to n_towers_per_link
0036   static constexpr int n_towers_cardPhi = 4;
0037   static constexpr int n_crystals_cardEta = (n_towers_Eta * n_towers_cardEta);
0038   static constexpr int n_crystals_cardPhi = (n_towers_Phi * n_towers_cardPhi);
0039 
0040   // outputs
0041   static constexpr int n_links_card = 4;      // 4 links per card
0042   static constexpr int n_clusters_link = 2;   // 2 clusters sent out in each link
0043   static constexpr int n_clusters_4link = 8;  // 8 clusters sent out in 4 links
0044   static constexpr int n_towers_per_link = 17;
0045 
0046   static constexpr int CRYSTALS_IN_TOWER_ETA = 5;
0047   static constexpr int CRYSTALS_IN_TOWER_PHI = 5;
0048 
0049   static constexpr int TOWER_IN_ETA = 3;  // number of towers in eta, in one 3x4 region (barrel)
0050   static constexpr int TOWER_IN_PHI = 4;  // number of towers in phi, in one 3x4 region (barrel)
0051 
0052   static constexpr int CRYSTAL_IN_ETA = 15;  // number of crystals in eta, in one 3x4 region (barrel)
0053   static constexpr int CRYSTAL_IN_PHI = 20;  // number of crystals in phi, in one 3x4 region (barrel)
0054 
0055   static constexpr float ECAL_eta_range = 1.4841;
0056   static constexpr float half_crystal_size = 0.00873;
0057 
0058   static constexpr float slideIsoPtThreshold = 80;
0059   static constexpr float a0_80 = 0.85, a1_80 = 0.0080, a0 = 0.21;                        // passes_iso
0060   static constexpr float b0 = 0.38, b1 = 1.9, b2 = 0.05;                                 // passes_looseTkiso
0061   static constexpr float c0_ss = 0.94, c1_ss = 0.052, c2_ss = 0.044;                     // passes_ss
0062   static constexpr float d0 = 0.96, d1 = 0.0003;                                         // passes_photon
0063   static constexpr float e0_looseTkss = 0.944, e1_looseTkss = 0.65, e2_looseTkss = 0.4;  // passes_looseTkss
0064   static constexpr float loose_ss_offset = 0.89;  // loosen ss requirement in two eta bins
0065   static constexpr float eta0_loose_ss = 0.4, eta1_loose_ss = 0.5, eta2_loose_ss = 0.75, eta3_loose_ss = 0.85;
0066   static constexpr float high_pt_threshold = 130.;  // apply iso and ss requirements below this threshold
0067   static constexpr float cut_500_MeV = 0.5;
0068 
0069   static constexpr float ECAL_LSB = 0.5;  // to convert from int to float (GeV) multiply by LSB
0070   static constexpr float HCAL_LSB = 0.5;
0071 
0072   static constexpr int N_CLUSTERS_PER_REGION = 4;  // number of clusters per ECAL region
0073   static constexpr int N_REGIONS_PER_CARD = 6;     // number of ECAL regions per card
0074 
0075   // GCT constants
0076   static constexpr int N_RCTCARDS_PHI = 8;
0077   static constexpr int N_RCTGCT_FIBERS = 4;
0078   static constexpr int N_RCTTOWERS_FIBER = 17;
0079   static constexpr int N_RCTCLUSTERS_FIBER = 2;
0080 
0081   static constexpr int N_GCTCARDS = 3;
0082   static constexpr int N_GCTCORR_FIBERS = 48;
0083   static constexpr int N_GCTTOWERS_FIBER = 17;
0084   static constexpr int N_GCTCLUSTERS_FIBER = 2;
0085 
0086   static constexpr int N_GCTINTERNAL_FIBERS = 64;
0087   static constexpr int N_GCTPOSITIVE_FIBERS = 32;
0088   static constexpr int N_GCTETA = 34;
0089   static constexpr int N_GCTPHI = 32;
0090 
0091   // for emulator: "top" of the GCT card in phi is tower idx 20, for GCT card #0:
0092   static constexpr int GCTCARD_0_TOWER_IPHI_OFFSET = 20;
0093   // same but for GCT cards #1 and 2 (cards wrap around phi = 180 degrees):
0094   static constexpr int GCTCARD_1_TOWER_IPHI_OFFSET = 44;
0095   static constexpr int GCTCARD_2_TOWER_IPHI_OFFSET = 68;
0096 
0097   static constexpr int N_GCTTOWERS_CLUSTER_ISO_ONESIDE = 5;  // window size of isolation sum (5x5 in towers)
0098 
0099   // Outputs to correlator
0100   // Fixed number of EG and PF clusters per RCT card in each output SLR
0101   static const int N_EG_CLUSTERS_PER_RCT_CARD = 9;
0102   static const int N_PF_CLUSTERS_PER_RCT_CARD = 12;
0103   // Height of one SLR region in phi in degrees
0104   static constexpr float PHI_RANGE_PER_SLR_DEGREES = 120;
0105 
0106   //////////////////////////////////////////////////////////////////////////
0107   // RCT: indexing helper functions
0108   //////////////////////////////////////////////////////////////////////////
0109 
0110   // Assert that the card index is within bounds. (Valid cc: 0 to 35, since there are 36 RCT cards)
0111   inline bool isValidCard(int cc) { return ((cc > -1) && (cc < 36)); }
0112 
0113   // RCT Cards: need to know their min/max crystal boundaries.
0114 
0115   // For a card (ranging from 0 to 35, since there are 36 cards), return the iEta of the crystal with max iEta.
0116   // This represents the card boundaries in eta (identical to getEtaMax_card in the original emulator)
0117   inline int getCard_iEtaMax(int cc) {
0118     assert(isValidCard(cc));
0119 
0120     int etamax = 0;
0121     if (cc % 2 == 0)                                            // Even card: negative eta
0122       etamax = (n_towers_cardEta * CRYSTALS_IN_TOWER_ETA - 1);  // First eta half. 5 crystals in eta in 1 tower.
0123     else                                                        // Odd card: positive eta
0124       etamax = (n_towers_Eta * CRYSTALS_IN_TOWER_ETA - 1);
0125     return etamax;
0126   }
0127 
0128   // Same as above but for minimum iEta.
0129   inline int getCard_iEtaMin(int cc) {
0130     int etamin = 0;
0131     if (cc % 2 == 0)  // Even card: negative eta
0132       etamin = (0);
0133     else  // Odd card: positive eta
0134       etamin = (n_towers_cardEta * CRYSTALS_IN_TOWER_ETA);
0135     return etamin;
0136   }
0137 
0138   // Same as above but for maximum iPhi.
0139   inline int getCard_iPhiMax(int cc) {
0140     int phimax = ((cc / 2) + 1) * 4 * CRYSTALS_IN_TOWER_PHI - 1;
0141     return phimax;
0142   }
0143 
0144   // Same as above but for minimum iPhi.
0145   inline int getCard_iPhiMin(int cc) {
0146     int phimin = (cc / 2) * 4 * CRYSTALS_IN_TOWER_PHI;
0147     return phimin;
0148   }
0149 
0150   // Given the RCT card number (0-35), get the crystal iEta of the "bottom left" corner
0151   inline int getCard_refCrystal_iEta(int cc) {
0152     if ((cc % 2) == 1) {  // if cc is odd (positive eta)
0153       return (17 * CRYSTALS_IN_TOWER_ETA);
0154     } else {  // if cc is even (negative eta) the bottom left corner is further in eta, hence +4
0155       return ((16 * CRYSTALS_IN_TOWER_ETA) + 4);
0156     }
0157   }
0158 
0159   // Given the RCT card number (0-35), get the global crystal iPhi of the "bottom left" corner (0- 71*5)
0160   inline int getCard_refCrystal_iPhi(int cc) {
0161     if ((cc % 2) == 1) {
0162       // if cc is odd: positive eta
0163       return int(cc / 2) * TOWER_IN_PHI * CRYSTALS_IN_TOWER_PHI;
0164     } else {
0165       // if cc is even, the bottom left corner is further in phi, hence the +4 and -1
0166       return (((int(cc / 2) * TOWER_IN_PHI) + 4) * CRYSTALS_IN_TOWER_PHI) - 1;
0167     }
0168   }
0169 
0170   // Towers: Go from real (eta, phi) to tower absolute ID
0171 
0172   /* 
0173   * For a real eta, get the tower absolute Eta index (possible values are 0-33, since there
0174   * are 34 towers in eta. (Adapted from getTower_absoluteEtaID)
0175   */
0176   inline int getTower_absEtaID(float eta) {
0177     float size_cell = 2 * ECAL_eta_range / n_towers_Eta;
0178     int etaID = int((eta + ECAL_eta_range) / size_cell);
0179     return etaID;
0180   }
0181 
0182   /* 
0183   * Same as above, but for phi.
0184   * Possible values range from 0-71 (Adapted from getTower_absolutePhiID)
0185   */
0186   inline int getTower_absPhiID(float phi) {
0187     float size_cell = 2 * M_PI / n_towers_Phi;
0188     int phiID = int((phi + M_PI) / size_cell);
0189     return phiID;
0190   }
0191 
0192   // Towers: Go from firmware specifics (RCT card, tower number in link, and link number (all firmware convention))
0193   //         to tower absolute ID.
0194 
0195   /* 
0196   * Get the global tower iEta (0-31) from the firmware card, tower number (0-16), and link (0-3). Respects the fact that
0197   * in the firmware, negative eta cards are "rotated" (link 0, tower 0) starts in the "top right" corner if we
0198   * look at a diagram of the barrel region.
0199   */
0200   inline int getAbsID_iEta_fromFirmwareCardTowerLink(int nCard, int nTower, int nLink) {
0201     // iEta only depends on the tower position in the link
0202     (void)nCard;
0203     (void)nLink;
0204     if ((nCard % 2) == 1) {  // if cc is odd (positive eta), e.g. nTower = 0 will correspond to absolute iEta ID 17.
0205       return n_towers_per_link + nTower;
0206     } else {  // if cc is even (negative eta): e.g. nTower = 0 will correspond to absolute iEta ID 16.
0207       return (16 - nTower);
0208     }
0209   }
0210 
0211   /*
0212   * Get the global tower iPhi (0-71) from the firmware card, tower number (0-16), and link (0-3).
0213   */
0214   inline int getAbsID_iPhi_fromFirmwareCardTowerLink(int nCard, int nTower, int nLink) {
0215     // iPhi only depends on the card and link number
0216     (void)nTower;
0217     if ((nCard % 2) == 1) {  // if cc is odd (positive eta),
0218       // e.g. cc=3, link #2, global iPhi = int(3/2) * 4 + 2 = 6
0219       return (int(nCard / 2) * TOWER_IN_PHI) + nLink;
0220     } else {  // if cc is even (negative eta)
0221       // e.g. cc=4, link #2, global iPhi = int(4/2) * 4 + (4 - 2 - 1)
0222       //                                 = 2*4 + 1
0223       //                                 = 9
0224       // minus one is because TOWER_IN_PHI is 4
0225       return (int(nCard / 2) * TOWER_IN_PHI) + (TOWER_IN_PHI - nLink - 1);
0226     }
0227   }
0228 
0229   // Towers: Go from absolute ID, back to real eta and phi.
0230 
0231   /*
0232   * From the tower absolute ID in eta (0-33), get the real eta of the tower center
0233   * Same as getTowerEta_fromAbsoluteID in previous CMSSW emulator
0234   */
0235   inline float getTowerEta_fromAbsID(int id) {
0236     float size_cell = 2 * ECAL_eta_range / n_towers_Eta;
0237     float eta = (id * size_cell) - ECAL_eta_range + 0.5 * size_cell;
0238     return eta;
0239   }
0240 
0241   /*
0242   * From the tower absolute ID in phi (0-71), get the real phi of the tower center
0243   * Same as getTowerPhi_fromAbsoluteID in previous CMSSW emulator
0244   */
0245   inline float getTowerPhi_fromAbsID(int id) {
0246     float size_cell = 2 * M_PI / n_towers_Phi;
0247     float phi = (id * size_cell) - M_PI + 0.5 * size_cell;
0248     return phi;
0249   }
0250 
0251   /* 
0252   * Get the RCT card region that a crystal is in, given the "local" iEta of the crystal 
0253   * 0 is region closest to eta = 0. Regions 0, 1, 2, 3, 4 are in the barrel, Region 5 is in overlap
0254   */
0255   inline int getRegionNumber(const int local_iEta) {
0256     int no = int(local_iEta / (TOWER_IN_ETA * CRYSTALS_IN_TOWER_ETA));
0257     assert(no < 6);
0258     return no;
0259   }
0260 
0261   /*******************************************************************/
0262   /* RCT classes and structs                                         */
0263   /*******************************************************************/
0264 
0265   /* 
0266    * Represents one input HCAL or ECAL hit.
0267    */
0268   class SimpleCaloHit {
0269   private:
0270     float pt_ = 0;
0271     float energy_ = 0.;
0272     ap_uint<10> et_uint_;
0273     GlobalVector position_;  // As opposed to GlobalPoint, so we can add them (for weighted average)
0274     HcalDetId id_hcal_;
0275     EBDetId id_;
0276 
0277   public:
0278     // tool functions
0279     inline void setPt() { pt_ = (position_.mag2() > 0) ? energy_ * sin(position_.theta()) : 0; };
0280     inline void setEnergy(float et) { energy_ = et / sin(position_.theta()); };
0281     inline void setEt_uint(ap_uint<10> et_uint) { et_uint_ = et_uint; }
0282     inline void setPosition(const GlobalVector& pos) { position_ = pos; };
0283     inline void setIdHcal(const HcalDetId& idhcal) { id_hcal_ = idhcal; };
0284     inline void setId(const EBDetId& id) { id_ = id; };
0285 
0286     inline float pt() const { return pt_; };
0287     inline float energy() const { return energy_; };
0288     inline ap_uint<10> et_uint() const { return et_uint_; };
0289     inline const GlobalVector& position() const { return position_; };
0290     inline const EBDetId& id() const { return id_; };
0291 
0292     /* 
0293        * Get crystal's iEta from real eta. (identical to getCrystal_etaID in L1EGammaCrystalsEmulatorProducer.cc)
0294        * This "global" iEta ranges from 0 to (33*5) since there are 34 towers in eta in the full detector, 
0295        * each with five crystals in eta.
0296        */
0297     int crystaliEta(void) const {
0298       float size_cell = 2 * ECAL_eta_range / (CRYSTALS_IN_TOWER_ETA * n_towers_Eta);
0299       int iEta = int((position().eta() + ECAL_eta_range) / size_cell);
0300       return iEta;
0301     }
0302 
0303     /* 
0304        * Get crystal's iPhi from real phi. (identical to getCrystal_phiID in L1EGammaCrystalsEmulatorProducer.cc)
0305        * This "global" iPhi ranges from 0 to (71*5) since there are 72 towers in phi in the full detector, each with five crystals in eta.
0306        */
0307     int crystaliPhi(void) const {
0308       float phi = position().phi();
0309       float size_cell = 2 * M_PI / (CRYSTALS_IN_TOWER_PHI * n_towers_Phi);
0310       int iPhi = int((phi + M_PI) / size_cell);
0311       return iPhi;
0312     }
0313 
0314     /*
0315        * Check if it falls within the boundary of a card.
0316        */
0317     bool isInCard(int cc) const {
0318       return (crystaliPhi() <= getCard_iPhiMax(cc) && crystaliPhi() >= getCard_iPhiMin(cc) &&
0319               crystaliEta() <= getCard_iEtaMax(cc) && crystaliEta() >= getCard_iEtaMin(cc));
0320     };
0321 
0322     /*
0323       * For a crystal with real eta, and falling in card cc, get its local iEta 
0324       * relative to the bottom left corner of the card (possible local iEta ranges from 0 to 17 * 5,
0325       * since in one card, there are 17 towers in eta, each with 5 crystals in eta.
0326       */
0327     int crystalLocaliEta(int cc) const { return abs(getCard_refCrystal_iEta(cc) - crystaliEta()); }
0328 
0329     /*
0330       * Same as above, but for iPhi (possible local iPhi ranges from 0 to (3*5), since in one card,
0331       * there are 4 towers in phi, each with 5 crystals in phi.
0332       */
0333     int crystalLocaliPhi(int cc) const { return abs(getCard_refCrystal_iPhi(cc) - crystaliPhi()); }
0334 
0335     /*
0336        * Print hit info
0337        */
0338     void printHitInfo(std::string description = "") const {
0339       std::cout << "[printHitInfo]: [" << description << "]"
0340                 << " hit with energy " << pt() << " at eta " << position().eta() << ", phi " << position().phi()
0341                 << std::endl;
0342     }
0343   };
0344 
0345   /*******************************************************************/
0346 
0347   /*
0348   * linkECAL class: represents one ECAL link (one tower: 5x5 crystals)
0349   */
0350 
0351   class linkECAL {
0352   private:
0353     ap_uint<10> crystalE[CRYSTALS_IN_TOWER_ETA][CRYSTALS_IN_TOWER_PHI];
0354 
0355   public:
0356     // constructor
0357     linkECAL() {}
0358 
0359     // Set members
0360     inline void zeroOut() {  // zero out the crystalE array
0361       for (int i = 0; i < CRYSTALS_IN_TOWER_ETA; i++) {
0362         for (int j = 0; j < CRYSTALS_IN_TOWER_PHI; j++) {
0363           crystalE[i][j] = 0;
0364         }
0365       }
0366     };
0367     inline void setCrystalE(int iEta, int iPhi, ap_uint<10> energy) {
0368       assert(iEta < CRYSTALS_IN_TOWER_ETA);
0369       assert(iPhi < CRYSTALS_IN_TOWER_PHI);
0370       crystalE[iEta][iPhi] = energy;
0371     };
0372     inline void addCrystalE(int iEta, int iPhi, ap_uint<10> energy) {
0373       assert(iEta < CRYSTALS_IN_TOWER_ETA);
0374       assert(iPhi < CRYSTALS_IN_TOWER_PHI);
0375       crystalE[iEta][iPhi] += energy;
0376     };
0377 
0378     // Access members
0379     inline ap_uint<10> getCrystalE(int iEta, int iPhi) {
0380       assert(iEta < 5);
0381       assert(iPhi < 5);
0382       return crystalE[iEta][iPhi];
0383     };
0384   };
0385 
0386   /*******************************************************************/
0387 
0388   /*
0389   * region3x4 class: represents one 3x4 ECAL region. The region stores no
0390   *                  information about which card it is located in.
0391   *                  idx: 0-4. Region 0 is the one closest to eta = 0, counting outwards in eta  
0392   */
0393 
0394   class region3x4 {
0395   private:
0396     int idx_ = -1;
0397     std::array<std::array<linkECAL, TOWER_IN_PHI>, TOWER_IN_ETA> linksECAL;  // 3x4 in towers
0398 
0399   public:
0400     // constructor
0401     region3x4() { idx_ = -1; }
0402 
0403     // set members
0404     inline void zeroOut() {
0405       for (int i = 0; i < TOWER_IN_ETA; i++) {
0406         for (int j = 0; j < TOWER_IN_PHI; j++) {
0407           linksECAL[i][j].zeroOut();
0408         }
0409       }
0410     };
0411     inline void setIdx(int idx) { idx_ = idx; };
0412 
0413     // get members
0414     inline float getIdx() const { return idx_; };
0415     inline linkECAL& getLinkECAL(int iEta, int iPhi) { return linksECAL[iEta][iPhi]; };
0416   };
0417 
0418   /*******************************************************************/
0419 
0420   /*
0421   * towerHCAL class: represents one HCAL tower
0422   */
0423 
0424   class towerHCAL {
0425   private:
0426     ap_uint<10> et;
0427     ap_uint<6> fb;
0428 
0429   public:
0430     // constructor
0431     towerHCAL() {
0432       et = 0;
0433       fb = 0;
0434     };
0435 
0436     // set members
0437     inline void zeroOut() {
0438       et = 0;
0439       fb = 0;
0440     };
0441     inline void addEt(ap_uint<10> newEt) { et += newEt; };
0442 
0443     // get members
0444     inline ap_uint<10> getEt() { return et; };
0445   };
0446 
0447   /*******************************************************************/
0448 
0449   /*
0450   * towers3x4 class: represents 3x4 array of HCAL towers. idx = 0, 1, ... 4 are the barrel gion
0451   */
0452 
0453   class towers3x4 {
0454   private:
0455     int idx_ = -1;
0456     std::array<std::array<towerHCAL, TOWER_IN_PHI>, TOWER_IN_ETA> towersHCAL;  // 3x4 in towers
0457 
0458   public:
0459     // constructor
0460     towers3x4() { idx_ = -1; };
0461 
0462     // set members
0463     inline void zeroOut() {
0464       for (int i = 0; i < TOWER_IN_ETA; i++) {
0465         for (int j = 0; j < TOWER_IN_PHI; j++) {
0466           towersHCAL[i][j].zeroOut();
0467         }
0468       }
0469     };
0470     inline void setIdx(int idx) { idx_ = idx; };
0471 
0472     // get members
0473     inline float getIdx() const { return idx_; };
0474     inline towerHCAL& getTowerHCAL(int iEta, int iPhi) { return towersHCAL[iEta][iPhi]; };
0475   };
0476 
0477   /*******************************************************************/
0478 
0479   /* 
0480    * card class: represents one RCT card. Each card has five 3x4 regions and one 2x4 region,
0481    *             which is represented by a 3x4 region with its third row zero'd out.
0482    *             idx 0-35: odd values of cardIdx span eta = 0 to eta = 1.41 
0483    *                       even values of cardIdx span eta = -1.41 to eta = 0
0484    *             The realEta and realPhi arrays store the (eta, phi) of the center of the towers.
0485    */
0486 
0487   class card {
0488   private:
0489     int idx_ = -1;
0490     std::array<region3x4, N_REGIONS_PER_CARD> card3x4Regions;
0491     std::array<towers3x4, N_REGIONS_PER_CARD> card3x4Towers;
0492 
0493   public:
0494     // constructor
0495     card() {
0496       idx_ = -1;
0497       for (int i = 0; i < N_REGIONS_PER_CARD; i++) {
0498         card3x4Regions[i].setIdx(i);
0499         card3x4Regions[i].zeroOut();
0500         card3x4Towers[i].setIdx(i);
0501         card3x4Towers[i].zeroOut();
0502       }
0503     };
0504 
0505     // set members
0506     inline void setIdx(int idx) { idx_ = idx; };
0507     inline void zeroOut() {
0508       for (int i = 0; i < N_REGIONS_PER_CARD; i++) {
0509         card3x4Regions[i].zeroOut();
0510         card3x4Towers[i].zeroOut();
0511       };
0512     };
0513 
0514     // get members
0515     inline float getIdx() const { return idx_; };
0516     inline region3x4& getRegion3x4(int idx) {
0517       assert(idx < N_REGIONS_PER_CARD);
0518       return card3x4Regions[idx];
0519     }
0520     inline towers3x4& getTowers3x4(int idx) {
0521       assert(idx < N_REGIONS_PER_CARD);
0522       return card3x4Towers[idx];
0523     }
0524   };
0525 
0526   /*******************************************************************/
0527 
0528   /*
0529    *  Crystal class for RCT
0530    */
0531 
0532   class crystal {
0533   public:
0534     ap_uint<10> energy;
0535 
0536     crystal() {
0537       energy = 0;
0538       //    timing = 0;
0539     }
0540 
0541     crystal(ap_uint<10> energy) {  // To-do: add timing information
0542       this->energy = energy;
0543       //    this->timing = 0;
0544     }
0545 
0546     crystal& operator=(const crystal& rhs) {
0547       energy = rhs.energy;
0548       //    timing = rhs.timing;
0549       return *this;
0550     }
0551   };
0552 
0553   /*
0554   * crystalMax class for RCT
0555   */
0556   class crystalMax {
0557   public:
0558     ap_uint<10> energy;
0559     uint8_t phiMax;
0560     uint8_t etaMax;
0561 
0562     crystalMax() {
0563       energy = 0;
0564       phiMax = 0;
0565       etaMax = 0;
0566     }
0567   };
0568 
0569   class ecaltp_t {
0570   public:
0571     ap_uint<10> energy;
0572     ap_uint<5> phi;
0573     ap_uint<5> eta;
0574   };
0575 
0576   class etaStrip_t {
0577   public:
0578     ecaltp_t cr0;
0579     ecaltp_t cr1;
0580     ecaltp_t cr2;
0581     ecaltp_t cr3;
0582     ecaltp_t cr4;
0583     ecaltp_t cr5;
0584     ecaltp_t cr6;
0585     ecaltp_t cr7;
0586     ecaltp_t cr8;
0587     ecaltp_t cr9;
0588     ecaltp_t cr10;
0589     ecaltp_t cr11;
0590     ecaltp_t cr12;
0591     ecaltp_t cr13;
0592     ecaltp_t cr14;
0593     ecaltp_t cr15;
0594     ecaltp_t cr16;
0595     ecaltp_t cr17;
0596     ecaltp_t cr18;
0597     ecaltp_t cr19;
0598   };
0599 
0600   class ecalRegion_t {
0601   public:
0602     etaStrip_t etaStrip0;
0603     etaStrip_t etaStrip1;
0604     etaStrip_t etaStrip2;
0605     etaStrip_t etaStrip3;
0606     etaStrip_t etaStrip4;
0607     etaStrip_t etaStrip5;
0608     etaStrip_t etaStrip6;
0609     etaStrip_t etaStrip7;
0610     etaStrip_t etaStrip8;
0611     etaStrip_t etaStrip9;
0612     etaStrip_t etaStrip10;
0613     etaStrip_t etaStrip11;
0614     etaStrip_t etaStrip12;
0615     etaStrip_t etaStrip13;
0616     etaStrip_t etaStrip14;
0617   };
0618 
0619   class etaStripPeak_t {
0620   public:
0621     ecaltp_t pk0;
0622     ecaltp_t pk1;
0623     ecaltp_t pk2;
0624     ecaltp_t pk3;
0625     ecaltp_t pk4;
0626     ecaltp_t pk5;
0627     ecaltp_t pk6;
0628     ecaltp_t pk7;
0629     ecaltp_t pk8;
0630     ecaltp_t pk9;
0631     ecaltp_t pk10;
0632     ecaltp_t pk11;
0633     ecaltp_t pk12;
0634     ecaltp_t pk13;
0635     ecaltp_t pk14;
0636   };
0637 
0638   class tower_t {
0639   public:
0640     ap_uint<16> data;
0641 
0642     tower_t() { data = 0; }
0643 
0644     tower_t(ap_uint<12> et, ap_uint<4> hoe) { data = (et) | (((ap_uint<16>)hoe) << 12); }
0645 
0646     ap_uint<12> et() { return (data & 0xFFF); }
0647     ap_uint<4> hoe() { return ((data >> 12) & 0xF); }
0648 
0649     float getEt() { return (float)et() * ECAL_LSB; }
0650 
0651     operator uint16_t() { return (uint16_t)data; }
0652 
0653     // Only for ECAL towers! Apply calibration and modify the et() value.
0654     void applyCalibration(float factor) {
0655       // Get the new pT as a float
0656       float newEt = getEt() * factor;
0657 
0658       // Convert the new pT to an unsigned int (16 bits so we can take the logical OR with the bit mask later)
0659       ap_uint<16> newEt_uint = (ap_uint<16>)(int)(newEt / ECAL_LSB);
0660       // Make sure the first four bits are zero
0661       newEt_uint = (newEt_uint & 0x0FFF);
0662 
0663       // Modify 'data'
0664       ap_uint<16> bitMask = 0xF000;  // last twelve digits are zero
0665       data = (data & bitMask);       // zero out the last twelve digits
0666       data = (data | newEt_uint);    // write in the new ET
0667     }
0668     /*
0669      * For towers: Calculate H/E ratio given the ECAL and HCAL energies and modify the hoe() value.
0670      */
0671     void addHoverEToTower(ap_uint<12> ECAL, ap_uint<12> HCAL) {
0672       ap_uint<4> hoeOut;
0673       ap_uint<1> hoeLSB = 0;
0674       ap_uint<4> hoe = 0;
0675       ap_uint<12> A;
0676       ap_uint<12> B;
0677 
0678       A = (ECAL > HCAL) ? ECAL : HCAL;
0679       B = (ECAL > HCAL) ? HCAL : ECAL;
0680 
0681       if (ECAL == 0 || HCAL == 0 || HCAL >= ECAL)
0682         hoeLSB = 0;
0683       else
0684         hoeLSB = 1;
0685       if (A > B) {
0686         if (A > 2 * B)
0687           hoe = 0x1;
0688         if (A > 4 * B)
0689           hoe = 0x2;
0690         if (A > 8 * B)
0691           hoe = 0x3;
0692         if (A > 16 * B)
0693           hoe = 0x4;
0694         if (A > 32 * B)
0695           hoe = 0x5;
0696         if (A > 64 * B)
0697           hoe = 0x6;
0698         if (A > 128 * B)
0699           hoe = 0x7;
0700       }
0701       hoeOut = hoeLSB | (hoe << 1);
0702       ap_uint<16> hoeOutLong =
0703           ((((ap_uint<16>)hoeOut) << 12) | 0x0000);  // e.g. 0b ____ 0000 0000 0000 where ___ are the hoe digits
0704       // Take the logical OR to preserve the saturation and tower ET bits
0705       ap_uint<16> bitMask = 0x0FFF;  // 0000 1111 1111 1111 : zero-out the HoE bits
0706       data = (data & bitMask);       // zero-out the HoE bits
0707       data = (data | hoeOutLong);    // write in the new HoE bits
0708     }
0709   };
0710 
0711   class clusterInfo {
0712   public:
0713     ap_uint<10> seedEnergy;
0714     ap_uint<15> energy;
0715     ap_uint<15> et5x5;
0716     ap_uint<15> et2x5;
0717     ap_uint<5> phiMax;
0718     ap_uint<5> etaMax;
0719     ap_uint<2> brems;
0720 
0721     clusterInfo() {
0722       seedEnergy = 0;
0723       energy = 0;
0724       et5x5 = 0;
0725       et2x5 = 0;
0726       phiMax = 0;
0727       etaMax = 0;
0728       brems = 0;
0729     }
0730   };
0731 
0732   //--------------------------------------------------------//
0733 
0734   /*
0735   * Cluster class for RCT.
0736   */
0737   class Cluster {
0738   public:
0739     ap_uint<28> data;
0740     int regionIdx;       // region is 0 through 4 in barrel, 5 for overlap in barrel+endcap
0741     float calib;         //  ECAL energy calibration factor
0742     ap_uint<2> brems;    // brems flag
0743     ap_uint<15> et5x5;   // et5x5 sum in towers around cluster
0744     ap_uint<15> et2x5;   // et2x5 sum in towers around cluster
0745     bool is_ss;          // shower shape flag
0746     bool is_looseTkss;   // loose Tk shower shape flag
0747     bool is_iso;         // isolation flag (not computed until GCT)
0748     bool is_looseTkiso;  // isolation flag (not computed until GCT)
0749 
0750     Cluster() {
0751       data = 0;
0752       regionIdx = -1;
0753       calib = 1.0;
0754       brems = 0;
0755       et5x5 = 0;
0756       et2x5 = 0;
0757       is_ss = false;
0758       is_looseTkss = false;
0759       is_iso = false;
0760       is_looseTkiso = false;
0761     }
0762 
0763     Cluster(ap_uint<12> clusterEnergy,
0764             ap_uint<5> towerEta,
0765             ap_uint<2> towerPhi,
0766             ap_uint<3> clusterEta,
0767             ap_uint<3> clusterPhi,
0768             ap_uint<3> satur,
0769             ap_uint<15> clusterEt5x5 = 0,
0770             ap_uint<15> clusterEt2x5 = 0,
0771             ap_uint<2> clusterBrems = 0,
0772             float clusterCalib = 1.0,
0773             bool cluster_is_ss = false,
0774             bool cluster_is_looseTkss = false,
0775             bool cluster_is_iso = false,
0776             bool cluster_is_looseTkiso = false,
0777             int clusterRegionIdx = 0) {
0778       data = (clusterEnergy) | (((ap_uint<32>)towerEta) << 12) | (((ap_uint<32>)towerPhi) << 17) |
0779              (((ap_uint<32>)clusterEta) << 19) | (((ap_uint<32>)clusterPhi) << 22) | (((ap_uint<32>)satur) << 25);
0780       regionIdx = clusterRegionIdx;
0781       et5x5 = clusterEt5x5;
0782       et2x5 = clusterEt2x5;
0783       brems = clusterBrems;
0784       calib = clusterCalib;
0785       is_ss = cluster_is_ss;
0786       is_looseTkss = cluster_is_looseTkss;
0787       is_iso = cluster_is_iso;
0788       is_looseTkiso = cluster_is_looseTkiso;
0789     }
0790 
0791     void setRegionIdx(int regIdx) { regionIdx = regIdx; }  // Newly added
0792 
0793     ap_uint<12> clusterEnergy() const { return (data & 0xFFF); }
0794     ap_uint<5> towerEta() const { return ((data >> 12) & 0x1F); }  // goes from 0 to 3 (need region for full info)
0795     ap_uint<2> towerPhi() const { return ((data >> 17) & 0x3); }
0796     ap_uint<3> clusterEta() const { return ((data >> 19) & 0x7); }
0797     ap_uint<3> clusterPhi() const { return ((data >> 22) & 0x7); }
0798     ap_uint<3> satur() const { return ((data >> 25) & 0x7); }
0799     ap_uint<15> uint_et2x5() const { return et2x5; }
0800     ap_uint<15> uint_et5x5() const { return et5x5; }
0801 
0802     operator uint32_t() const { return (ap_uint<28>)data; }
0803     int region() const { return regionIdx; }  // Newly added
0804     int getBrems() const { return (int)brems; }
0805     float getCalib() const { return (float)calib; }
0806     float getPt() const { return ((float)clusterEnergy() * ECAL_LSB); }  // Return pT as a float
0807     float getEt5x5() const { return ((float)et5x5 * ECAL_LSB); }         // Return ET5x5 as a float
0808     float getEt2x5() const { return ((float)et2x5 * ECAL_LSB); }         // Return ET2x5 as a float
0809 
0810     int towerEtaInCard() { return ((int)(region() * TOWER_IN_ETA) + towerEta()); }
0811 
0812     bool getIsSS() { return is_ss; }
0813     bool getIsLooseTkss() { return is_looseTkss; }
0814     bool getIsIso() { return is_iso; }
0815     bool getIsLooseTkIso() { return is_looseTkiso; }
0816 
0817     /*
0818       * Apply calibration (float) to the pT in-place.
0819       */
0820     void applyCalibration(float factor) {
0821       float newPt = getPt() * factor;
0822       // Convert the new pT to an unsigned int, 28 bits to take the logical OR with the mask
0823       ap_uint<28> newPt_uint = (ap_uint<28>)(int)(newPt / ECAL_LSB);
0824       // Make sure that the new pT only takes up the last twelve bits
0825       newPt_uint = (newPt_uint & 0x0000FFF);
0826 
0827       // Modify 'data'
0828       ap_uint<28> bitMask = 0xFFFF000;  // last twelve digits are zero
0829       data = (data & bitMask);          // zero out the last twelve digits
0830       data = (data | newPt_uint);       // write in the new ET
0831     }
0832 
0833     // Get Cluster crystal iEta from card number, region, tower eta, and cluster eta indices
0834     const int crystaliEtaFromCardRegionInfo(int cc) {
0835       int crystalEta_in_card =
0836           ((region() * TOWER_IN_ETA * CRYSTALS_IN_TOWER_ETA) + (towerEta() * CRYSTALS_IN_TOWER_ETA) + clusterEta());
0837       if ((cc % 2) == 1) {
0838         return (getCard_refCrystal_iEta(cc) + crystalEta_in_card);
0839       } else {
0840         // if card is even (negative eta)
0841         return (getCard_refCrystal_iEta(cc) - crystalEta_in_card);
0842       }
0843     }
0844 
0845     // Get Cluster crystal iPhi from card number, region, tower eta, and cluster phi indices
0846     const int crystaliPhiFromCardRegionInfo(int cc) {
0847       int crystalPhi_in_card = (towerPhi() * CRYSTALS_IN_TOWER_PHI) + clusterPhi();
0848       if ((cc % 2) == 1) {
0849         // if card is odd (positive eta)
0850         return (getCard_refCrystal_iPhi(cc) + crystalPhi_in_card);
0851       } else {
0852         // if card is even (negative eta)
0853         return (getCard_refCrystal_iPhi(cc) - crystalPhi_in_card);
0854       }
0855     }
0856 
0857     // Get real eta
0858     const float realEta(int cc) {
0859       float size_cell = 2 * ECAL_eta_range / (CRYSTALS_IN_TOWER_ETA * n_towers_Eta);
0860       return crystaliEtaFromCardRegionInfo(cc) * size_cell - ECAL_eta_range + half_crystal_size;
0861     }
0862 
0863     // Get real phi
0864     const float realPhi(int cc) {
0865       float size_cell = 2 * M_PI / (CRYSTALS_IN_TOWER_PHI * n_towers_Phi);
0866       return crystaliPhiFromCardRegionInfo(cc) * size_cell - M_PI + half_crystal_size;
0867     }
0868 
0869     // Print info
0870     void printClusterInfo(int cc, std::string description = "") {
0871       std::cout << "[Print Cluster class info:] [" << description << "]: "
0872                 << "card " << cc << ", "
0873                 << "et (float): " << getPt() << ", "
0874                 << "eta: " << realEta(cc) << ", "
0875                 << "phi: " << realPhi(cc) << std::endl;
0876     }
0877   };
0878 
0879   /*
0880   * Compare the ET of two clusters (pass this to std::sort to get clusters sorted in decreasing ET).
0881   */
0882   inline bool compareClusterET(const Cluster& lhs, const Cluster& rhs) {
0883     return (lhs.clusterEnergy() > rhs.clusterEnergy());
0884   }
0885 
0886   /*******************************************************************/
0887   /* RCT helper functions                                            */
0888   /*******************************************************************/
0889   ecalRegion_t initStructure(crystal temporary[CRYSTAL_IN_ETA][CRYSTAL_IN_PHI]);
0890   ecaltp_t bestOf2(const ecaltp_t ecaltp0, const ecaltp_t ecaltp1);
0891   ecaltp_t getPeakBin20N(const etaStrip_t etaStrip);
0892   crystalMax getPeakBin15N(const etaStripPeak_t etaStrip);
0893   void getECALTowersEt(crystal tempX[CRYSTAL_IN_ETA][CRYSTAL_IN_PHI], ap_uint<12> towerEt[12]);
0894   clusterInfo getClusterPosition(const ecalRegion_t ecalRegion);
0895   Cluster packCluster(ap_uint<15>& clusterEt, ap_uint<5>& etaMax_t, ap_uint<5>& phiMax_t);
0896   void removeClusterFromCrystal(crystal temp[CRYSTAL_IN_ETA][CRYSTAL_IN_PHI],
0897                                 ap_uint<5> seed_eta,
0898                                 ap_uint<5> seed_phi,
0899                                 ap_uint<2> brems);
0900   clusterInfo getBremsValuesPos(crystal tempX[CRYSTAL_IN_ETA][CRYSTAL_IN_PHI],
0901                                 ap_uint<5> seed_eta,
0902                                 ap_uint<5> seed_phi);
0903   clusterInfo getBremsValuesNeg(crystal tempX[CRYSTAL_IN_ETA][CRYSTAL_IN_PHI],
0904                                 ap_uint<5> seed_eta,
0905                                 ap_uint<5> seed_phi);
0906   clusterInfo getClusterValues(crystal tempX[CRYSTAL_IN_ETA][CRYSTAL_IN_PHI], ap_uint<5> seed_eta, ap_uint<5> seed_phi);
0907   Cluster getClusterFromRegion3x4(crystal temp[CRYSTAL_IN_ETA][CRYSTAL_IN_PHI]);
0908   void stitchClusterOverRegionBoundary(std::vector<Cluster>& cluster_list, int towerEtaUpper, int towerEtaLower, int cc);
0909 
0910   /*******************************************************************/
0911   /* Cluster flags                                                   */
0912   /*******************************************************************/
0913   inline bool passes_iso(float pt, float iso) {
0914     bool is_iso = true;
0915     if (pt > high_pt_threshold)
0916       is_iso = true;
0917     else if (pt < slideIsoPtThreshold) {
0918       if (!((a0_80 - a1_80 * pt) > iso))
0919         is_iso = false;
0920     } else {
0921       if (iso > a0)
0922         is_iso = false;
0923     }
0924     return is_iso;
0925   }
0926 
0927   inline bool passes_looseTkiso(float pt, float iso) {
0928     bool is_iso;
0929     if (pt > high_pt_threshold)
0930       is_iso = true;
0931     else
0932       is_iso = (b0 + b1 * std::exp(-b2 * pt) > iso);
0933     return is_iso;
0934   }
0935 
0936   inline bool passes_ss(float pt, float eta, float ss) {
0937     bool is_ss;
0938     if (pt > high_pt_threshold)
0939       is_ss = true;
0940     else if ((abs(eta) > eta0_loose_ss && abs(eta) < eta1_loose_ss) ||
0941              (abs(eta) > eta2_loose_ss && abs(eta) < eta3_loose_ss))  // temporary adjustment
0942       is_ss = ((loose_ss_offset + c1_ss * std::exp(-c2_ss * pt)) <= ss);
0943     else
0944       is_ss = ((c0_ss + c1_ss * std::exp(-c2_ss * pt)) <= ss);
0945     return is_ss;
0946   }
0947 
0948   inline bool passes_looseTkss(float pt, float eta, float ss) {
0949     bool is_ss;
0950     if (pt > high_pt_threshold)
0951       is_ss = true;
0952     else if ((abs(eta) > eta0_loose_ss && abs(eta) < eta1_loose_ss) ||
0953              (abs(eta) > eta2_loose_ss && abs(eta) < eta3_loose_ss))  // temporary adjustment
0954       is_ss = ((loose_ss_offset - e1_looseTkss * std::exp(-e2_looseTkss * pt)) <= ss);
0955     else
0956       is_ss = ((e0_looseTkss - e1_looseTkss * std::exp(-e2_looseTkss * pt)) <= ss);
0957     return is_ss;
0958   }
0959 
0960   /*******************************************************************/
0961   /* GCT classes.                                                    */
0962   /*******************************************************************/
0963 
0964   class RCTcluster_t {
0965   public:
0966     ap_uint<12> et;
0967     ap_uint<5> towEta;  // goes from 0 to 17 (vs. class Cluster in the RCT emulator)
0968     ap_uint<2> towPhi;
0969     ap_uint<3> crEta;
0970     ap_uint<3> crPhi;
0971 
0972     ap_uint<12> iso;
0973     ap_uint<15> et2x5;
0974     ap_uint<15> et5x5;
0975     bool is_ss;
0976     bool is_looseTkss;
0977     bool is_iso;
0978     bool is_looseTkiso;
0979     ap_uint<2> brems;
0980 
0981     int nGCTCard;
0982   };
0983 
0984   class RCTtower_t {
0985   public:
0986     ap_uint<12> et;
0987     ap_uint<4> hoe;
0988     // For CMSSW outputs, not firmware
0989     ap_uint<12> ecalEt;
0990     ap_uint<12> hcalEt;
0991   };
0992 
0993   class RCTtoGCTfiber_t {
0994   public:
0995     RCTtower_t RCTtowers[N_RCTTOWERS_FIBER];
0996     RCTcluster_t RCTclusters[N_RCTCLUSTERS_FIBER];
0997   };
0998 
0999   class RCTcard_t {
1000   public:
1001     RCTtoGCTfiber_t RCTtoGCTfiber[N_RCTGCT_FIBERS];
1002   };
1003 
1004   class GCTcard_t {
1005   public:
1006     RCTcard_t RCTcardEtaPos[N_RCTCARDS_PHI];
1007     RCTcard_t RCTcardEtaNeg[N_RCTCARDS_PHI];
1008   };
1009 
1010   class GCTcluster_t {
1011   public:
1012     bool isPositiveEta;  // is cluster in positive eta side or not
1013     ap_uint<12> et;
1014     ap_uint<6> towEta;
1015     ap_uint<7> towPhi;
1016     ap_uint<3> crEta;
1017     ap_uint<3> crPhi;
1018     ap_uint<12> iso;
1019 
1020     ap_uint<15> et2x5;
1021     ap_uint<15> et5x5;
1022     bool is_ss;
1023     bool is_looseTkss;
1024     bool is_iso;
1025     bool is_looseTkiso;
1026 
1027     unsigned int hoe;       // not defined
1028     unsigned int hoe_flag;  // not defined
1029     unsigned int fb;        // not defined
1030     unsigned int timing;    // not defined
1031     ap_uint<2>
1032         brems;  // 0 if no brems applied, 1 or 2 if brems applied (one for + direction, one for - direction: check firmware)
1033 
1034     float relIso;  // for analyzer only, not firmware
1035     int nGCTCard;  // for analyzer only, not firmware
1036 
1037     inline float etFloat() const { return ((float)et * ECAL_LSB); }        // Return energy as a float
1038     inline float isoFloat() const { return ((float)iso * ECAL_LSB); }      // Return energy as a float
1039     inline float et2x5Float() const { return ((float)et2x5 * ECAL_LSB); }  // Return energy as a float
1040     inline float et5x5Float() const { return ((float)et5x5 * ECAL_LSB); }  // Return energy as a float
1041     inline float relIsoFloat() const { return relIso; }                    // relIso is a float already
1042 
1043     inline int etInt() const { return et; }
1044     inline int isoInt() const { return iso; }
1045 
1046     inline int standaloneWP() const { return (is_iso && is_ss); }
1047     inline int looseL1TkMatchWP() const { return (is_looseTkiso && is_looseTkss); }
1048     inline int photonWP() const { return 1; }  // NOTE: NO PHOTON WP
1049 
1050     inline int passesShowerShape() const { return is_ss; }
1051 
1052     /*
1053        * Initialize from RCTcluster_t.
1054        */
1055     void initFromRCTCluster(int iRCTcardIndex, bool isPosEta, const RCTcluster_t& rctCluster) {
1056       isPositiveEta = isPosEta;
1057 
1058       et = rctCluster.et;
1059       towEta = rctCluster.towEta;
1060       if (isPositiveEta) {
1061         towPhi = rctCluster.towPhi + (iRCTcardIndex * 4);
1062       } else {
1063         towPhi = (3 - rctCluster.towPhi) + (iRCTcardIndex * 4);
1064       }
1065       crEta = rctCluster.crEta;
1066       if (isPositiveEta) {
1067         crPhi = rctCluster.crPhi;
1068       } else {
1069         crPhi = (4 - rctCluster.crPhi);
1070       }
1071       et2x5 = rctCluster.et2x5;
1072       et5x5 = rctCluster.et5x5;
1073       is_ss = rctCluster.is_ss;
1074       is_looseTkss = rctCluster.is_looseTkss;
1075       iso = 0;                // initialize: no info from RCT, so set it to null
1076       relIso = 0;             // initialize: no info from RCT, so set it to null
1077       is_iso = false;         // initialize: no info from RCT, so set it to false
1078       is_looseTkiso = false;  // initialize: no info from RCT, so set it to false
1079       hoe = 0;                // initialize: no info from RCT, so set it to null
1080       hoe_flag = 0;           // initialize: no info from RCT, so set it to null
1081       fb = 0;                 // initialize: no info from RCT, so set it to null
1082       timing = 0;             // initialize: no info from RCT, so set it to null
1083       brems = rctCluster.brems;
1084       nGCTCard = rctCluster.nGCTCard;
1085     }
1086 
1087     /*
1088        * Get GCT cluster c's iEta (global iEta convention), (0-33*5). Called in realEta().
1089        */
1090     int globalClusteriEta(void) const {
1091       // First get the "iEta/iPhi" in the GCT card. i.e. in the diagram where the barrel is split up into three GCT cards,
1092       // (iEta, iPhi) = (0, 0) is the top left corner of the GCT card.
1093       int iEta_in_gctCard;
1094       if (!isPositiveEta) {
1095         // Negative eta: c.towEta and c.crEta count outwards from the real eta = 0 center line, so to convert to the barrel diagram global iEta
1096         // (global iEta = 0 from LHS of page), do (17*5 - 1) minus the GCT value.
1097         // e.g. If in GCT, a negative card's cluster had iEta = 84, this would be global iEta = 0.
1098         iEta_in_gctCard =
1099             ((N_GCTTOWERS_FIBER * CRYSTALS_IN_TOWER_ETA - 1) - ((towEta * CRYSTALS_IN_TOWER_ETA) + crEta));
1100       } else {
1101         // c.towEta and c.crEta count outwards from the real eta = 0 center line, so for positive
1102         // eta we need to add the 17*5 offset so that positive eta 0+epsilon starts at 17*5.
1103         // e.g. If in GCT, a positive card's cluster had iEta = 0, this would be global iEta = 85.
1104         // e.g. If in GCT, a positive card's cluster had iEta = 84, this would be global iEta = 169.
1105         iEta_in_gctCard = ((N_GCTTOWERS_FIBER * CRYSTALS_IN_TOWER_ETA) + ((towEta * CRYSTALS_IN_TOWER_ETA) + crEta));
1106       }
1107       // Last, convert to the global iEta/iPhi in the barrel region. For eta these two indices are the same (but this is not true for phi).
1108       int iEta_in_barrel = iEta_in_gctCard;
1109       return iEta_in_barrel;
1110     }
1111 
1112     /* 
1113        * Get GCT cluster iPhi (global convention, (0-71*5)). Called in realPhi().
1114        * Use with getPhi_fromCrystaliPhi from Phase2L1RCT.h to convert from GCT cluster to real phi.
1115        * If returnGlobalGCTiPhi is true (Default value) then return the iPhi in the entire GCT barrel. Otherwise
1116        * just return the iPhi in the current GCT card.
1117        */
1118     int globalClusteriPhi(bool returnGlobalGCTiPhi = true) const {
1119       int iPhi_in_gctCard = ((towPhi * CRYSTALS_IN_TOWER_PHI) + crPhi);
1120       // If we should return the global GCT iPhi, get the iPhi offset due to the number of the GCT card
1121       int iPhi_card_offset = 0;
1122       if (returnGlobalGCTiPhi) {
1123         if (nGCTCard == 0)
1124           iPhi_card_offset = GCTCARD_0_TOWER_IPHI_OFFSET * CRYSTALS_IN_TOWER_PHI;
1125         else if (nGCTCard == 1)
1126           iPhi_card_offset = GCTCARD_1_TOWER_IPHI_OFFSET * CRYSTALS_IN_TOWER_PHI;
1127         else if (nGCTCard == 2)
1128           iPhi_card_offset = GCTCARD_2_TOWER_IPHI_OFFSET * CRYSTALS_IN_TOWER_PHI;
1129       }
1130       // Detector wraps around in phi: modulo number of crystals in phi (n_towers_Phi = 72)
1131       int iPhi_in_barrel = (iPhi_card_offset + iPhi_in_gctCard) % (n_towers_Phi * CRYSTALS_IN_TOWER_PHI);
1132       return iPhi_in_barrel;
1133     }
1134 
1135     /*
1136        * Each cluster falls in a tower: get this tower iEta in the GCT card (same as global) given the cluster's info. 
1137        */
1138     int globalToweriEta(void) const { return (int)(globalClusteriEta() / 5); }
1139 
1140     /*
1141       * Each cluster falls in a tower: get this tower iPhi in global given the cluster's info. 
1142       */
1143     int globalToweriPhi(void) const {
1144       bool getGlobalIndex = true;
1145       return (int)(globalClusteriPhi(getGlobalIndex) / 5);
1146     }
1147 
1148     /*
1149        * Get tower iPhi IN GCT CARD
1150        */
1151     int inCardToweriPhi(void) const {
1152       bool getGlobalIndex = false;
1153       return (int)(globalClusteriPhi(getGlobalIndex) / 5);
1154     }
1155 
1156     /*
1157        * Get tower iEta IN GCT CARD (conveniently the same as global eta)
1158        */
1159     int inCardToweriEta(void) const { return (int)(globalClusteriEta() / 5); }
1160 
1161     /*
1162        * Get GCT cluster's real eta from global iEta (0-33*5).
1163        */
1164     float realEta(void) const {
1165       float size_cell = 2 * ECAL_eta_range / (CRYSTALS_IN_TOWER_ETA * n_towers_Eta);
1166       return globalClusteriEta() * size_cell - ECAL_eta_range + half_crystal_size;
1167     }
1168 
1169     /* 
1170        * Get GCT cluster's real eta from global iPhi (0-71*5).
1171        */
1172     float realPhi(void) const {
1173       float size_cell = 2 * M_PI / (CRYSTALS_IN_TOWER_PHI * n_towers_Phi);
1174       return globalClusteriPhi() * size_cell - M_PI + half_crystal_size;
1175     }
1176 
1177     /* 
1178        * Return the 4-vector.
1179        */
1180     reco::Candidate::PolarLorentzVector p4(void) const {
1181       return reco::Candidate::PolarLorentzVector(etFloat(), realEta(), realPhi(), 0.);
1182     }
1183 
1184     /*
1185        *  Compute relative isolation and set its flags in-place, assuming that the isolation is already computed.
1186        */
1187     void setRelIsoAndFlags(void) {
1188       float relativeIsolationAsFloat = 0;
1189       if (et > 0) {
1190         relativeIsolationAsFloat = (isoFloat() / etFloat());
1191       } else {
1192         relativeIsolationAsFloat = 0;
1193       }
1194       relIso = relativeIsolationAsFloat;
1195       is_iso = passes_iso(etFloat(), relIso);
1196       is_looseTkiso = passes_looseTkiso(isoFloat(), relIso);
1197     }
1198 
1199     /*
1200        * Create a l1tp2::CaloCrystalCluster object.
1201        */
1202     l1tp2::CaloCrystalCluster createCaloCrystalCluster(void) const {
1203       l1tp2::CaloCrystalCluster caloCrystalCluster(
1204           p4(),
1205           etFloat(),      // convert to float
1206           0,              // supposed to be H over E in the constructor but we do not store/use this information
1207           relIsoFloat(),  // for consistency with the old emulator, in this field save (iso energy sum)/(cluster energy)
1208           0,              // DetId seedCrystal
1209           0,              // puCorrPt
1210           0,              // brems: not propagated to output (0, 1, or 2 as computed in firmware)
1211           0,              // et2x2 (not calculated)
1212           et2x5Float(),   // et2x5 (save float)
1213           0,              // et3x5 (not calculated)
1214           et5x5Float(),   // et5x5 (save float)
1215           standaloneWP(),  // standalone WP
1216           false,           // electronWP98: not computed
1217           false,           // photonWP80: not computed
1218           false,           // electronWP90: not computed
1219           false,           // looseL1TkMatchWP: not computed
1220           false            // stage2effMatch: not computed
1221       );
1222 
1223       // Flags
1224       std::map<std::string, float> params;
1225       params["standaloneWP_showerShape"] = is_ss;
1226       params["standaloneWP_isolation"] = is_iso;
1227       params["trkMatchWP_showerShape"] = is_looseTkss;
1228       params["trkMatchWP_isolation"] = is_looseTkiso;
1229       caloCrystalCluster.setExperimentalParams(params);
1230 
1231       return caloCrystalCluster;
1232     }
1233 
1234     /*
1235       * Create a l1t::EGamma object.
1236       */
1237     l1t::EGamma createL1TEGamma(void) const {
1238       // n.b. No photon WP, photonWP() always returns true
1239       int quality =
1240           (standaloneWP() * std::pow(2, 0)) + (looseL1TkMatchWP() * std::pow(2, 1)) + (photonWP() * std::pow(2, 2));
1241 
1242       // The constructor zeros out everyhing except the p4()
1243       l1t::EGamma eg = l1t::EGamma(p4(), etInt(), globalClusteriEta(), globalClusteriPhi(), quality, isoInt());
1244 
1245       // Write in fields that were zerod out
1246       eg.setRawEt(etInt());                // et as int
1247       eg.setTowerIEta(globalToweriEta());  // 0-33 in barrel
1248       eg.setTowerIPhi(globalToweriPhi());  // 0-71 in barrel
1249       eg.setIsoEt(isoInt());               // raw isolation sum as int
1250       eg.setShape(passesShowerShape());    // write shower shape flag to this field
1251       return eg;
1252     }
1253 
1254     /*
1255      * Create a l1tp2::DigitizedClusterCorrelator object, with corrTowPhiOffset specifying the offset necessary to correct the tower phi to the region
1256      * unique to each GCT card.
1257      */
1258     l1tp2::DigitizedClusterCorrelator createDigitizedClusterCorrelator(const int corrTowPhiOffset) const {
1259       return l1tp2::DigitizedClusterCorrelator(
1260           etFloat(),  // technically we are just multiplying and then dividing again by the LSB
1261           globalClusteriEta(),
1262           ((towPhi - corrTowPhiOffset) * CRYSTALS_IN_TOWER_PHI) +
1263               crPhi,  // cannot use globalClusteriPhi() helper function because correlator offset is different than GCT offset
1264           hoe,
1265           hoe_flag,
1266           iso,
1267           (is_iso) | (is_looseTkiso << 1),  // 2 bits: e.g. 0b10 means is_looseTkiso was true, and is_iso was false
1268           fb,
1269           timing,
1270           (is_ss) | (is_looseTkss << 1),  // 2 bits (same as iso flags) is_ss in lowest bit, is_looseTkss in higher bit
1271           brems,
1272           nGCTCard);
1273     }
1274 
1275     /*
1276      * Create a l1tp2::DigitizedClusterGT object
1277      */
1278     l1tp2::DigitizedClusterGT createDigitizedClusterGT(bool isValid) const {
1279       // Constructor arguments take phi, then eta
1280       return l1tp2::DigitizedClusterGT(isValid, etFloat(), realPhi(), realEta());
1281     }
1282 
1283     /*
1284       * Print GCT cluster information.
1285       */
1286     void printGCTClusterInfo(std::string description = "") {
1287       std::cout << "[PrintGCTClusterInfo:] [" << description << "]: "
1288                 << "et (float): " << etFloat() << ", "
1289                 << "eta: " << realEta() << ", "
1290                 << "phi: " << realPhi() << ", "
1291                 << "isPositiveEta " << isPositiveEta << ", "
1292                 << "towEta: " << towEta << ", "
1293                 << "towPhi: " << towPhi << ", "
1294                 << "crEta: " << crEta << ", "
1295                 << "crPhi: " << crPhi << ", "
1296                 << "iso (GeV): " << isoFloat() << ", "
1297                 << "rel iso (unitless float): " << relIsoFloat() << ", "
1298                 << "et2x5 (GeV): " << et2x5Float() << ", "
1299                 << "et5x5 (GeV): " << et5x5Float() << ", "
1300                 << "is_ss: " << is_ss << ", "
1301                 << "is_looseTkss" << is_looseTkss << ", "
1302                 << "is_iso: " << is_iso << ", "
1303                 << "is_looseTkiso: " << is_looseTkiso << ", "
1304                 << "brems: " << brems << std::endl;
1305     }
1306   };
1307 
1308   class GCTtower_t {
1309   public:
1310     ap_uint<12> et;
1311     ap_uint<4> hoe;
1312     ap_uint<2> fb;  // not defined yet in emulator
1313     // For CMSSW outputs, not firmware
1314     ap_uint<12> ecalEt;
1315     ap_uint<12> hcalEt;
1316 
1317     inline float totalEtFloat() const {
1318       return ((float)et * ECAL_LSB);
1319     }  // Return total energy as a float (assuming the energy uses the ECAL LSB convention)
1320     inline float ecalEtFloat() const { return ((float)ecalEt * ECAL_LSB); }  // Return ECAL energy as a float
1321     inline float hcalEtFloat() const {
1322       return ((float)hcalEt * HCAL_LSB);
1323     }  // Return HCAL energy as a float, use HCAL LSB
1324 
1325     /*
1326        * Initialize from RCTtower_t.
1327        */
1328     void initFromRCTTower(const RCTtower_t& rctTower) {
1329       et = rctTower.et;
1330       hoe = rctTower.hoe;
1331       ecalEt = rctTower.ecalEt;
1332       hcalEt = rctTower.hcalEt;
1333     }
1334 
1335     /*
1336       * Correlator fiber convention -> Global GCT convention
1337       * Get tower's global (iEta) from the GCTCorrFiber index [0, 64) and the tower's postion in the fiber [0, 17).
1338       * Recall that GCTCorrFiber is [0, 32) for negative eta and [32, 64) for positive eta. The tower's position in the fiber [0, 17)
1339       * always counts outwards from real eta = 0.
1340       * Use in conjunction with (float) getTowerEta_fromAbsID(int id) from Phase2L1RCT.h to get a tower's real eta.
1341       */
1342     int globalToweriEta(unsigned int nGCTCard, unsigned int gctCorrFiberIdx, unsigned int posInFiber) {
1343       (void)nGCTCard;                                                        // not needed
1344       bool isTowerInPositiveEta = (gctCorrFiberIdx < N_GCTPOSITIVE_FIBERS);  // N_GCTPOSITIVE_FIBERS = 32
1345       int global_toweriEta;
1346       if (isTowerInPositiveEta) {
1347         global_toweriEta = (N_GCTTOWERS_FIBER + posInFiber);  // N_GCTTOWERS_FIBER = 17
1348       } else {
1349         // e.g. For negative eta, posInFiber = 0 is at real eta = 0, and global tower iEta is 17 - 1 - 0 = 16
1350         // posInFiber = 16 is at real eta = -1.4841, and global tower iEta is 17 - 1 - 16 = 0.
1351         global_toweriEta = (N_GCTTOWERS_FIBER - 1 - posInFiber);
1352       }
1353       return global_toweriEta;
1354     }
1355 
1356     /*
1357        * Correlator fiber convention -> Global GCT convention
1358       * Get tower's global (iPhi) from the GCT card number (0, 1, 2), and the GCTCorrFiber index [0, 64).
1359       * GCTCorrFiber is [0, 32) for negative eta and [32, 64) for positive eta. In the phi direction, fiber index #0 has the same phi 
1360       * as fiber index #32, so only the (fiber index modulo 32) matters for the phi direction. 
1361       * The tower's position in the fiber doesn't matter; in each fiber the phi is the same. 
1362       * Use in conjunction with (float) getTowerPhi_fromAbsID(int id) from Phase2L1RCT.h to get a tower's real phi.
1363       */
1364     int globalToweriPhi(unsigned int nGCTCard, unsigned int gctCorrFiberIdx, unsigned int posInFiber) {
1365       (void)posInFiber;                                                           // not needed
1366       unsigned int effectiveFiberIdx = (gctCorrFiberIdx % N_GCTPOSITIVE_FIBERS);  // N_GCTPOSITIVE_FIBERS = 32
1367       int toweriPhi_card_offset = 0;
1368       if (nGCTCard == 0)
1369         toweriPhi_card_offset = GCTCARD_0_TOWER_IPHI_OFFSET;
1370       else if (nGCTCard == 1)
1371         toweriPhi_card_offset = GCTCARD_1_TOWER_IPHI_OFFSET;
1372       else if (nGCTCard == 2)
1373         toweriPhi_card_offset = GCTCARD_2_TOWER_IPHI_OFFSET;
1374 
1375       //  as explained above, effectiveFiberIdx is [0, 32). n_towers_Phi = 72:
1376       int global_tower_iPhi = (toweriPhi_card_offset + effectiveFiberIdx) % (n_towers_Phi);
1377       return global_tower_iPhi;
1378     }
1379 
1380     /*
1381      * For fulltowers that are indexed by GCT local index: eta
1382      */
1383     int globalToweriEtaFromGCTcardiEta(int gctCard_tower_iEta) {
1384       int global_iEta = gctCard_tower_iEta;
1385       return global_iEta;
1386     }
1387 
1388     /*
1389      * For fulltowers that are indexed by GCT local index: phi. Very similar to globalToweriPhi function but keep them separate for clarity.
1390      */
1391     int globalToweriPhiFromGCTcardiPhi(unsigned int nGCTCard, int gctCard_tower_iPhi) {
1392       assert(nGCTCard <= 2);  // Make sure the card number is valid
1393       int toweriPhi_card_offset = 0;
1394       if (nGCTCard == 0)
1395         toweriPhi_card_offset = GCTCARD_0_TOWER_IPHI_OFFSET;
1396       else if (nGCTCard == 1)
1397         toweriPhi_card_offset = GCTCARD_1_TOWER_IPHI_OFFSET;
1398       else if (nGCTCard == 2)
1399         toweriPhi_card_offset = GCTCARD_2_TOWER_IPHI_OFFSET;
1400 
1401       int global_iPhi = (toweriPhi_card_offset + gctCard_tower_iPhi) % (n_towers_Phi);  //   n_towers_Phi = 72
1402       return global_iPhi;
1403     }
1404 
1405     /* 
1406        * Method to create a l1tp2::CaloTower object from the fiber and tower-in-fiber indices. 
1407        * nGCTCard (0, 1, 2) is needed to determine the absolute eta/phi.
1408        * iFiber and iTowerInFiber are the indices of the tower in the card, e.g. GCTinternal.GCTCorrFiber[iFiber].GCTtowers[iTowerInFiber]
1409        */
1410     l1tp2::CaloTower createCaloTowerFromFiberIdx(int nGCTCard, int iFiber, int iTowerInFiber) {
1411       l1tp2::CaloTower l1CaloTower;
1412       l1CaloTower.setEcalTowerEt(ecalEtFloat());  // float: ECAL divide by 8.0
1413       l1CaloTower.setHcalTowerEt(hcalEtFloat());  // float: HCAL multiply by LSB
1414       int global_tower_iEta = globalToweriEta(nGCTCard, iFiber, iTowerInFiber);
1415       int global_tower_iPhi = globalToweriPhi(nGCTCard, iFiber, iTowerInFiber);
1416       l1CaloTower.setTowerIEta(global_tower_iEta);
1417       l1CaloTower.setTowerIPhi(global_tower_iPhi);
1418       l1CaloTower.setTowerEta(getTowerEta_fromAbsID(global_tower_iEta));
1419       l1CaloTower.setTowerPhi(getTowerPhi_fromAbsID(global_tower_iPhi));
1420       return l1CaloTower;
1421     }
1422 
1423     /*
1424      * Method to create a l1tp2::CaloTower object from the global tower ieta and iphi.
1425      */
1426     l1tp2::CaloTower createFullTowerFromCardIdx(int nGCTCard, int gctCard_tower_iEta, int gctCard_tower_iPhi) {
1427       l1tp2::CaloTower l1CaloTower;
1428       // Store total Et (HCAL+ECAL) in the ECAL Et member
1429       l1CaloTower.setEcalTowerEt(totalEtFloat());
1430       l1CaloTower.setHcalTowerEt(ecalEtFloat());
1431       int global_tower_iEta = globalToweriEtaFromGCTcardiEta(gctCard_tower_iEta);
1432       int global_tower_iPhi = globalToweriPhiFromGCTcardiPhi(nGCTCard, gctCard_tower_iPhi);
1433       l1CaloTower.setTowerIEta(global_tower_iEta);
1434       l1CaloTower.setTowerIPhi(global_tower_iPhi);
1435       l1CaloTower.setTowerEta(getTowerEta_fromAbsID(global_tower_iEta));
1436       l1CaloTower.setTowerPhi(getTowerPhi_fromAbsID(global_tower_iPhi));
1437       return l1CaloTower;
1438     }
1439 
1440     /*
1441      * Method to create a l1tp2::DigitizedTowerCorrelator, from the GCT card number, the fiber index *inside the GCT card* (excluding overlap region),
1442      * and the index of the tower inside the fiber.
1443      */
1444     l1tp2::DigitizedTowerCorrelator createDigitizedTowerCorrelator(unsigned int indexCard,
1445                                                                    unsigned int indexFiber,
1446                                                                    unsigned int indexTower) {
1447       return l1tp2::DigitizedTowerCorrelator(totalEtFloat(), hoe, fb, indexCard, indexFiber, indexTower);
1448     }
1449 
1450     /*
1451      * Print GCTtower_t tower information.
1452      */
1453     void printGCTTowerInfoFromGlobalIdx(int global_tower_iEta, int global_tower_iPhi, std::string description = "") {
1454       std::cout << "[Print GCTtower_t class info from global idx:] [" << description << "]: "
1455                 << "total et (float): " << totalEtFloat() << ", "
1456                 << "ecal et (float): " << ecalEtFloat() << ", "
1457                 << "hcal et (float): " << hcalEtFloat() << ", "
1458                 << "fb: " << fb << ", "
1459                 << "global tower ieta: " << global_tower_iEta << ", "
1460                 << "global tower iphi: " << global_tower_iPhi << ", "
1461                 << "eta: " << getTowerEta_fromAbsID(global_tower_iEta) << ", "
1462                 << "phi: " << getTowerPhi_fromAbsID(global_tower_iPhi) << std::endl;
1463     }
1464   };
1465 
1466   class GCTCorrfiber_t {
1467   public:
1468     GCTtower_t GCTtowers[N_GCTTOWERS_FIBER];
1469     GCTcluster_t GCTclusters[N_GCTCLUSTERS_FIBER];
1470   };
1471 
1472   class GCTtoCorr_t {
1473   public:
1474     GCTCorrfiber_t GCTCorrfiber[N_GCTCORR_FIBERS];
1475   };
1476 
1477   class GCTinternal_t {
1478   public:
1479     GCTCorrfiber_t GCTCorrfiber[N_GCTINTERNAL_FIBERS];
1480 
1481     void computeClusterIsolationInPlace(int nGCTCard) {
1482       for (unsigned int iFiber = 0; iFiber < N_GCTINTERNAL_FIBERS; iFiber++) {
1483         for (unsigned int iCluster = 0; iCluster < N_GCTCLUSTERS_FIBER; iCluster++) {
1484           // We will only save clusters with > 0 GeV, so only need to do this for clusters with >0 energy
1485           if (GCTCorrfiber[iFiber].GCTclusters[iCluster].et == 0) {
1486             GCTCorrfiber[iFiber].GCTclusters[iCluster].iso = 0;
1487             continue;
1488           }
1489 
1490           ap_uint<12> uint_isolation = 0;
1491 
1492           // do not add the GCT card off-set, so we remain in the gct local card iEta/iPhi
1493           int toweriEta_in_GCT_card = GCTCorrfiber[iFiber].GCTclusters[iCluster].inCardToweriEta();
1494           int toweriPhi_in_GCT_card = GCTCorrfiber[iFiber].GCTclusters[iCluster].inCardToweriPhi();
1495 
1496           // If cluster is in the overlap region, do not compute isolation
1497           bool inOverlapWithAnotherGCTCard = (((toweriPhi_in_GCT_card >= 0) && (toweriPhi_in_GCT_card < 4)) ||
1498                                               ((toweriPhi_in_GCT_card >= 28) && (toweriPhi_in_GCT_card < 32)));
1499           if (inOverlapWithAnotherGCTCard) {
1500             GCTCorrfiber[iFiber].GCTclusters[iCluster].iso = 0;
1501             continue;
1502           }
1503 
1504           // Size 5x5 in towers: include the overlap-region-between-GCT-cards-if-applicable. In eta direction, the min and max towers (inclusive!) are:
1505           int isoWindow_toweriEta_in_GCT_card_min = std::max(0, toweriEta_in_GCT_card - 2);
1506           int isoWindow_toweriEta_in_GCT_card_max = std::min(toweriEta_in_GCT_card + 2, N_GCTETA - 1);  // N_GCTETA = 34
1507           // e.g. if our window is centered at tower_iEta = 5, we want to sum towers_iEta 3, 4, (5), 6, 7, inclusive
1508           // e.g. if our window is near the boundary, tower_iEta = 32, we want to sum towers_iEta 30, 31, (32), 33
1509           // inclusive (but there are only N_GCTETA = 34 towers, so we stop at tower_iEta = 33)
1510 
1511           // in phi direction, the min and max towers (inclusive!) are:
1512           int isoWindow_toweriPhi_in_GCT_card_min = std::max(0, toweriPhi_in_GCT_card - 2);
1513           int isoWindow_toweriPhi_in_GCT_card_max = std::min(toweriPhi_in_GCT_card + 2, N_GCTPHI - 1);
1514 
1515           // Keep track of the number of towers we summed over
1516           int nTowersSummed = 0;
1517 
1518           // First add any nearby clusters to the isolation
1519           for (unsigned int candFiber = 0; candFiber < N_GCTINTERNAL_FIBERS; candFiber++) {
1520             for (unsigned int candCluster = 0; candCluster < N_GCTCLUSTERS_FIBER; candCluster++) {
1521               // Do not double-count the cluster we are calculating the isolation for
1522               if (!((candFiber == iFiber) && (candCluster == iCluster))) {
1523                 // Only consider clusters with et > 0 for isolation sum
1524                 if (GCTCorrfiber[candFiber].GCTclusters[candCluster].et > 0) {
1525                   // Get the candidate cluster's tower iEta and iPhi in GCT card
1526                   int candidate_toweriEta = GCTCorrfiber[candFiber].GCTclusters[candCluster].inCardToweriEta();
1527                   int candidate_toweriPhi = GCTCorrfiber[candFiber].GCTclusters[candCluster].inCardToweriPhi();
1528 
1529                   // If the tower that the candidate cluster is in, is within a 5x5 window, add the candidate cluster energy's to the isolation as a proxy for the ECAL energy
1530                   if (((candidate_toweriEta >= isoWindow_toweriEta_in_GCT_card_min) &&
1531                        (candidate_toweriEta <= isoWindow_toweriEta_in_GCT_card_max)) &&
1532                       ((candidate_toweriPhi >= isoWindow_toweriPhi_in_GCT_card_min) &&
1533                        (candidate_toweriPhi <= isoWindow_toweriPhi_in_GCT_card_max))) {
1534                     uint_isolation += GCTCorrfiber[candFiber].GCTclusters[candCluster].et;
1535                   }
1536                 }
1537               }
1538             }
1539           }
1540 
1541           //  From "tower index in GCT card", get which fiber it is in (out of 64 fibers), and which tower it is inside the fiber (out of 17 towers)
1542           for (int iEta = isoWindow_toweriEta_in_GCT_card_min; iEta <= isoWindow_toweriEta_in_GCT_card_max; iEta++) {
1543             for (int iPhi = isoWindow_toweriPhi_in_GCT_card_min; iPhi <= isoWindow_toweriPhi_in_GCT_card_max; iPhi++) {
1544               nTowersSummed += 1;
1545 
1546               int indexInto64Fibers;
1547               int indexInto17TowersInFiber;
1548 
1549               bool isTowerInPositiveEta = (iEta >= N_GCTTOWERS_FIBER);
1550               if (isTowerInPositiveEta) {
1551                 // phi index is simple (e.g. if real phi = +80 degrees, iPhi in GCT = 31)
1552                 indexInto64Fibers = iPhi;
1553                 // if real eta = 1.47, iEta in GCT card = 33. If real eta = 0.0, iEta in GCT = 17, so iEta in fiber = 17%17 = 0.
1554                 indexInto17TowersInFiber = (iEta % 17);
1555               } else {
1556                 // add offset (e.g. if real phi = +80 degrees, iPhi in GCT = 31, and my index into GCT fibers 31 + 32 = 63)
1557                 indexInto64Fibers = (iPhi + N_GCTPOSITIVE_FIBERS);
1558                 // e.g.  if real eta = 0, iEta innew GCT card = 16, i.e. our index into the GCT fiber is 16-16 = 0
1559                 indexInto17TowersInFiber = (16 - iEta);
1560               }
1561 
1562               ap_uint<12> ecalEt = GCTCorrfiber[indexInto64Fibers].GCTtowers[indexInto17TowersInFiber].ecalEt;
1563               uint_isolation += ecalEt;
1564             }
1565           }
1566 
1567           // Scale the isolation sum up if we summed over fewer than (5x5) = 25 towers
1568           float scaleFactor =
1569               ((float)(N_GCTTOWERS_CLUSTER_ISO_ONESIDE * N_GCTTOWERS_CLUSTER_ISO_ONESIDE) / (float)nTowersSummed);
1570 
1571           uint_isolation = (ap_uint<12>)(((float)uint_isolation) * scaleFactor);
1572 
1573           // Set the iso in the cluster
1574           GCTCorrfiber[iFiber].GCTclusters[iCluster].iso = uint_isolation;
1575         }
1576       }
1577     }
1578 
1579     void setIsolationInfo(void) {
1580       for (unsigned int iFiber = 0; iFiber < N_GCTINTERNAL_FIBERS; iFiber++) {
1581         for (unsigned int iCluster = 0; iCluster < N_GCTCLUSTERS_FIBER; iCluster++) {
1582           // update the cluster's isolation information
1583           GCTCorrfiber[iFiber].GCTclusters[iCluster].setRelIsoAndFlags();
1584         }
1585       }
1586     }
1587   };
1588 
1589   class GCTintTowers_t {
1590   public:
1591     GCTtower_t GCTtower[N_GCTETA][N_GCTPHI];
1592 
1593     // Write contents to output CMSSW collection. Note the use of the GCTtower_t method that creates the
1594     // l1tp2::CaloTower object from the global eta/phi.
1595     void writeToPFOutput(int nGCTCard, std::unique_ptr<l1tp2::CaloTowerCollection> const& gctFullTowers) {
1596       for (unsigned int iEta = 0; iEta < N_GCTETA; iEta++) {
1597         for (unsigned int iPhi = 0; iPhi < N_GCTPHI; iPhi++) {
1598           GCTtower_t thisFullTower = GCTtower[iEta][iPhi];
1599           gctFullTowers->push_back(thisFullTower.createFullTowerFromCardIdx(nGCTCard, iEta, iPhi));
1600         }
1601       }
1602     }
1603   };
1604 
1605   /* For each GCT card (3 of them in total, for barrel + endcap), list the sixteen                
1606     * RCT cards that fall in them. The first eight are in positive eta, the next                   
1607     * eight are in negative eta (see figure of one GCT card). The RCT cards themselves             
1608     * run from 0 to 35 (see RCT figure).                                                          
1609     * Hard-coded because the GCT cards wrap around the barrel region.                            
1610     * Used only to convert the RCT emulator outputs to the GCT emulator inputs.                   
1611     */
1612   static const unsigned int GCTcardtoRCTcardnumber[N_GCTCARDS][N_RCTCARDS_PHI * 2] = {
1613       // GCT Card 0
1614       {11, 13, 15, 17, 19, 21, 23, 25, 10, 12, 14, 16, 18, 20, 22, 24},
1615 
1616       // GCT Card 1
1617       {23, 25, 27, 29, 31, 33, 35, 1, 22, 24, 26, 28, 30, 32, 34, 0},
1618 
1619       // GCT Card 2
1620       {35, 1, 3, 5, 7, 9, 11, 13, 34, 0, 2, 4, 6, 8, 10, 12}};
1621 
1622   /*
1623    * Helper function to monitor l1tp2::CaloTower members.
1624    */
1625   inline void printl1tp2TowerInfo(l1tp2::CaloTower thisTower, std::string description = "") {
1626     std::cout << "[Print l1tp2::CaloTower info:] [" << description << "]: "
1627               << ".ecalTowerEta() (float): " << thisTower.ecalTowerEt() << ", "
1628               << ".hcalTowerEta() (float): " << thisTower.hcalTowerEt() << ", "
1629               << ".towerIEta(): " << thisTower.towerIEta() << ", "
1630               << ".towerIPhi(): " << thisTower.towerIPhi() << ", "
1631               << ".towerEta() " << thisTower.towerEta() << ", "
1632               << ".towerPhi() " << thisTower.towerPhi() << std::endl;
1633   }
1634 
1635   void algo_top(const GCTcard_t& GCTcard,
1636                 GCTtoCorr_t& GCTtoCorr,
1637                 unsigned int nGCTCard,
1638                 std::unique_ptr<l1tp2::CaloCrystalClusterCollection> const& gctClusters,
1639                 std::unique_ptr<l1tp2::CaloTowerCollection> const& gctTowers,
1640                 std::unique_ptr<l1tp2::CaloTowerCollection> const& gctFullTowers,
1641                 std::unique_ptr<l1t::EGammaBxCollection> const& gctEGammas,
1642                 std::unique_ptr<l1tp2::DigitizedClusterCorrelatorCollection> const& gctDigitizedClustersCorrelator,
1643                 std::unique_ptr<l1tp2::DigitizedTowerCorrelatorCollection> const& gctDigitizedTowersCorrelator,
1644                 std::unique_ptr<l1tp2::DigitizedClusterGTCollection> const& gctDigitizedClustersGT,
1645                 l1tp2::ParametricCalibration calib_);
1646 
1647   GCTinternal_t getClustersTowers(const GCTcard_t& GCTcard, unsigned int nGCTCard);
1648 
1649   void doProximityAndBremsStitching(const RCTcard_t (&inputCards)[N_RCTCARDS_PHI],
1650                                     RCTcard_t (&outputCards)[N_RCTCARDS_PHI],
1651                                     int iStartingCard,
1652                                     bool isPositiveEta);
1653 
1654   GCTcard_t getClustersCombined(const GCTcard_t& GCTcard, unsigned int nGCTCard);
1655 
1656   GCTintTowers_t getFullTowers(const GCTinternal_t& GCTinternal);
1657 
1658   void writeToCorrelatorAndGTOutputs(
1659       const GCTinternal_t& GCTinternal,
1660       GCTtoCorr_t& GCTtoCorrOutput,
1661       std::unique_ptr<l1tp2::CaloCrystalClusterCollection> const& gctClustersOutput,
1662       std::unique_ptr<l1tp2::CaloTowerCollection> const& gctTowersOutput,
1663       std::unique_ptr<l1t::EGammaBxCollection> const& gctEGammas,
1664       std::unique_ptr<l1tp2::DigitizedClusterCorrelatorCollection> const& gctDigitizedClustersCorrelator,
1665       std::unique_ptr<l1tp2::DigitizedTowerCorrelatorCollection> const& gctDigitizedTowersCorrelator,
1666       std::unique_ptr<l1tp2::DigitizedClusterGTCollection> const& gctDigitizedClustersGT,
1667       int nGCTCard,
1668       int fiberStart,
1669       int fiberEnd,
1670       int corrFiberIndexOffset,
1671       int corrTowPhiOffset);
1672 
1673   /*******************************************************************************************/
1674   /* Interface to correlator helper functions (defined in Phase2L1CaloBarrelToCorrelator.h)  */
1675   /*******************************************************************************************/
1676   float deltaPhiInDegrees(float phi1, float phi2, const float c);
1677   float wrappedPhiInDegrees(float phi);
1678 
1679   /*
1680    * Generic function to compare hadronic (had) or EM digi clusters.
1681    */
1682   template <typename T>
1683   bool compareGCTDigiClusterET(const T& lhs, const T& rhs) {
1684     return (lhs.ptFloat() > rhs.ptFloat());
1685   }
1686 
1687   /*
1688   * Generic function to sort a vector of hadronic or EM digi clusters in each SLR in descending pT,
1689   * then if there are fewer than nMaxClustersInSLR clusters in the vector, pad the vector with null clusters,
1690   * or truncate the vector if there are more than nMaxClustersInSLR clusters.
1691   */
1692   template <typename T>
1693   void sortAndPadSLR(std::vector<T>& thisSLR, const int nMaxClustersInSLR) {
1694     // input is a vector and can be sorted
1695     std::sort(thisSLR.begin(), thisSLR.end(), compareGCTDigiClusterET<T>);
1696     int nClusters = thisSLR.size();
1697 
1698     // If there are fewer than the designated number of clusters, pad with zeros
1699     if (nClusters < nMaxClustersInSLR) {
1700       for (int i = 0; i < (nMaxClustersInSLR - nClusters); i++) {
1701         T zeroCluster;
1702         thisSLR.push_back(zeroCluster);
1703       }
1704     }
1705     // If there are more than the designated number of clusters, truncate the vector
1706     else if (nClusters > nMaxClustersInSLR) {
1707       // Get the iterator to the sixth element and delete til the end of the vector
1708       thisSLR.erase(thisSLR.begin() + nMaxClustersInSLR, thisSLR.end());
1709     }
1710   }
1711 
1712 }  // namespace p2eg
1713 
1714 #endif