Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-12-05 03:12:11

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