Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-11 23:28:02

0001 #ifndef L1TRIGGER_L1CALOTRIGGER_PHASE2L1CALOJETEMULATOR_H
0002 #define L1TRIGGER_L1CALOTRIGGER_PHASE2L1CALOJETEMULATOR_H
0003 
0004 #include "DataFormats/L1TCalorimeterPhase2/interface/Phase2L1CaloJet.h"
0005 #include <cstdlib>
0006 
0007 static constexpr int nBarrelEta = 34;
0008 static constexpr int nBarrelPhi = 72;
0009 static constexpr int nHgcalEta = 36;
0010 static constexpr int nHgcalPhi = 72;
0011 static constexpr int nHfEta = 24;
0012 static constexpr int nHfPhi = 72;
0013 static constexpr int nSTEta = 6;
0014 static constexpr int nSTEta_hf = 4;
0015 static constexpr int nSTPhi = 24;
0016 static constexpr int nJets = 10;
0017 
0018 namespace gctobj {
0019 
0020   class towerMax {
0021   public:
0022     float energy;
0023     int phi;
0024     int eta;
0025     float energyMax;
0026     int phiMax;
0027     int etaMax;
0028     int phiCenter;
0029     int etaCenter;
0030 
0031     towerMax() {
0032       energy = 0;
0033       phi = 0;
0034       eta = 0;
0035       energyMax = 0;
0036       phiMax = 0;
0037       etaMax = 0;
0038       phiCenter = 0;
0039       etaCenter = 0;
0040     }
0041   };
0042 
0043   class jetInfo {
0044   public:
0045     float seedEnergy;
0046     float energy;
0047     float tauEt;
0048     int phi;
0049     int eta;
0050     float energyMax;
0051     int phiMax;
0052     int etaMax;
0053     int phiCenter;
0054     int etaCenter;
0055 
0056     jetInfo() {
0057       seedEnergy = 0;
0058       energy = 0;
0059       tauEt = 0;
0060       phi = 0;
0061       eta = 0;
0062       energyMax = 0;
0063       phiMax = 0;
0064       etaMax = 0;
0065       phiCenter = 0;
0066       etaCenter = 0;
0067     }
0068   };
0069 
0070   typedef struct {
0071     float et;
0072     int eta;
0073     int phi;
0074     float towerEt;
0075     int towerEta;
0076     int towerPhi;
0077     int centerEta;
0078     int centerPhi;
0079   } GCTsupertower_t;
0080 
0081   typedef struct {
0082     GCTsupertower_t cr[nSTPhi];
0083   } etaStrip_t;
0084 
0085   typedef struct {
0086     GCTsupertower_t etaStrip[nSTEta];
0087   } hgcalRegion_t;
0088 
0089   typedef struct {
0090     GCTsupertower_t pk[nSTEta];
0091   } etaStripPeak_t;
0092 
0093   typedef struct {
0094     float et;
0095     float hoe;
0096   } GCTtower_t;
0097 
0098   typedef struct {
0099     GCTtower_t GCTtower[nBarrelEta / 2][nBarrelPhi];
0100   } GCTintTowers_t;
0101 
0102   inline int getPeakBinOf3(float et0, float et1, float et2) {
0103     int x;
0104     float temp;
0105     if (et0 > et1) {
0106       x = 0;
0107       temp = et0;
0108     } else {
0109       x = 1;
0110       temp = et1;
0111     }
0112     if (et2 > temp) {
0113       x = 2;
0114     }
0115     return x;
0116   }
0117 
0118   inline int getEtCenterOf3(float et0, float et1, float et2) {
0119     float etSum = et0 + et1 + et2;
0120     float iEtSum = 0.5 * et0 + 1.5 * et1 + 2.5 * et2;
0121     int iAve = 0xEEF;
0122     if (iEtSum <= etSum)
0123       iAve = 0;
0124     else if (iEtSum <= 2 * etSum)
0125       iAve = 1;
0126     else
0127       iAve = 2;
0128     return iAve;
0129   }
0130 
0131   inline void makeST(const float GCTintTowers[nBarrelEta / 2][nBarrelPhi],
0132                      GCTsupertower_t supertower_return[nSTEta][nSTPhi]) {
0133     float et_sumEta[nSTEta][nSTPhi][3];
0134     float stripEta[nSTEta][nSTPhi][3];
0135     float stripPhi[nSTEta][nSTPhi][3];
0136 
0137     float ex_et[nBarrelEta / 2 + 1][nBarrelPhi];
0138     for (int j = 0; j < nBarrelPhi; j++) {
0139       ex_et[nBarrelEta / 2][j] = 0;
0140       for (int i = 0; i < nBarrelEta / 2; i++) {
0141         ex_et[i][j] = GCTintTowers[i][j];
0142       }
0143     }
0144 
0145     int index_i = 0;
0146     int index_j = 0;
0147     for (int i = 0; i < nBarrelEta / 2 + 1; i += 3) {  // 17+1 to divide into 6 super towers
0148       index_j = 0;
0149       for (int j = 0; j < nBarrelPhi; j += 3) {  // 72 phi to 24 super towers
0150         stripEta[index_i][index_j][0] = ex_et[i][j] + ex_et[i][j + 1] + ex_et[i][j + 2];
0151         stripEta[index_i][index_j][1] = ex_et[i + 1][j] + ex_et[i + 1][j + 1] + ex_et[i + 1][j + 2];
0152         stripEta[index_i][index_j][2] = ex_et[i + 2][j] + ex_et[i + 2][j + 1] + ex_et[i + 2][j + 2];
0153         stripPhi[index_i][index_j][0] = ex_et[i][j] + ex_et[i + 1][j] + ex_et[i + 2][j];
0154         stripPhi[index_i][index_j][1] = ex_et[i][j + 1] + ex_et[i + 1][j + 1] + ex_et[i + 2][j + 1];
0155         stripPhi[index_i][index_j][2] = ex_et[i][j + 2] + ex_et[i + 1][j + 2] + ex_et[i + 2][j + 2];
0156         for (int k = 0; k < 3; k++) {
0157           et_sumEta[index_i][index_j][k] = ex_et[i + k][j] + ex_et[i + k][j + 1] + ex_et[i + k][j + 2];
0158         }
0159         index_j++;
0160       }
0161       index_i++;
0162     }
0163     for (int i = 0; i < nSTEta; i++) {
0164       for (int j = 0; j < nSTPhi; j++) {
0165         GCTsupertower_t temp;
0166         float supertowerEt = et_sumEta[i][j][0] + et_sumEta[i][j][1] + et_sumEta[i][j][2];
0167         temp.et = supertowerEt;
0168         temp.eta = i;
0169         temp.phi = j;
0170         int peakEta = getPeakBinOf3(stripEta[i][j][0], stripEta[i][j][1], stripEta[i][j][2]);
0171         temp.towerEta = peakEta;
0172         int peakPhi = getPeakBinOf3(stripPhi[i][j][0], stripPhi[i][j][1], stripPhi[i][j][2]);
0173         temp.towerPhi = peakPhi;
0174         float peakEt = ex_et[i * 3 + peakEta][j * 3 + peakPhi];
0175         temp.towerEt = peakEt;
0176         int cEta = getEtCenterOf3(stripEta[i][j][0], stripEta[i][j][1], stripEta[i][j][2]);
0177         temp.centerEta = cEta;
0178         int cPhi = getEtCenterOf3(stripPhi[i][j][0], stripPhi[i][j][1], stripPhi[i][j][2]);
0179         temp.centerPhi = cPhi;
0180         supertower_return[i][j] = temp;
0181       }
0182     }
0183   }
0184 
0185   inline void makeST_hgcal(const float hgcalTowers[nHgcalEta / 2][nHgcalPhi],
0186                            GCTsupertower_t supertower_return[nSTEta][nSTPhi]) {
0187     float et_sumEta[nSTEta][nSTPhi][3];
0188     float stripEta[nSTEta][nSTPhi][3];
0189     float stripPhi[nSTEta][nSTPhi][3];
0190 
0191     int index_i = 0;
0192     int index_j = 0;
0193     for (int i = 0; i < nHgcalEta / 2; i += 3) {  // 18 eta to 6 super towers
0194       index_j = 0;
0195       for (int j = 0; j < nHgcalPhi; j += 3) {  // 72 phi to 24 super towers
0196         stripEta[index_i][index_j][0] = hgcalTowers[i][j] + hgcalTowers[i][j + 1] + hgcalTowers[i][j + 2];
0197         stripEta[index_i][index_j][1] = hgcalTowers[i + 1][j] + hgcalTowers[i + 1][j + 1] + hgcalTowers[i + 1][j + 2];
0198         stripEta[index_i][index_j][2] = hgcalTowers[i + 2][j] + hgcalTowers[i + 2][j + 1] + hgcalTowers[i + 2][j + 2];
0199         stripPhi[index_i][index_j][0] = hgcalTowers[i][j] + hgcalTowers[i + 1][j] + hgcalTowers[i + 2][j];
0200         stripPhi[index_i][index_j][1] = hgcalTowers[i][j + 1] + hgcalTowers[i + 1][j + 1] + hgcalTowers[i + 2][j + 1];
0201         stripPhi[index_i][index_j][2] = hgcalTowers[i][j + 2] + hgcalTowers[i + 1][j + 2] + hgcalTowers[i + 2][j + 2];
0202         for (int k = 0; k < 3; k++) {
0203           et_sumEta[index_i][index_j][k] =
0204               hgcalTowers[i + k][j] + hgcalTowers[i + k][j + 1] + hgcalTowers[i + k][j + 2];
0205         }
0206         index_j++;
0207       }
0208       index_i++;
0209     }
0210 
0211     for (int i = 0; i < nSTEta; i++) {
0212       for (int j = 0; j < nSTPhi; j++) {
0213         GCTsupertower_t temp;
0214         temp.et = 0;
0215         temp.eta = 0;
0216         temp.phi = 0;
0217         temp.towerEta = 0;
0218         temp.towerPhi = 0;
0219         temp.towerEt = 0;
0220         temp.centerEta = 0;
0221         temp.centerPhi = 0;
0222         float supertowerEt = et_sumEta[i][j][0] + et_sumEta[i][j][1] + et_sumEta[i][j][2];
0223         temp.et = supertowerEt;
0224         temp.eta = i;
0225         temp.phi = j;
0226         int peakEta = getPeakBinOf3(stripEta[i][j][0], stripEta[i][j][1], stripEta[i][j][2]);
0227         temp.towerEta = peakEta;
0228         int peakPhi = getPeakBinOf3(stripPhi[i][j][0], stripPhi[i][j][1], stripPhi[i][j][2]);
0229         temp.towerPhi = peakPhi;
0230         float peakEt = hgcalTowers[i * 3 + peakEta][j * 3 + peakPhi];
0231         temp.towerEt = peakEt;
0232         int cEta = getEtCenterOf3(stripEta[i][j][0], stripEta[i][j][1], stripEta[i][j][2]);
0233         temp.centerEta = cEta;
0234         int cPhi = getEtCenterOf3(stripPhi[i][j][0], stripPhi[i][j][1], stripPhi[i][j][2]);
0235         temp.centerPhi = cPhi;
0236         supertower_return[i][j] = temp;
0237       }
0238     }
0239   }
0240 
0241   inline void makeST_hf(const float hfTowers[nHfEta / 2][nHfPhi], GCTsupertower_t supertower_return[nSTEta][nSTPhi]) {
0242     float et_sumEta[nSTEta][nSTPhi][3];
0243     float stripEta[nSTEta][nSTPhi][3];
0244     float stripPhi[nSTEta][nSTPhi][3];
0245 
0246     int index_i = 0;  // 5th and 6th ST to be set 0
0247     int index_j = 0;
0248     for (int i = 0; i < nHfEta / 2; i += 3) {  // 12 eta to 4 super towers
0249       index_j = 0;
0250       for (int j = 0; j < nHfPhi; j += 3) {  // 72 phi to 24 super towers
0251         stripEta[index_i][index_j][0] = hfTowers[i][j] + hfTowers[i][j + 1] + hfTowers[i][j + 2];
0252         stripEta[index_i][index_j][1] = hfTowers[i + 1][j] + hfTowers[i + 1][j + 1] + hfTowers[i + 1][j + 2];
0253         stripEta[index_i][index_j][2] = hfTowers[i + 2][j] + hfTowers[i + 2][j + 1] + hfTowers[i + 2][j + 2];
0254         stripPhi[index_i][index_j][0] = hfTowers[i][j] + hfTowers[i + 1][j] + hfTowers[i + 2][j];
0255         stripPhi[index_i][index_j][1] = hfTowers[i][j + 1] + hfTowers[i + 1][j + 1] + hfTowers[i + 2][j + 1];
0256         stripPhi[index_i][index_j][2] = hfTowers[i][j + 2] + hfTowers[i + 1][j + 2] + hfTowers[i + 2][j + 2];
0257         for (int k = 0; k < 3; k++) {
0258           et_sumEta[index_i][index_j][k] = hfTowers[i + k][j] + hfTowers[i + k][j + 1] + hfTowers[i + k][j + 2];
0259         }
0260         index_j++;
0261       }
0262       index_i++;
0263     }
0264 
0265     for (int i = 0; i < nSTEta; i++) {
0266       for (int j = 0; j < nSTPhi; j++) {
0267         GCTsupertower_t temp;
0268         temp.et = 0;
0269         temp.eta = 0;
0270         temp.phi = 0;
0271         temp.towerEta = 0;
0272         temp.towerPhi = 0;
0273         temp.towerEt = 0;
0274         temp.centerEta = 0;
0275         temp.centerPhi = 0;
0276         if (i < 4) {
0277           float supertowerEt = et_sumEta[i][j][0] + et_sumEta[i][j][1] + et_sumEta[i][j][2];
0278           temp.et = supertowerEt;
0279           temp.eta = i;
0280           temp.phi = j;
0281           int peakEta = getPeakBinOf3(stripEta[i][j][0], stripEta[i][j][1], stripEta[i][j][2]);
0282           temp.towerEta = peakEta;
0283           int peakPhi = getPeakBinOf3(stripPhi[i][j][0], stripPhi[i][j][1], stripPhi[i][j][2]);
0284           temp.towerPhi = peakPhi;
0285           float peakEt = hfTowers[i * 3 + peakEta][j * 3 + peakPhi];
0286           temp.towerEt = peakEt;
0287           int cEta = getEtCenterOf3(stripEta[i][j][0], stripEta[i][j][1], stripEta[i][j][2]);
0288           temp.centerEta = cEta;
0289           int cPhi = getEtCenterOf3(stripPhi[i][j][0], stripPhi[i][j][1], stripPhi[i][j][2]);
0290           temp.centerPhi = cPhi;
0291         }
0292         supertower_return[i][j] = temp;
0293       }
0294     }
0295   }
0296 
0297   inline GCTsupertower_t bestOf2(const GCTsupertower_t& calotp0, const GCTsupertower_t& calotp1) {
0298     GCTsupertower_t x;
0299     x = (calotp0.et > calotp1.et) ? calotp0 : calotp1;
0300     return x;
0301   }
0302 
0303   inline GCTsupertower_t getPeakBin24N(const etaStrip_t& etaStrip) {
0304     GCTsupertower_t best01 = bestOf2(etaStrip.cr[0], etaStrip.cr[1]);
0305     GCTsupertower_t best23 = bestOf2(etaStrip.cr[2], etaStrip.cr[3]);
0306     GCTsupertower_t best45 = bestOf2(etaStrip.cr[4], etaStrip.cr[5]);
0307     GCTsupertower_t best67 = bestOf2(etaStrip.cr[6], etaStrip.cr[7]);
0308     GCTsupertower_t best89 = bestOf2(etaStrip.cr[8], etaStrip.cr[9]);
0309     GCTsupertower_t best1011 = bestOf2(etaStrip.cr[10], etaStrip.cr[11]);
0310     GCTsupertower_t best1213 = bestOf2(etaStrip.cr[12], etaStrip.cr[13]);
0311     GCTsupertower_t best1415 = bestOf2(etaStrip.cr[14], etaStrip.cr[15]);
0312     GCTsupertower_t best1617 = bestOf2(etaStrip.cr[16], etaStrip.cr[17]);
0313     GCTsupertower_t best1819 = bestOf2(etaStrip.cr[18], etaStrip.cr[19]);
0314     GCTsupertower_t best2021 = bestOf2(etaStrip.cr[20], etaStrip.cr[21]);
0315     GCTsupertower_t best2223 = bestOf2(etaStrip.cr[22], etaStrip.cr[23]);
0316 
0317     GCTsupertower_t best0123 = bestOf2(best01, best23);
0318     GCTsupertower_t best4567 = bestOf2(best45, best67);
0319     GCTsupertower_t best891011 = bestOf2(best89, best1011);
0320     GCTsupertower_t best12131415 = bestOf2(best1213, best1415);
0321     GCTsupertower_t best16171819 = bestOf2(best1617, best1819);
0322     GCTsupertower_t best20212223 = bestOf2(best2021, best2223);
0323 
0324     GCTsupertower_t best01234567 = bestOf2(best0123, best4567);
0325     GCTsupertower_t best89101112131415 = bestOf2(best891011, best12131415);
0326     GCTsupertower_t best16to23 = bestOf2(best16171819, best20212223);
0327 
0328     GCTsupertower_t best0to15 = bestOf2(best01234567, best89101112131415);
0329     GCTsupertower_t bestOf24 = bestOf2(best0to15, best16to23);
0330 
0331     return bestOf24;
0332   }
0333 
0334   inline towerMax getPeakBin6N(const etaStripPeak_t& etaStrip) {
0335     towerMax x;
0336 
0337     GCTsupertower_t best01 = bestOf2(etaStrip.pk[0], etaStrip.pk[1]);
0338     GCTsupertower_t best23 = bestOf2(etaStrip.pk[2], etaStrip.pk[3]);
0339     GCTsupertower_t best45 = bestOf2(etaStrip.pk[4], etaStrip.pk[5]);
0340 
0341     GCTsupertower_t best0123 = bestOf2(best01, best23);
0342 
0343     GCTsupertower_t bestOf6 = bestOf2(best0123, best45);
0344 
0345     x.energy = bestOf6.et;
0346     x.phi = bestOf6.phi;
0347     x.eta = bestOf6.eta;
0348     x.energyMax = bestOf6.towerEt;
0349     x.etaMax = bestOf6.towerEta;
0350     x.phiMax = bestOf6.towerPhi;
0351     x.etaCenter = bestOf6.centerEta;
0352     x.phiCenter = bestOf6.centerPhi;
0353     return x;
0354   }
0355 
0356   inline jetInfo getJetPosition(GCTsupertower_t temp[nSTEta][nSTPhi]) {
0357     etaStripPeak_t etaStripPeak;
0358     jetInfo jet;
0359 
0360     for (int i = 0; i < nSTEta; i++) {
0361       etaStrip_t test;
0362       for (int j = 0; j < nSTPhi; j++) {
0363         test.cr[j] = temp[i][j];
0364       }
0365       etaStripPeak.pk[i] = getPeakBin24N(test);
0366     }
0367 
0368     towerMax peakIn6;
0369     peakIn6 = getPeakBin6N(etaStripPeak);
0370 
0371     jet.seedEnergy = peakIn6.energy;
0372     jet.energy = 0;
0373     jet.tauEt = 0;
0374     jet.eta = peakIn6.eta;
0375     jet.phi = peakIn6.phi;
0376     jet.energyMax = peakIn6.energyMax;
0377     jet.etaMax = peakIn6.etaMax;        // overwritten in getJetValues
0378     jet.phiMax = peakIn6.phiMax;        // overwritten in getJetValues
0379     jet.etaCenter = peakIn6.etaCenter;  // overwritten in getJetValues
0380     jet.phiCenter = peakIn6.phiCenter;  // overwritten in getJetValues
0381 
0382     return jet;
0383   }
0384 
0385   inline jetInfo getJetValues(GCTsupertower_t tempX[nSTEta][nSTPhi], int seed_eta, int seed_phi) {
0386     float temp[nSTEta + 2][nSTPhi + 2];
0387     float eta_slice[3] = {0.f, 0.f, 0.f};
0388     jetInfo jet_tmp;
0389 
0390     for (int i = 0; i < nSTEta + 2; i++) {
0391       for (int k = 0; k < nSTPhi + 2; k++) {
0392         temp[i][k] = 0;
0393       }
0394     }
0395 
0396     for (int i = 0; i < nSTEta; i++) {
0397       for (int k = 0; k < nSTPhi; k++) {
0398         temp[i + 1][k + 1] = tempX[i][k].et;
0399       }
0400     }
0401 
0402     int seed_eta1, seed_phi1;
0403 
0404     seed_eta1 = seed_eta;  //to start from corner
0405     seed_phi1 = seed_phi;  //to start from corner
0406     float tmp1, tmp2, tmp3;
0407 
0408     for (int j = 0; j < nSTEta; j++) {
0409       for (int k = 0; k < nSTPhi; k++) {
0410         if (j == seed_eta1 && k == seed_phi1) {
0411           for (int m = 0; m < 3; m++) {
0412             tmp1 = temp[j + m][k];
0413             tmp2 = temp[j + m][k + 1];
0414             tmp3 = temp[j + m][k + 2];
0415             eta_slice[m] = tmp1 + tmp2 + tmp3;
0416           }
0417         }
0418       }
0419     }
0420 
0421     jet_tmp.energy = eta_slice[0] + eta_slice[1] + eta_slice[2];
0422     jet_tmp.tauEt = eta_slice[1];  //set tau Pt to be sum of ST energies in center eta slice */
0423     // To find the jet centre: note that seed supertower is always (1, 1)
0424     jet_tmp.etaCenter =
0425         3 * seed_eta + tempX[seed_eta][seed_phi].centerEta;  //this is the ET weighted eta centre of the ST
0426     jet_tmp.phiCenter =
0427         3 * seed_phi + tempX[seed_eta][seed_phi].centerPhi;  //this is the ET weighted phi centre of the ST
0428     jet_tmp.etaMax = 3 * seed_eta + tempX[seed_eta][seed_phi].towerEta;
0429     jet_tmp.phiMax = 3 * seed_phi + tempX[seed_eta][seed_phi].towerPhi;
0430 
0431     // set the used supertower ets to 0
0432     for (int i = 0; i < nSTEta; i++) {
0433       if (i + 1 >= seed_eta && i <= seed_eta + 1) {
0434         for (int k = 0; k < nSTPhi; k++) {
0435           if (k + 1 >= seed_phi && k <= seed_phi + 1)
0436             tempX[i][k].et = 0;
0437         }
0438       }
0439     }
0440 
0441     return jet_tmp;
0442   }
0443 
0444   inline jetInfo getRegion(GCTsupertower_t temp[nSTEta][nSTPhi]) {
0445     jetInfo jet_tmp, jet;
0446     jet_tmp = getJetPosition(temp);
0447     int seed_phi = jet_tmp.phi;
0448     int seed_eta = jet_tmp.eta;
0449     float seed_energy = jet_tmp.seedEnergy;
0450     jet = getJetValues(temp, seed_eta, seed_phi);
0451     if (seed_energy > 10.) {  // suppress <= 10 GeV ST as seed
0452       jet_tmp.energy = jet.energy;
0453       jet_tmp.tauEt = jet.tauEt;
0454     } else {
0455       jet_tmp.energy = 0.;
0456       jet_tmp.tauEt = 0.;
0457     }
0458     jet_tmp.etaCenter = jet.etaCenter;  // this is the ET weighted eta centre of the ST
0459     jet_tmp.phiCenter = jet.phiCenter;  // this is the ET weighted phi centre of the ST
0460     jet_tmp.etaMax = jet.etaMax;        // this is the leading tower eta in the ST
0461     jet_tmp.phiMax = jet.phiMax;        // this is the leading tower phi in the ST
0462     return jet_tmp;
0463   }
0464 
0465   inline bool compareByEt(l1tp2::Phase2L1CaloJet i, l1tp2::Phase2L1CaloJet j) { return (i.jetEt() > j.jetEt()); };
0466 
0467 }  // namespace gctobj
0468 
0469 #endif