CMS 3D CMS Logo

Phase2L1CaloJetEmulator.h
Go to the documentation of this file.
1 #ifndef L1TRIGGER_L1CALOTRIGGER_PHASE2L1CALOJETEMULATOR_H
2 #define L1TRIGGER_L1CALOTRIGGER_PHASE2L1CALOJETEMULATOR_H
3 
5 #include <cstdlib>
6 
7 static constexpr int nBarrelEta = 34;
8 static constexpr int nBarrelPhi = 72;
9 static constexpr int nHgcalEta = 36;
10 static constexpr int nHgcalPhi = 72;
11 static constexpr int nHfEta = 24;
12 static constexpr int nHfPhi = 72;
13 static constexpr int nSTEta = 6;
14 static constexpr int nSTEta_hf = 4;
15 static constexpr int nSTPhi = 24;
16 static constexpr int nJets = 10;
17 
18 namespace gctobj {
19 
20  class towerMax {
21  public:
22  float energy;
23  int phi;
24  int eta;
25  float energyMax;
26  int phiMax;
27  int etaMax;
28  int phiCenter;
29  int etaCenter;
30 
32  energy = 0;
33  phi = 0;
34  eta = 0;
35  energyMax = 0;
36  phiMax = 0;
37  etaMax = 0;
38  phiCenter = 0;
39  etaCenter = 0;
40  }
41  };
42 
43  class jetInfo {
44  public:
45  float seedEnergy;
46  float energy;
47  float tauEt;
48  int phi;
49  int eta;
50  float energyMax;
51  int phiMax;
52  int etaMax;
53  int phiCenter;
54  int etaCenter;
55 
56  jetInfo() {
57  seedEnergy = 0;
58  energy = 0;
59  tauEt = 0;
60  phi = 0;
61  eta = 0;
62  energyMax = 0;
63  phiMax = 0;
64  etaMax = 0;
65  phiCenter = 0;
66  etaCenter = 0;
67  }
68  };
69 
70  typedef struct {
71  float et;
72  int eta;
73  int phi;
74  float towerEt;
75  int towerEta;
76  int towerPhi;
77  int centerEta;
78  int centerPhi;
80 
81  typedef struct {
83  } etaStrip_t;
84 
85  typedef struct {
87  } hgcalRegion_t;
88 
89  typedef struct {
92 
93  typedef struct {
94  float et;
95  float hoe;
96  } GCTtower_t;
97 
98  typedef struct {
100  } GCTintTowers_t;
101 
102  inline int getPeakBinOf3(float et0, float et1, float et2) {
103  int x;
104  float temp;
105  if (et0 > et1) {
106  x = 0;
107  temp = et0;
108  } else {
109  x = 1;
110  temp = et1;
111  }
112  if (et2 > temp) {
113  x = 2;
114  }
115  return x;
116  }
117 
118  inline int getEtCenterOf3(float et0, float et1, float et2) {
119  float etSum = et0 + et1 + et2;
120  float iEtSum = 0.5 * et0 + 1.5 * et1 + 2.5 * et2;
121  int iAve = 0xEEF;
122  if (iEtSum <= etSum)
123  iAve = 0;
124  else if (iEtSum <= 2 * etSum)
125  iAve = 1;
126  else
127  iAve = 2;
128  return iAve;
129  }
130 
131  inline void makeST(const float GCTintTowers[nBarrelEta / 2][nBarrelPhi],
132  GCTsupertower_t supertower_return[nSTEta][nSTPhi]) {
133  float et_sumEta[nSTEta][nSTPhi][3];
134  float stripEta[nSTEta][nSTPhi][3];
135  float stripPhi[nSTEta][nSTPhi][3];
136 
137  float ex_et[nBarrelEta / 2 + 1][nBarrelPhi];
138  for (int j = 0; j < nBarrelPhi; j++) {
139  ex_et[nBarrelEta / 2][j] = 0;
140  for (int i = 0; i < nBarrelEta / 2; i++) {
141  ex_et[i][j] = GCTintTowers[i][j];
142  }
143  }
144 
145  int index_i = 0;
146  int index_j = 0;
147  for (int i = 0; i < nBarrelEta / 2 + 1; i += 3) { // 17+1 to divide into 6 super towers
148  index_j = 0;
149  for (int j = 0; j < nBarrelPhi; j += 3) { // 72 phi to 24 super towers
150  stripEta[index_i][index_j][0] = ex_et[i][j] + ex_et[i][j + 1] + ex_et[i][j + 2];
151  stripEta[index_i][index_j][1] = ex_et[i + 1][j] + ex_et[i + 1][j + 1] + ex_et[i + 1][j + 2];
152  stripEta[index_i][index_j][2] = ex_et[i + 2][j] + ex_et[i + 2][j + 1] + ex_et[i + 2][j + 2];
153  stripPhi[index_i][index_j][0] = ex_et[i][j] + ex_et[i + 1][j] + ex_et[i + 2][j];
154  stripPhi[index_i][index_j][1] = ex_et[i][j + 1] + ex_et[i + 1][j + 1] + ex_et[i + 2][j + 1];
155  stripPhi[index_i][index_j][2] = ex_et[i][j + 2] + ex_et[i + 1][j + 2] + ex_et[i + 2][j + 2];
156  for (int k = 0; k < 3; k++) {
157  et_sumEta[index_i][index_j][k] = ex_et[i + k][j] + ex_et[i + k][j + 1] + ex_et[i + k][j + 2];
158  }
159  index_j++;
160  }
161  index_i++;
162  }
163  for (int i = 0; i < nSTEta; i++) {
164  for (int j = 0; j < nSTPhi; j++) {
166  float supertowerEt = et_sumEta[i][j][0] + et_sumEta[i][j][1] + et_sumEta[i][j][2];
167  temp.et = supertowerEt;
168  temp.eta = i;
169  temp.phi = j;
170  int peakEta = getPeakBinOf3(stripEta[i][j][0], stripEta[i][j][1], stripEta[i][j][2]);
171  temp.towerEta = peakEta;
172  int peakPhi = getPeakBinOf3(stripPhi[i][j][0], stripPhi[i][j][1], stripPhi[i][j][2]);
173  temp.towerPhi = peakPhi;
174  float peakEt = ex_et[i * 3 + peakEta][j * 3 + peakPhi];
175  temp.towerEt = peakEt;
176  int cEta = getEtCenterOf3(stripEta[i][j][0], stripEta[i][j][1], stripEta[i][j][2]);
177  temp.centerEta = cEta;
178  int cPhi = getEtCenterOf3(stripPhi[i][j][0], stripPhi[i][j][1], stripPhi[i][j][2]);
179  temp.centerPhi = cPhi;
180  supertower_return[i][j] = temp;
181  }
182  }
183  }
184 
185  inline void makeST_hgcal(const float hgcalTowers[nHgcalEta / 2][nHgcalPhi],
186  GCTsupertower_t supertower_return[nSTEta][nSTPhi]) {
187  float et_sumEta[nSTEta][nSTPhi][3];
188  float stripEta[nSTEta][nSTPhi][3];
189  float stripPhi[nSTEta][nSTPhi][3];
190 
191  int index_i = 0;
192  int index_j = 0;
193  for (int i = 0; i < nHgcalEta / 2; i += 3) { // 18 eta to 6 super towers
194  index_j = 0;
195  for (int j = 0; j < nHgcalPhi; j += 3) { // 72 phi to 24 super towers
196  stripEta[index_i][index_j][0] = hgcalTowers[i][j] + hgcalTowers[i][j + 1] + hgcalTowers[i][j + 2];
197  stripEta[index_i][index_j][1] = hgcalTowers[i + 1][j] + hgcalTowers[i + 1][j + 1] + hgcalTowers[i + 1][j + 2];
198  stripEta[index_i][index_j][2] = hgcalTowers[i + 2][j] + hgcalTowers[i + 2][j + 1] + hgcalTowers[i + 2][j + 2];
199  stripPhi[index_i][index_j][0] = hgcalTowers[i][j] + hgcalTowers[i + 1][j] + hgcalTowers[i + 2][j];
200  stripPhi[index_i][index_j][1] = hgcalTowers[i][j + 1] + hgcalTowers[i + 1][j + 1] + hgcalTowers[i + 2][j + 1];
201  stripPhi[index_i][index_j][2] = hgcalTowers[i][j + 2] + hgcalTowers[i + 1][j + 2] + hgcalTowers[i + 2][j + 2];
202  for (int k = 0; k < 3; k++) {
203  et_sumEta[index_i][index_j][k] =
204  hgcalTowers[i + k][j] + hgcalTowers[i + k][j + 1] + hgcalTowers[i + k][j + 2];
205  }
206  index_j++;
207  }
208  index_i++;
209  }
210 
211  for (int i = 0; i < nSTEta; i++) {
212  for (int j = 0; j < nSTPhi; j++) {
214  temp.et = 0;
215  temp.eta = 0;
216  temp.phi = 0;
217  temp.towerEta = 0;
218  temp.towerPhi = 0;
219  temp.towerEt = 0;
220  temp.centerEta = 0;
221  temp.centerPhi = 0;
222  float supertowerEt = et_sumEta[i][j][0] + et_sumEta[i][j][1] + et_sumEta[i][j][2];
223  temp.et = supertowerEt;
224  temp.eta = i;
225  temp.phi = j;
226  int peakEta = getPeakBinOf3(stripEta[i][j][0], stripEta[i][j][1], stripEta[i][j][2]);
227  temp.towerEta = peakEta;
228  int peakPhi = getPeakBinOf3(stripPhi[i][j][0], stripPhi[i][j][1], stripPhi[i][j][2]);
229  temp.towerPhi = peakPhi;
230  float peakEt = hgcalTowers[i * 3 + peakEta][j * 3 + peakPhi];
231  temp.towerEt = peakEt;
232  int cEta = getEtCenterOf3(stripEta[i][j][0], stripEta[i][j][1], stripEta[i][j][2]);
233  temp.centerEta = cEta;
234  int cPhi = getEtCenterOf3(stripPhi[i][j][0], stripPhi[i][j][1], stripPhi[i][j][2]);
235  temp.centerPhi = cPhi;
236  supertower_return[i][j] = temp;
237  }
238  }
239  }
240 
241  inline void makeST_hf(const float hfTowers[nHfEta / 2][nHfPhi], GCTsupertower_t supertower_return[nSTEta][nSTPhi]) {
242  float et_sumEta[nSTEta][nSTPhi][3];
243  float stripEta[nSTEta][nSTPhi][3];
244  float stripPhi[nSTEta][nSTPhi][3];
245 
246  int index_i = 0; // 5th and 6th ST to be set 0
247  int index_j = 0;
248  for (int i = 0; i < nHfEta / 2; i += 3) { // 12 eta to 4 super towers
249  index_j = 0;
250  for (int j = 0; j < nHfPhi; j += 3) { // 72 phi to 24 super towers
251  stripEta[index_i][index_j][0] = hfTowers[i][j] + hfTowers[i][j + 1] + hfTowers[i][j + 2];
252  stripEta[index_i][index_j][1] = hfTowers[i + 1][j] + hfTowers[i + 1][j + 1] + hfTowers[i + 1][j + 2];
253  stripEta[index_i][index_j][2] = hfTowers[i + 2][j] + hfTowers[i + 2][j + 1] + hfTowers[i + 2][j + 2];
254  stripPhi[index_i][index_j][0] = hfTowers[i][j] + hfTowers[i + 1][j] + hfTowers[i + 2][j];
255  stripPhi[index_i][index_j][1] = hfTowers[i][j + 1] + hfTowers[i + 1][j + 1] + hfTowers[i + 2][j + 1];
256  stripPhi[index_i][index_j][2] = hfTowers[i][j + 2] + hfTowers[i + 1][j + 2] + hfTowers[i + 2][j + 2];
257  for (int k = 0; k < 3; k++) {
258  et_sumEta[index_i][index_j][k] = hfTowers[i + k][j] + hfTowers[i + k][j + 1] + hfTowers[i + k][j + 2];
259  }
260  index_j++;
261  }
262  index_i++;
263  }
264 
265  for (int i = 0; i < nSTEta; i++) {
266  for (int j = 0; j < nSTPhi; j++) {
268  temp.et = 0;
269  temp.eta = 0;
270  temp.phi = 0;
271  temp.towerEta = 0;
272  temp.towerPhi = 0;
273  temp.towerEt = 0;
274  temp.centerEta = 0;
275  temp.centerPhi = 0;
276  if (i < 4) {
277  float supertowerEt = et_sumEta[i][j][0] + et_sumEta[i][j][1] + et_sumEta[i][j][2];
278  temp.et = supertowerEt;
279  temp.eta = i;
280  temp.phi = j;
281  int peakEta = getPeakBinOf3(stripEta[i][j][0], stripEta[i][j][1], stripEta[i][j][2]);
282  temp.towerEta = peakEta;
283  int peakPhi = getPeakBinOf3(stripPhi[i][j][0], stripPhi[i][j][1], stripPhi[i][j][2]);
284  temp.towerPhi = peakPhi;
285  float peakEt = hfTowers[i * 3 + peakEta][j * 3 + peakPhi];
286  temp.towerEt = peakEt;
287  int cEta = getEtCenterOf3(stripEta[i][j][0], stripEta[i][j][1], stripEta[i][j][2]);
288  temp.centerEta = cEta;
289  int cPhi = getEtCenterOf3(stripPhi[i][j][0], stripPhi[i][j][1], stripPhi[i][j][2]);
290  temp.centerPhi = cPhi;
291  }
292  supertower_return[i][j] = temp;
293  }
294  }
295  }
296 
297  inline GCTsupertower_t bestOf2(const GCTsupertower_t& calotp0, const GCTsupertower_t& calotp1) {
299  x = (calotp0.et > calotp1.et) ? calotp0 : calotp1;
300  return x;
301  }
302 
303  inline GCTsupertower_t getPeakBin24N(const etaStrip_t& etaStrip) {
304  GCTsupertower_t best01 = bestOf2(etaStrip.cr[0], etaStrip.cr[1]);
305  GCTsupertower_t best23 = bestOf2(etaStrip.cr[2], etaStrip.cr[3]);
306  GCTsupertower_t best45 = bestOf2(etaStrip.cr[4], etaStrip.cr[5]);
307  GCTsupertower_t best67 = bestOf2(etaStrip.cr[6], etaStrip.cr[7]);
308  GCTsupertower_t best89 = bestOf2(etaStrip.cr[8], etaStrip.cr[9]);
309  GCTsupertower_t best1011 = bestOf2(etaStrip.cr[10], etaStrip.cr[11]);
310  GCTsupertower_t best1213 = bestOf2(etaStrip.cr[12], etaStrip.cr[13]);
311  GCTsupertower_t best1415 = bestOf2(etaStrip.cr[14], etaStrip.cr[15]);
312  GCTsupertower_t best1617 = bestOf2(etaStrip.cr[16], etaStrip.cr[17]);
313  GCTsupertower_t best1819 = bestOf2(etaStrip.cr[18], etaStrip.cr[19]);
314  GCTsupertower_t best2021 = bestOf2(etaStrip.cr[20], etaStrip.cr[21]);
315  GCTsupertower_t best2223 = bestOf2(etaStrip.cr[22], etaStrip.cr[23]);
316 
317  GCTsupertower_t best0123 = bestOf2(best01, best23);
318  GCTsupertower_t best4567 = bestOf2(best45, best67);
319  GCTsupertower_t best891011 = bestOf2(best89, best1011);
320  GCTsupertower_t best12131415 = bestOf2(best1213, best1415);
321  GCTsupertower_t best16171819 = bestOf2(best1617, best1819);
322  GCTsupertower_t best20212223 = bestOf2(best2021, best2223);
323 
324  GCTsupertower_t best01234567 = bestOf2(best0123, best4567);
325  GCTsupertower_t best89101112131415 = bestOf2(best891011, best12131415);
326  GCTsupertower_t best16to23 = bestOf2(best16171819, best20212223);
327 
328  GCTsupertower_t best0to15 = bestOf2(best01234567, best89101112131415);
329  GCTsupertower_t bestOf24 = bestOf2(best0to15, best16to23);
330 
331  return bestOf24;
332  }
333 
334  inline towerMax getPeakBin6N(const etaStripPeak_t& etaStrip) {
335  towerMax x;
336 
337  GCTsupertower_t best01 = bestOf2(etaStrip.pk[0], etaStrip.pk[1]);
338  GCTsupertower_t best23 = bestOf2(etaStrip.pk[2], etaStrip.pk[3]);
339  GCTsupertower_t best45 = bestOf2(etaStrip.pk[4], etaStrip.pk[5]);
340 
341  GCTsupertower_t best0123 = bestOf2(best01, best23);
342 
343  GCTsupertower_t bestOf6 = bestOf2(best0123, best45);
344 
345  x.energy = bestOf6.et;
346  x.phi = bestOf6.phi;
347  x.eta = bestOf6.eta;
348  x.energyMax = bestOf6.towerEt;
349  x.etaMax = bestOf6.towerEta;
350  x.phiMax = bestOf6.towerPhi;
351  x.etaCenter = bestOf6.centerEta;
352  x.phiCenter = bestOf6.centerPhi;
353  return x;
354  }
355 
357  etaStripPeak_t etaStripPeak;
358  jetInfo jet;
359 
360  for (int i = 0; i < nSTEta; i++) {
362  for (int j = 0; j < nSTPhi; j++) {
363  test.cr[j] = temp[i][j];
364  }
365  etaStripPeak.pk[i] = getPeakBin24N(test);
366  }
367 
368  towerMax peakIn6;
369  peakIn6 = getPeakBin6N(etaStripPeak);
370 
371  jet.seedEnergy = peakIn6.energy;
372  jet.energy = 0;
373  jet.tauEt = 0;
374  jet.eta = peakIn6.eta;
375  jet.phi = peakIn6.phi;
376  jet.energyMax = peakIn6.energyMax;
377  jet.etaMax = peakIn6.etaMax; // overwritten in getJetValues
378  jet.phiMax = peakIn6.phiMax; // overwritten in getJetValues
379  jet.etaCenter = peakIn6.etaCenter; // overwritten in getJetValues
380  jet.phiCenter = peakIn6.phiCenter; // overwritten in getJetValues
381 
382  return jet;
383  }
384 
385  inline jetInfo getJetValues(GCTsupertower_t tempX[nSTEta][nSTPhi], int seed_eta, int seed_phi) {
386  float temp[nSTEta + 2][nSTPhi + 2];
387  float eta_slice[3] = {0.f, 0.f, 0.f};
388  jetInfo jet_tmp;
389 
390  for (int i = 0; i < nSTEta + 2; i++) {
391  for (int k = 0; k < nSTPhi + 2; k++) {
392  temp[i][k] = 0;
393  }
394  }
395 
396  for (int i = 0; i < nSTEta; i++) {
397  for (int k = 0; k < nSTPhi; k++) {
398  temp[i + 1][k + 1] = tempX[i][k].et;
399  }
400  }
401 
402  int seed_eta1, seed_phi1;
403 
404  seed_eta1 = seed_eta; //to start from corner
405  seed_phi1 = seed_phi; //to start from corner
406  float tmp1, tmp2, tmp3;
407 
408  for (int j = 0; j < nSTEta; j++) {
409  for (int k = 0; k < nSTPhi; k++) {
410  if (j == seed_eta1 && k == seed_phi1) {
411  for (int m = 0; m < 3; m++) {
412  tmp1 = temp[j + m][k];
413  tmp2 = temp[j + m][k + 1];
414  tmp3 = temp[j + m][k + 2];
415  eta_slice[m] = tmp1 + tmp2 + tmp3;
416  }
417  }
418  }
419  }
420 
421  jet_tmp.energy = eta_slice[0] + eta_slice[1] + eta_slice[2];
422  jet_tmp.tauEt = eta_slice[1]; //set tau Pt to be sum of ST energies in center eta slice */
423  // To find the jet centre: note that seed supertower is always (1, 1)
424  jet_tmp.etaCenter =
425  3 * seed_eta + tempX[seed_eta][seed_phi].centerEta; //this is the ET weighted eta centre of the ST
426  jet_tmp.phiCenter =
427  3 * seed_phi + tempX[seed_eta][seed_phi].centerPhi; //this is the ET weighted phi centre of the ST
428  jet_tmp.etaMax = 3 * seed_eta + tempX[seed_eta][seed_phi].towerEta;
429  jet_tmp.phiMax = 3 * seed_phi + tempX[seed_eta][seed_phi].towerPhi;
430 
431  // set the used supertower ets to 0
432  for (int i = 0; i < nSTEta; i++) {
433  if (i + 1 >= seed_eta && i <= seed_eta + 1) {
434  for (int k = 0; k < nSTPhi; k++) {
435  if (k + 1 >= seed_phi && k <= seed_phi + 1)
436  tempX[i][k].et = 0;
437  }
438  }
439  }
440 
441  return jet_tmp;
442  }
443 
444  inline jetInfo getRegion(GCTsupertower_t temp[nSTEta][nSTPhi], float TTseedThreshold) {
445  jetInfo jet_tmp, jet;
446  jet_tmp = getJetPosition(temp);
447  int seed_phi = jet_tmp.phi;
448  int seed_eta = jet_tmp.eta;
449  float seed_energy = jet_tmp.seedEnergy;
450  float seed_tower_energy = jet_tmp.energyMax;
451  jet = getJetValues(temp, seed_eta, seed_phi);
452  if (seed_energy > 10. &&
453  seed_tower_energy >
454  TTseedThreshold) { // suppress <= 10 GeV ST as ST seed and <=5 GeV (3 GeV) as max TT in ST barrel/HF (endcap)
455  jet_tmp.energy = jet.energy;
456  jet_tmp.tauEt = jet.tauEt;
457  } else {
458  jet_tmp.energy = 0.;
459  jet_tmp.tauEt = 0.;
460  }
461  jet_tmp.etaCenter = jet.etaCenter; // this is the ET weighted eta centre of the ST
462  jet_tmp.phiCenter = jet.phiCenter; // this is the ET weighted phi centre of the ST
463  jet_tmp.etaMax = jet.etaMax; // this is the leading tower eta in the ST
464  jet_tmp.phiMax = jet.phiMax; // this is the leading tower phi in the ST
465  return jet_tmp;
466  }
467 
468  inline bool compareByEt(l1tp2::Phase2L1CaloJet i, l1tp2::Phase2L1CaloJet j) { return (i.jetEt() > j.jetEt()); };
469 
470 } // namespace gctobj
471 
472 #endif
static constexpr int nSTEta_hf
GCTsupertower_t cr[nSTPhi]
static constexpr int nBarrelPhi
static constexpr int nSTEta
GCTsupertower_t getPeakBin24N(const etaStrip_t &etaStrip)
static constexpr int nHgcalEta
GCTsupertower_t pk[nSTEta]
static constexpr int nJets
int getEtCenterOf3(float et0, float et1, float et2)
GCTsupertower_t bestOf2(const GCTsupertower_t &calotp0, const GCTsupertower_t &calotp1)
static constexpr int nHfPhi
void makeST_hf(const float hfTowers[nHfEta/2][nHfPhi], GCTsupertower_t supertower_return[nSTEta][nSTPhi])
static constexpr int nSTPhi
jetInfo getJetPosition(GCTsupertower_t temp[nSTEta][nSTPhi])
int getPeakBinOf3(float et0, float et1, float et2)
void makeST(const float GCTintTowers[nBarrelEta/2][nBarrelPhi], GCTsupertower_t supertower_return[nSTEta][nSTPhi])
static constexpr int nHgcalPhi
void makeST_hgcal(const float hgcalTowers[nHgcalEta/2][nHgcalPhi], GCTsupertower_t supertower_return[nSTEta][nSTPhi])
static constexpr int nBarrelEta
static constexpr int nHfEta
float x
towerMax getPeakBin6N(const etaStripPeak_t &etaStrip)
jetInfo getJetValues(GCTsupertower_t tempX[nSTEta][nSTPhi], int seed_eta, int seed_phi)
jetInfo getRegion(GCTsupertower_t temp[nSTEta][nSTPhi], float TTseedThreshold)
bool compareByEt(l1tp2::Phase2L1CaloJet i, l1tp2::Phase2L1CaloJet j)