Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-05 05:04:36

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