Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:20:07

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