File indexing completed on 2024-05-10 02:20:59
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) {
0148 index_j = 0;
0149 for (int j = 0; j < nBarrelPhi; j += 3) {
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) {
0194 index_j = 0;
0195 for (int j = 0; j < nHgcalPhi; j += 3) {
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;
0247 int index_j = 0;
0248 for (int i = 0; i < nHfEta / 2; i += 3) {
0249 index_j = 0;
0250 for (int j = 0; j < nHfPhi; j += 3) {
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;
0378 jet.phiMax = peakIn6.phiMax;
0379 jet.etaCenter = peakIn6.etaCenter;
0380 jet.phiCenter = peakIn6.phiCenter;
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;
0405 seed_phi1 = seed_phi;
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];
0423
0424 jet_tmp.etaCenter =
0425 3 * seed_eta + tempX[seed_eta][seed_phi].centerEta;
0426 jet_tmp.phiCenter =
0427 3 * seed_phi + tempX[seed_eta][seed_phi].centerPhi;
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
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], float TTseedThreshold) {
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 float seed_tower_energy = jet_tmp.energyMax;
0451 jet = getJetValues(temp, seed_eta, seed_phi);
0452 if (seed_energy > 10. &&
0453 seed_tower_energy >
0454 TTseedThreshold) {
0455 jet_tmp.energy = jet.energy;
0456 jet_tmp.tauEt = jet.tauEt;
0457 } else {
0458 jet_tmp.energy = 0.;
0459 jet_tmp.tauEt = 0.;
0460 }
0461 jet_tmp.etaCenter = jet.etaCenter;
0462 jet_tmp.phiCenter = jet.phiCenter;
0463 jet_tmp.etaMax = jet.etaMax;
0464 jet_tmp.phiMax = jet.phiMax;
0465 return jet_tmp;
0466 }
0467
0468 inline bool compareByEt(l1tp2::Phase2L1CaloJet i, l1tp2::Phase2L1CaloJet j) { return (i.jetEt() > j.jetEt()); };
0469
0470 }
0471
0472 #endif