ATLAS Offline Software
Loading...
Searching...
No Matches
gFEXJwoJAlgo.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3*/
4//***************************************************************************
5// gFEXJwoJAlg - Jets without jets algorithm for gFEX
6// -------------------
7// begin : 21 12 2023
8// email : cecilia.tosciri@cern.ch
9//***************************************************************************
10
11
12
13#include "gFEXJwoJAlgo.h"
15
16#include <cmath> //for std::sqrt
17
18namespace LVL1 {
19
20 // default constructor for persistency
21gFEXJwoJAlgo::gFEXJwoJAlgo(const std::string& type, const std::string& name, const IInterface* parent):
22 AthAlgTool(type, name, parent)
23 {
24 declareInterface<IgFEXJwoJAlgo>(this);
25 }
26
27
29
30 ATH_CHECK(m_DBToolKey.initialize());
31
32 return StatusCode::SUCCESS;
33
34}
35
36
37void gFEXJwoJAlgo::setAlgoConstant(int aFPGA_A, int bFPGA_A,
38 int aFPGA_B, int bFPGA_B,
39 int aFPGA_C, int bFPGA_C,
40 int gXE_seedThrA, int gXE_seedThrB, int gXE_seedThrC) {
41 m_aFPGA_A = aFPGA_A;
42 m_bFPGA_A = bFPGA_A;
43 m_aFPGA_B = aFPGA_B;
44 m_bFPGA_B = bFPGA_B;
45 m_aFPGA_C = aFPGA_C;
46 m_bFPGA_C = bFPGA_C;
47 m_gBlockthresholdA = gXE_seedThrA;
48 m_gBlockthresholdB = gXE_seedThrB;
49 m_gBlockthresholdC = gXE_seedThrC;
50}
51
52std::vector<std::unique_ptr<gFEXJwoJTOB>> gFEXJwoJAlgo::jwojAlgo(const gTowersType& Atwr, int pucA_JWJ,
53 const gTowersType& Btwr, int pucB_JWJ,
54 const gTowersType& Ctwr, int pucC_JWJ,
55 std::array<int32_t, 4> & outTOB) const {
56
57
59 if (!myDBTool.isValid()) {
60 ATH_MSG_ERROR("Could not retrieve DB tool " << m_DBToolKey);
61 throw std::runtime_error("Could not retrieve DB tool");
62 }
63
64 std::string fwVersion = myDBTool->get_FWVersion();
65 int major = std::stoi(fwVersion);
66 bool SumETfast = (major >= 1);
67 bool metRho = (major >= 2);
68
69 // find gBlocks
70 gTowersType AgBlk;
71 gTowersType Ascaled;
72
73 gTowersType BgBlk;
74 gTowersType Bscaled;
75
76 gTowersType CgBlk;
77 gTowersType Cscaled;
78
79 gTowersType hasSeed;
80
81 gBlockAB(Atwr, AgBlk, hasSeed, m_gBlockthresholdA);
82 gBlockAB(Btwr, BgBlk, hasSeed, m_gBlockthresholdB);
83 gBlockAB(Ctwr, CgBlk, hasSeed, m_gBlockthresholdC);
84
85
86 // switch to 10 bit number
87 // DMS -- do we eventaully need to check for overflows here?
88
89 for(int irow = 0; irow<FEXAlgoSpaceDefs::ABCrows; irow++){
90 for(int jcolumn = 0; jcolumn<FEXAlgoSpaceDefs::ABcolumns; jcolumn++){
91 Ascaled[irow][jcolumn] = Atwr[irow][jcolumn] >> 2;
92 AgBlk[irow][jcolumn] = AgBlk[irow][jcolumn] >> 2;
93
94 Bscaled[irow][jcolumn] = Btwr[irow][jcolumn] >> 2;
95 BgBlk[irow][jcolumn] = BgBlk[irow][jcolumn] >> 2;
96
97 Cscaled[irow][jcolumn] = Ctwr[irow][jcolumn] >> 2;
98 CgBlk[irow][jcolumn] = CgBlk[irow][jcolumn] >> 2;
99
100 }
101 }
102
103
104 //FPGA A observables
105 int A_MHT_x = 0x0;
106 int A_MHT_y = 0x0;
107 int A_MST_x = 0x0;
108 int A_MST_y = 0x0;
109 int A_MET_x = 0x0;
110 int A_MET_y = 0x0;
111
112 int A_eth = 0x0;
113 int A_ets = 0x0;
114 int A_etw = 0x0;
115
116 //FPGA B observables
117 int B_MHT_x = 0x0;
118 int B_MHT_y = 0x0;
119 int B_MST_x = 0x0;
120 int B_MST_y = 0x0;
121 int B_MET_x = 0x0;
122 int B_MET_y = 0x0;
123
124 int B_eth = 0x0;
125 int B_ets = 0x0;
126 int B_etw = 0x0;
127
128 //FPGA C observables
129 int C_MHT_x = 0x0;
130 int C_MHT_y = 0x0;
131 int C_MST_x = 0x0;
132 int C_MST_y = 0x0;
133 int C_MET_x = 0x0;
134 int C_MET_y = 0x0;
135
136 int C_eth = 0x0;
137 int C_ets = 0x0;
138 int C_etw = 0x0;
139
140 //Global observables
141 int MHT_x = 0x0;
142 int MHT_y = 0x0;
143 int MST_x = 0x0;
144 int MST_y = 0x0;
145 int MET_x = 0x0;
146 int MET_y = 0x0;
147
148 int ETH = 0x0;
149 int ETS = 0x0;
150 int ETW = 0x0;
151
152 int total_sumEt = 0x0;
153 int MET = 0x0;
154
155
156 // will need to hard code etFPGA ,a's and b's
157 int etBprime = 0;
158
159 if (metRho) metFPGA_rho(0, Ascaled, pucA_JWJ, AgBlk, m_gBlockthresholdA, m_aFPGA_A, m_bFPGA_A, A_MHT_x, A_MHT_y, A_MST_x, A_MST_y, A_MET_x, A_MET_y);
160 else metFPGA(0, Ascaled, AgBlk, m_gBlockthresholdA, m_aFPGA_A, m_bFPGA_A, A_MHT_x, A_MHT_y, A_MST_x, A_MST_y, A_MET_x, A_MET_y);
161 if (SumETfast) etFastFPGA(0, Ascaled, AgBlk, m_gBlockthresholdA, m_aFPGA_A, etBprime, A_eth, A_ets, A_etw);
162 else etFPGA(0, Ascaled, AgBlk, m_gBlockthresholdA, m_aFPGA_A, etBprime, A_eth, A_ets, A_etw);
163
164 if (metRho) metFPGA_rho(1, Bscaled, pucB_JWJ, BgBlk, m_gBlockthresholdB, m_aFPGA_B, m_bFPGA_B, B_MHT_x, B_MHT_y, B_MST_x, B_MST_y, B_MET_x, B_MET_y);
165 else metFPGA(1, Bscaled, BgBlk, m_gBlockthresholdB, m_aFPGA_B, m_bFPGA_B, B_MHT_x, B_MHT_y, B_MST_x, B_MST_y, B_MET_x, B_MET_y);
166 if (SumETfast) etFastFPGA(1, Bscaled, BgBlk, m_gBlockthresholdB, m_aFPGA_B, etBprime, B_eth, B_ets, B_etw);
167 else etFPGA(1, Bscaled, BgBlk, m_gBlockthresholdB, m_aFPGA_B, etBprime, B_eth, B_ets, B_etw);
168
169 if (metRho) metFPGA_rho(2, Cscaled, pucC_JWJ, CgBlk, m_gBlockthresholdC, m_aFPGA_C, m_bFPGA_C, C_MHT_x, C_MHT_y, C_MST_x, C_MST_y, C_MET_x, C_MET_y);
170 else metFPGA(2, Cscaled, CgBlk, m_gBlockthresholdC, m_aFPGA_C, m_bFPGA_C, C_MHT_x, C_MHT_y, C_MST_x, C_MST_y, C_MET_x, C_MET_y);
171 if (SumETfast) etFastFPGA(2, Cscaled, CgBlk, m_gBlockthresholdC, m_aFPGA_C, etBprime, C_eth, C_ets, C_etw);
172 else etFPGA(2, Cscaled, CgBlk, m_gBlockthresholdC, m_aFPGA_C, etBprime, C_eth, C_ets, C_etw);
173
174 metTotal(A_MHT_x, A_MHT_y, B_MHT_x, B_MHT_y, C_MHT_x, C_MHT_y, MHT_x, MHT_y);
175 metTotal(A_MST_x, A_MST_y, B_MST_x, B_MST_y, C_MST_x, C_MST_y, MST_x, MST_y);
176 metTotal(A_MET_x, A_MET_y, B_MET_x, B_MET_y, C_MET_x, C_MET_y, MET_x, MET_y);
177
178 etTotal(A_eth, B_eth, C_eth, ETH);
179 etTotal(A_ets, B_ets, C_ets, ETS);
180 etTotal(A_etw, B_etw, C_etw, ETW);
181 total_sumEt = ETW;
182
183 // components should all be less than 12 bits at this point with 200 MeV LSB
184 int MET2 = MET_x * MET_x + MET_y * MET_y;
185
186 if (MET2 > 0x0FFFFFF) {
187 MET = 0x000FFF;
188 } else {
189 // repeat the byte stream converter calculation here -- not what the hardware actually does
190 MET = std::sqrt(MET2);
191
192
193 // best guess at current hardware. Note that this is computed in the bytestream converter
194 // take most signficant 12 bits
195 //int MET12 = MET2 >> 12;
196 // simulate the look up -- only 6 most signficant bits currently set -- to be checked
197 //MET = ( (int)(std::sqrt(MET12)) << 6) & 0x00000FC0 ;
198 }
199
200
201 //Define a vector to be filled with all the TOBs of one event
202 std::vector<std::unique_ptr<gFEXJwoJTOB>> tobs_v;
203 tobs_v.resize(4);
204
205
206 // fill in TOBs
207 // The order of the TOBs is given according to the TOB ID (TODO: check how it's done in fw)
208
209 // First TOB is (MET, SumEt)
210 outTOB[0] = (total_sumEt & 0x00000FFF) << 0; //set the Quantity2 to the corresponding slot (LSB)
211 outTOB[0] = outTOB[0] | (MET & 0x00000FFF) << 12;//Quantity 1 (in bit number 12)
212 if (total_sumEt != 0) outTOB[0] = outTOB[0] | 0x00000001 << 24;//Status bit for Quantity 2 (0 if quantity is null)
213 if (MET != 0) outTOB[0] = outTOB[0] | 0x00000001 << 25;//Status bit for Quantity 1 (0 if quantity is null)
214 outTOB[0] = outTOB[0] | (1 & 0x0000001F) << 26;//TOB ID is 1 for scalar values (5 bits starting at 26)
215
216 // std::cout << "DMS MET " << std::hex << MET << " total_sumEt " << total_sumEt << std::endl << std::dec;
217
218// Second TOB is (MET_x, MET_y)
219 outTOB[1] = (MET_y & 0x00000FFF) << 0; //set the Quantity2 to the corresponding slot (LSB)
220 outTOB[1] = outTOB[1] | (MET_x & 0x00000FFF) << 12;//Quantity 1 (in bit number 12)
221 if (MET_y != 0) outTOB[1] = outTOB[1] | 0x00000001 << 24;//Status bit for Quantity 2 (0 if quantity is null)
222 if (MET_x != 0) outTOB[1] = outTOB[1] | 0x00000001 << 25;//Status bit for Quantity 1 (0 if quantity is null)
223 outTOB[1] = outTOB[1] | (2 & 0x0000001F) << 26;//TOB ID is 2 for MET_x, MET_y (5 bits starting at 26)
224
225// Third TOB is hard components (MHT_x, MHT_y)
226 outTOB[2] = (MHT_y & 0x00000FFF) << 0; //set the Quantity2 to the corresponding slot (LSB)
227 outTOB[2] = outTOB[2] | (MHT_x & 0x00000FFF) << 12;//Quantity 1 (in bit number 12)
228 if (MHT_y != 0) outTOB[2] = outTOB[2] | 0x00000001 << 24;//Status bit for Quantity 2 (0 if quantity is null)
229 if (MHT_x != 0) outTOB[2] = outTOB[2] | 0x00000001 << 25;//Status bit for Quantity 1 (0 if quantity is null)
230 outTOB[2] = outTOB[2] | (3 & 0x0000001F) << 26;//TOB ID is 3 for hard components (5 bits starting at 26)
231
232 // Fourth TOB is hard components (MST_x, MST_y)
233 outTOB[3] = (MST_y & 0x00000FFF) << 0; //set the Quantity2 to the corresponding slot (LSB)
234 outTOB[3] = outTOB[3] | (MST_x & 0x00000FFF) << 12;//Quantity 1 (in bit number 12)
235 if (MST_y != 0) outTOB[3] = outTOB[3] | 0x00000001 << 24;//Status bit for Quantity 2 (0 if quantity is null)
236 if (MST_x != 0) outTOB[3] = outTOB[3] | 0x00000001 << 25;//Status bit for Quantity 1 (0 if quantity is null)
237 outTOB[3] = outTOB[3] | (4 & 0x0000001F) << 26;//TOB ID is 4 for soft components (5 bits starting at 26)
238
239
240 tobs_v[0] = std::make_unique<gFEXJwoJTOB>();
241 tobs_v[0]->setWord(outTOB[0]);
242 tobs_v[0]->setQuantity1(MET);
243 tobs_v[0]->setQuantity2(total_sumEt);
244 tobs_v[0]->setSaturation(0); //Always 0 for now, need a threshold?
245 tobs_v[0]->setTobID(1);
246 if( MET != 0 ) tobs_v[0]->setStatus1(1);
247 else tobs_v[0]->setStatus1(0);
248 if(total_sumEt!= 0) tobs_v[0]->setStatus2(1);
249 else tobs_v[0]->setStatus2(0);
250
251 tobs_v[1] = std::make_unique<gFEXJwoJTOB>();
252 tobs_v[1]->setWord(outTOB[1]);
253 tobs_v[1]->setQuantity1(MET_x);
254 tobs_v[1]->setQuantity2(MET_y);
255 tobs_v[1]->setSaturation(0); //Always 0 for now, need a threshold?
256 tobs_v[1]->setTobID(2);
257 if( MET_x != 0 ) tobs_v[1]->setStatus1(1);
258 else tobs_v[1]->setStatus1(0);
259 if(MET_y!= 0) tobs_v[1]->setStatus2(1);
260 else tobs_v[1]->setStatus2(0);
261
262 tobs_v[2] = std::make_unique<gFEXJwoJTOB>();
263 tobs_v[2]->setWord(outTOB[2]);
264 tobs_v[2]->setQuantity1(MHT_x);
265 tobs_v[2]->setQuantity2(MHT_y);
266 tobs_v[2]->setSaturation(0); //Always 0 for now, need a threshold?
267 tobs_v[2]->setTobID(3);
268 if( MHT_x != 0 ) tobs_v[2]->setStatus1(1);
269 else tobs_v[2]->setStatus1(0);
270 if(MHT_y!= 0) tobs_v[2]->setStatus2(1);
271 else tobs_v[2]->setStatus2(0);
272
273 tobs_v[3] = std::make_unique<gFEXJwoJTOB>();
274 tobs_v[3]->setWord(outTOB[3]);
275 tobs_v[3]->setQuantity1(MST_x);
276 tobs_v[3]->setQuantity2(MST_y);
277 tobs_v[3]->setSaturation(0); //Always 0 for now, need a threshold?
278 tobs_v[3]->setTobID(4);
279 if( MST_x != 0 ) tobs_v[3]->setStatus1(1);
280 else tobs_v[2]->setStatus1(0);
281 if(MST_y!= 0) tobs_v[3]->setStatus2(1);
282 else tobs_v[3]->setStatus2(0);
283
284
285 return tobs_v;
286
287}
288
289
290
291void gFEXJwoJAlgo::gBlockAB(const gTowersType & twrs, gTowersType & gBlkSum, gTowersType & hasSeed, int seedThreshold) const {
292
293 int rows = twrs.size();
294 int cols = twrs[0].size();
295 for( int irow = 0; irow < rows; irow++ ){
296 for(int jcolumn = 0; jcolumn<cols; jcolumn++){
297 // zero jet sum here
298 gBlkSum[irow][jcolumn] = 0;
299 int krowUp = (irow + 1)%32;
300 int krowDn = (irow - 1 +32)%32;
301 if( (jcolumn == 0) || (jcolumn == 6) ) {
302 //left edge case
303 gBlkSum[irow][jcolumn] =
304 twrs[irow][jcolumn] + twrs[krowUp][jcolumn] + twrs[krowDn][jcolumn] +
305 twrs[irow][jcolumn+1] + twrs[krowUp][jcolumn+1] + twrs[krowDn][jcolumn+1];
306 } else if( (jcolumn == 5) || (jcolumn == 11)) {
307 // right edge case
308 gBlkSum[irow][jcolumn] =
309 twrs[irow][jcolumn] + twrs[krowUp][jcolumn] + twrs[krowDn][jcolumn] +
310 twrs[irow][jcolumn-1] + twrs[krowUp][jcolumn-1] + twrs[krowDn][jcolumn-1];
311 } else{
312 // normal case
313 gBlkSum[irow][jcolumn] =
314 twrs[irow][jcolumn] + twrs[krowUp][jcolumn] + twrs[krowDn][jcolumn] +
315 twrs[irow][jcolumn-1] + twrs[krowUp][jcolumn-1] + twrs[krowDn][jcolumn-1] +
316 twrs[irow][jcolumn+1] + twrs[krowUp][jcolumn+1] + twrs[krowDn][jcolumn+1];
317 }
318
319 if( gBlkSum[irow][jcolumn] > seedThreshold) {
320 hasSeed[irow][jcolumn] = 1;
321 } else {
322 hasSeed[irow][jcolumn] = 0;
323 }
324
325 if ( gBlkSum[irow][jcolumn] < 0 )
326 gBlkSum[irow][jcolumn] = 0;
327
328 // was bits 11+3 downto 3, now is 11 downto 0
329 if ( gBlkSum[irow][jcolumn] > FEXAlgoSpaceDefs::gBlockMax ) {
330 gBlkSum[irow][jcolumn] = FEXAlgoSpaceDefs::gBlockMax;
331 }
332 }
333 }
334}
335void gFEXJwoJAlgo::metFPGA_rho(int FPGAnum, const gTowersType& twrs, int puc_jwj,
336 const gTowersType & gBlkSum, int gBlockthreshold,
337 int aFPGA, int bFPGA,
338 int & MHT_x, int & MHT_y,
339 int & MST_x, int & MST_y,
340 int & MET_x, int & MET_y) const {
341 gBlockthreshold = gBlockthreshold * 200 / 800;
342
343 int64_t h_tx_hi = 0;
344 int64_t h_ty_hi = 0;
345 int64_t h_tx_lw = 0;
346 int64_t h_ty_lw = 0;
347
348 int64_t e_tx_hi = 0;
349 int64_t e_ty_hi = 0;
350 int64_t e_tx_lw = 0;
351 int64_t e_ty_lw = 0;
352
353
354 int64_t RHO_SUM_OF_COS_h_tx_hi = 0;
355 int64_t RHO_SUM_OF_SIN_h_ty_hi = 0;
356 int64_t RHO_SUM_OF_COS_h_tx_lw = 0;
357 int64_t RHO_SUM_OF_SIN_h_ty_lw = 0;
358
359 int64_t RHO_SUM_OF_COS_e_tx_hi = 0;
360 int64_t RHO_SUM_OF_SIN_e_ty_hi = 0;
361 int64_t RHO_SUM_OF_COS_e_tx_lw = 0;
362 int64_t RHO_SUM_OF_SIN_e_ty_lw = 0;
363
364
365 for( int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
366 for(int jcolumn = 6; jcolumn<12; jcolumn++){
367 if( FPGAnum == 2){
368 int frow = 2*(irow/2) + 1;
369
370 if(gBlkSum[irow][jcolumn] > gBlockthreshold){
371 h_tx_hi += (twrs[irow][jcolumn])*(cosLUT(frow, 5));
372 h_ty_hi += (twrs[irow][jcolumn])*(sinLUT(frow, 5));
373 RHO_SUM_OF_COS_h_tx_hi += (cosLUT(frow, 5));
374 RHO_SUM_OF_SIN_h_ty_hi += (sinLUT(frow, 5));
375
376 } else {
377 e_tx_hi += (twrs[irow][jcolumn])*(cosLUT(frow, 5));
378 e_ty_hi += (twrs[irow][jcolumn])*(sinLUT(frow, 5));
379 RHO_SUM_OF_COS_e_tx_hi += (cosLUT(frow, 5));
380 RHO_SUM_OF_SIN_e_ty_hi += (sinLUT(frow, 5));
381 }
382
383 } else {
384
385 if(gBlkSum[irow][jcolumn] > gBlockthreshold){
386 h_tx_hi += (twrs[irow][jcolumn])*(cosLUT(irow, 5));
387 h_ty_hi += (twrs[irow][jcolumn])*(sinLUT(irow, 5));
388 RHO_SUM_OF_COS_h_tx_hi += (cosLUT(irow, 5));
389 RHO_SUM_OF_SIN_h_ty_hi += (sinLUT(irow, 5));
390 } else {
391 e_tx_hi += (twrs[irow][jcolumn])*(cosLUT(irow, 5));
392 e_ty_hi += (twrs[irow][jcolumn])*(sinLUT(irow, 5));
393 RHO_SUM_OF_COS_e_tx_hi += (cosLUT(irow, 5));
394 RHO_SUM_OF_SIN_e_ty_hi += (sinLUT(irow, 5));
395 }
396 }
397 }
398
399 for(int jcolumn = 0; jcolumn<6; jcolumn++){
400 if( FPGAnum == 2){
401 int frow = 2*(irow/2) + 1;
402
403 if(gBlkSum[irow][jcolumn] > gBlockthreshold){
404 h_tx_lw += (twrs[irow][jcolumn])*(cosLUT(frow, 5));
405 h_ty_lw += (twrs[irow][jcolumn])*(sinLUT(frow, 5));
406 RHO_SUM_OF_COS_h_tx_lw += (cosLUT(frow, 5));
407 RHO_SUM_OF_SIN_h_ty_lw += (sinLUT(frow, 5));
408 } else{
409 e_tx_lw += (twrs[irow][jcolumn])*(cosLUT(frow, 5));
410 e_ty_lw += (twrs[irow][jcolumn])*(sinLUT(frow, 5));
411 RHO_SUM_OF_COS_e_tx_lw += (cosLUT(frow, 5));
412 RHO_SUM_OF_SIN_e_ty_lw += (sinLUT(frow, 5));
413 }
414 } else {
415
416 if(gBlkSum[irow][jcolumn] > gBlockthreshold){
417 h_tx_lw += (twrs[irow][jcolumn])*(cosLUT(irow, 5));
418 h_ty_lw += (twrs[irow][jcolumn])*(sinLUT(irow, 5));
419 RHO_SUM_OF_COS_h_tx_lw += (cosLUT(irow, 5));
420 RHO_SUM_OF_SIN_h_ty_lw += (sinLUT(irow, 5));
421 } else {
422 e_tx_lw += (twrs[irow][jcolumn])*(cosLUT(irow, 5));
423 e_ty_lw += (twrs[irow][jcolumn])*(sinLUT(irow, 5));
424 RHO_SUM_OF_COS_e_tx_lw += (cosLUT(irow, 5));
425 RHO_SUM_OF_SIN_e_ty_lw += (sinLUT(irow, 5));
426 }
427 }
428 }
429 }
430
431 // REMEMBER TO DO BIT ADJUSTMENTS FOR SUBTRACTION
432
433 long int fMHT_x = (h_tx_hi + h_tx_lw) ;
434 long int fMHT_y = (h_ty_hi + h_ty_lw) ;
435 long int fMST_x = (e_tx_hi + e_tx_lw) ;
436 long int fMST_y = (e_ty_hi + e_ty_lw) ;
437
438 long int RHO_MULTIPLIED_BY_SUM_OF_COS_HARD_RESULT_hi = ( puc_jwj * (RHO_SUM_OF_COS_h_tx_hi) ) >> 4 ; // could be >> 2, or >> 14
439 long int RHO_MULTIPLIED_BY_SUM_OF_SIN_HARD_RESULT_hi = ( puc_jwj * (RHO_SUM_OF_SIN_h_ty_hi) ) >> 4 ;
440 long int RHO_MULTIPLIED_BY_SUM_OF_COS_SOFT_RESULT_hi = ( puc_jwj * (RHO_SUM_OF_COS_e_tx_hi) ) >> 4 ;
441 long int RHO_MULTIPLIED_BY_SUM_OF_SIN_SOFT_RESULT_hi = ( puc_jwj * (RHO_SUM_OF_SIN_e_ty_hi) ) >> 4 ;
442
443
444 long int RHO_MULTIPLIED_BY_SUM_OF_COS_HARD_RESULT_lw = ( puc_jwj * (RHO_SUM_OF_COS_h_tx_lw) ) >> 4 ;
445 long int RHO_MULTIPLIED_BY_SUM_OF_SIN_HARD_RESULT_lw = ( puc_jwj * (RHO_SUM_OF_SIN_h_ty_lw) ) >> 4 ;
446 long int RHO_MULTIPLIED_BY_SUM_OF_COS_SOFT_RESULT_lw = ( puc_jwj * (RHO_SUM_OF_COS_e_tx_lw) ) >> 4 ;
447 long int RHO_MULTIPLIED_BY_SUM_OF_SIN_SOFT_RESULT_lw = ( puc_jwj * (RHO_SUM_OF_SIN_e_ty_lw) ) >> 4 ;
448
449 long int RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_h_tx_hi = (h_tx_hi - RHO_MULTIPLIED_BY_SUM_OF_COS_HARD_RESULT_hi) ;
450 long int RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_h_ty_hi = (h_ty_hi - RHO_MULTIPLIED_BY_SUM_OF_SIN_HARD_RESULT_hi) ;
451 long int RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_e_tx_hi = (e_tx_hi - RHO_MULTIPLIED_BY_SUM_OF_COS_SOFT_RESULT_hi) ;
452 long int RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_e_ty_hi = (e_ty_hi - RHO_MULTIPLIED_BY_SUM_OF_SIN_SOFT_RESULT_hi) ;
453
454 long int RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_h_tx_lw = (h_tx_lw - RHO_MULTIPLIED_BY_SUM_OF_COS_HARD_RESULT_lw) ;
455 long int RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_h_ty_lw = (h_ty_lw - RHO_MULTIPLIED_BY_SUM_OF_SIN_HARD_RESULT_lw) ;
456 long int RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_e_tx_lw = (e_tx_lw - RHO_MULTIPLIED_BY_SUM_OF_COS_SOFT_RESULT_lw) ;
457 long int RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_e_ty_lw = (e_ty_lw - RHO_MULTIPLIED_BY_SUM_OF_SIN_SOFT_RESULT_lw) ;
458
459 MHT_x = (RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_h_tx_hi + RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_h_tx_lw) >> 3;
460 MHT_y = (RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_h_ty_hi + RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_h_ty_lw) >> 3;
461 MST_x = (RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_e_tx_hi + RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_e_tx_lw) >> 3;
462 MST_y = (RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_e_ty_hi + RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_e_ty_lw) >> 3;
463
464 fMHT_x = (RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_h_tx_hi + RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_h_tx_lw) ;
465 fMHT_y = (RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_h_ty_hi + RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_h_ty_lw) ;
466 fMST_x = (RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_e_tx_hi + RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_e_tx_lw) ;
467 fMST_y = (RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_e_ty_hi + RHO_SUBTRACTED_BEFORE_FINAL_MULTIPLY_e_ty_lw) ;
468
469 long int fMET_x = ( aFPGA * (fMHT_x) + bFPGA * (fMST_x) ) >> 13 ;
470 long int fMET_y = ( aFPGA * (fMHT_y) + bFPGA * (fMST_y) ) >> 13 ;
471
472 MET_x = fMET_x;
473 MET_y = fMET_y;
474
475}
476
477
478void gFEXJwoJAlgo::metFPGA(int FPGAnum, const gTowersType& twrs,
479 const gTowersType & gBlkSum, int gBlockthreshold,
480 int aFPGA, int bFPGA,
481 int & MHT_x, int & MHT_y,
482 int & MST_x, int & MST_y,
483 int & MET_x, int & MET_y) const {
484
485 gBlockthreshold = gBlockthreshold * 200 / 800; //gBlockthreshold is provided in counts with a resolution of 200 MeV, but here needs to be applied with a resolution of 800 GeV
486 // in the RTL code these are 19+ 5 = 24 bits
487 int64_t h_tx_hi = 0;
488 int64_t h_ty_hi = 0;
489 int64_t h_tx_lw = 0;
490 int64_t h_ty_lw = 0;
491
492 int64_t e_tx_hi = 0;
493 int64_t e_ty_hi = 0;
494 int64_t e_tx_lw = 0;
495 int64_t e_ty_lw = 0;
496
497 for( int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
498 for(int jcolumn = 6; jcolumn<12; jcolumn++){
499 if( FPGAnum == 2){
500 int frow = 2*(irow/2) + 1;
501
502 if(gBlkSum[irow][jcolumn] > gBlockthreshold){
503 h_tx_hi += (twrs[irow][jcolumn])*(cosLUT(frow, 5));
504 h_ty_hi += (twrs[irow][jcolumn])*(sinLUT(frow, 5));
505 } else {
506 e_tx_hi += (twrs[irow][jcolumn])*(cosLUT(frow, 5));
507 e_ty_hi += (twrs[irow][jcolumn])*(sinLUT(frow, 5));
508 }
509
510 } else {
511
512 if(gBlkSum[irow][jcolumn] > gBlockthreshold){
513 h_tx_hi += (twrs[irow][jcolumn])*(cosLUT(irow, 5));
514 h_ty_hi += (twrs[irow][jcolumn])*(sinLUT(irow, 5));
515 } else {
516 e_tx_hi += (twrs[irow][jcolumn])*(cosLUT(irow, 5));
517 e_ty_hi += (twrs[irow][jcolumn])*(sinLUT(irow, 5));
518 }
519 }
520 }
521
522 for(int jcolumn = 0; jcolumn<6; jcolumn++){
523 if( FPGAnum == 2){
524 int frow = 2*(irow/2) + 1;
525
526 if(gBlkSum[irow][jcolumn] > gBlockthreshold){
527 h_tx_lw += (twrs[irow][jcolumn])*(cosLUT(frow, 5));
528 h_ty_lw += (twrs[irow][jcolumn])*(sinLUT(frow, 5));
529 } else{
530 e_tx_lw += (twrs[irow][jcolumn])*(cosLUT(frow, 5));
531 e_ty_lw += (twrs[irow][jcolumn])*(sinLUT(frow, 5));
532 }
533 } else {
534
535 if(gBlkSum[irow][jcolumn] > gBlockthreshold){
536 h_tx_lw += (twrs[irow][jcolumn])*(cosLUT(irow, 5));
537 h_ty_lw += (twrs[irow][jcolumn])*(sinLUT(irow, 5));
538 } else {
539 e_tx_lw += (twrs[irow][jcolumn])*(cosLUT(irow, 5));
540 e_ty_lw += (twrs[irow][jcolumn])*(sinLUT(irow, 5));
541 }
542 }
543 }
544 }
545
546 // According to https://gitlab.cern.ch/atlas-l1calo/gfex/firmware/-/issues/406#note_5662344
547 // there is no division -- so LSB is indeed 25 MeV
548
549 //but later changed to 200 MeV so divide by 8
550
551 long int fMHT_x = (h_tx_hi + h_tx_lw) ;
552 long int fMHT_y = (h_ty_hi + h_ty_lw) ;
553 long int fMST_x = (e_tx_hi + e_tx_lw) ;
554 long int fMST_y = (e_ty_hi + e_ty_lw) ;
555
556 MHT_x = (h_tx_hi + h_tx_lw) >> 3;
557 MHT_y = (h_ty_hi + h_ty_lw) >> 3;
558 MST_x = (e_tx_hi + e_tx_lw) >> 3;
559 MST_y = (e_ty_hi + e_ty_lw) >> 3;
560
561 // a and b coffecients are 10 bits
562 // multiplication has an addtional 2^10
563 // constant JWJ_OW : integer := 35;--Out width
564 // values are 35 bits long and top 16 bits are taken -- so divide by 2^19
565 // 2^10/2^19 = 1/2^9 = 1/512
566
567 long int fMET_x = ( aFPGA * (fMHT_x) + bFPGA * (fMST_x) ) >> 13 ;
568 long int fMET_y = ( aFPGA * (fMHT_y) + bFPGA * (fMST_y) ) >> 13 ;
569
570 MET_x = fMET_x;
571 MET_y = fMET_y;
572
573}
574
575void gFEXJwoJAlgo::etFPGA(int FPGAnum, const gTowersType& twrs, gTowersType &gBlkSum,
576 int gBlockthreshold, int A, int B, int &eth, int &ets, int &etw) const {
577
578 gBlockthreshold = gBlockthreshold * 200 / 800; //gBlockthreshold is provided in counts with a resolution of 200 MeV, but here needs to be applied with a resolution of 800 GeV
579
580 int64_t ethard_hi = 0;
581 int64_t etsoft_hi = 0;
582 int64_t ethard_lo = 0;
583 int64_t etsoft_lo = 0;
584
585 int64_t ethard = 0.0;
586 int64_t etsoft = 0.0;
587
588 int multiplicitiveFactor = 0;
589
590 if(FPGAnum < 2 ) {
591 multiplicitiveFactor = cosLUT(0, 5);
592 } else{
593 multiplicitiveFactor = cosLUT(1, 5);
594 }
595
596// firmware treats upper and lower columns differnetly
597
598 for( int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
599 for(int jcolumn = 0; jcolumn<6; jcolumn++){
600 if(gBlkSum[irow][jcolumn] > gBlockthreshold){
601 ethard_lo = ethard_lo + twrs[irow][jcolumn]*multiplicitiveFactor;
602 } else {
603 etsoft_lo = etsoft_lo + twrs[irow][jcolumn]*multiplicitiveFactor;
604 }
605 }
606 }
607
608 for( int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
609 for(int jcolumn = 6; jcolumn<12; jcolumn++){
610 if(gBlkSum[irow][jcolumn] > gBlockthreshold){
611 ethard_hi = ethard_hi + twrs[irow][jcolumn]*multiplicitiveFactor;
612 } else {
613 etsoft_hi = etsoft_hi + twrs[irow][jcolumn]*multiplicitiveFactor;
614 }
615 }
616 }
617
618 ethard = ethard_hi + ethard_lo;
619 etsoft = etsoft_hi + etsoft_lo;
620
621
622 int64_t etsum_hi = ethard_hi*A + etsoft_hi*B ;
623 if ( etsum_hi < 0 ) etsum_hi = 0;
624
625 int64_t etsum_lo = ethard_lo*A + etsoft_lo*B ;
626 if ( etsum_lo < 0 ) etsum_lo = 0;
627
628 int64_t etsum = etsum_hi + etsum_lo;
629
630
631 // convert 200 MeV LSB here
632 eth = ethard>>3;
633 ets = etsoft>>3;
634 etw = (etsum >>13 ) ;
635
636 if( etw < 0 ) etw = 0;
637 // max value is 15 bits with 800 MeV LSB -- so 17 bits here
638 if( etw > 0X001FFFF ) etw = 0X001FFFF ;
639
640
641 if(msgLvl(MSG::DEBUG)) {
642 std::cout << "DMS FPGA gTEJWOJ " << std::hex << FPGAnum << "et sum hard " << eth << "etsum soft" << ets << " A " << A << " B " << B << " weighted term " << etw << std::endl << std::dec;
643 }
644}
645
646void gFEXJwoJAlgo::etFastFPGA(int FPGAnum, const gTowersType& twrs, gTowersType &gBlkSum,
647 int gBlockthreshold, int A, int B, int &eth, int &ets, int &etw) const {
648
649 gBlockthreshold = gBlockthreshold * 200 / 800; //gBlockthreshold is provided in counts with a resolution of 200 MeV, but here needs to be applied with a resolution of 800 GeV
650 int64_t ethard_hi = 0;
651 int64_t etsoft_hi = 0;
652 int64_t ethard_lo = 0;
653 int64_t etsoft_lo = 0;
654
655 int64_t ethard = 0.0;
656 int64_t etsoft = 0.0;
657
658 // firmware treats upper and lower columns differently
659 for( int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
660 for(int jcolumn = 0; jcolumn<6; jcolumn++){
661 if(gBlkSum[irow][jcolumn] > gBlockthreshold){
662 ethard_lo = ethard_lo + twrs[irow][jcolumn];
663 }
664 else {
665 etsoft_lo = etsoft_lo + twrs[irow][jcolumn];
666 }
667 }
668 }
669
670 for( int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
671 for(int jcolumn = 6; jcolumn<12; jcolumn++){
672 if(gBlkSum[irow][jcolumn] > gBlockthreshold){
673 ethard_hi = ethard_hi + twrs[irow][jcolumn];
674 }
675 else {
676 etsoft_hi = etsoft_hi + twrs[irow][jcolumn];
677 }
678 }
679 }
680
681 ethard = ethard_hi + ethard_lo;
682 etsoft = etsoft_hi + etsoft_lo;
683
684 // convert 200 MeV LSB here
685 eth = ethard;
686 ets = etsoft; // Keep for reference -- not used in etFast
687 etw = ethard; // For EtFast the weighted term is just the hard term
688
689 // 16 bits signed, set max and min
690 if( etw < -32768 ) etw = -32768;
691 if( etw > 32767 ) etw = 32767;
692
693 if(msgLvl(MSG::DEBUG)) {
694 std::cout << "DMS FPGA gTEJWOJ " << std::hex << FPGAnum << "et sum hard " << eth << "etsum soft" << ets << " A " << A << " B " << B << " weighted term " << etw << std::endl << std::dec;
695 }
696}
697
698void gFEXJwoJAlgo::metTotal(int A_MET_x, int A_MET_y,
699 int B_MET_x, int B_MET_y,
700 int C_MET_x, int C_MET_y,
701 int & MET_x, int & MET_y) const {
702
703
704 MET_x = A_MET_x + B_MET_x + C_MET_x;
705 MET_y = A_MET_y + B_MET_y+ C_MET_y;
706
707 // Truncation of the result, as the individual quantities are 16 bits, while the TOB field is 12 bits
708 // MET_x = MET_x >> 4;
709 // MET_y = MET_y >> 4;
710
711 if (MET_x < -0x000800) MET_x = -0x000800; //-2048
712 if (MET_y < -0x000800) MET_y = -0x000800; //-2048
713
714 if (MET_x > 0x0007FF) MET_x = 0x0007FF; //2047
715 if (MET_y > 0x0007FF) MET_y = 0x0007FF; //2047
716
717}
718
720 int B_ET,
721 int C_ET,
722 int & ET ) const {
723
724 // leave at 200 MeV scale
725 if (A_ET < 0 ) A_ET = 0;
726 if (B_ET < 0 ) B_ET = 0;
727 if (C_ET < 0 ) C_ET = 0;
728
729 ET = (A_ET + B_ET + C_ET);
730
731 // main value of ET is always positive
732 if( ET > 0x0000FFF) ET = 0x0000FFF;
733
734}
735
736//----------------------------------------------------------------------------------
737// bitwise simulation of sine LUT in firmware
738//----------------------------------------------------------------------------------
739float gFEXJwoJAlgo::sinLUT(unsigned int phiIDX, unsigned int aw) const
740{
741 float c = static_cast<float>(phiIDX)/std::pow(2,aw);
742 float rad = (2*M_PI) *c;
743 float rsin = std::sin(rad);
744 int isin = std::round(rsin*(std::pow(2,aw) - 1));
745 return isin;
746
747}
748
749//----------------------------------------------------------------------------------
750// bitwise simulation cosine LUT in firmware
751//----------------------------------------------------------------------------------
752float gFEXJwoJAlgo::cosLUT(unsigned int phiIDX, unsigned int aw) const
753{
754 float c = static_cast<float>(phiIDX)/std::pow(2,aw);
755 float rad = (2*M_PI) *c;
756 float rcos = std::cos(rad);
757 int icos = std::round(rcos*(std::pow(2,aw) - 1));
758 return icos;
759}
760
761} // namespace LVL1
#define M_PI
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_ERROR(x)
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
bool msgLvl(const MSG::Level lvl) const
static constexpr int ABCrows
static constexpr int gBlockMax
static constexpr int ABcolumns
void etFastFPGA(int FPGAnum, const gTowersType &twrs, gTowersType &gBlkSum, int gBlockthreshold, int A, int B, int &eth, int &ets, int &etw) const
void metFPGA(int FPGAnum, const gTowersType &twrs, const gTowersType &gBlkSum, int gBlockthreshold, int aFPGA, int bFPGA, int &MHT_x, int &MHT_y, int &MST_x, int &MST_y, int &MET_x, int &MET_y) const
SG::ReadCondHandleKey< gFEXDBCondData > m_DBToolKey
void etFPGA(int FPGAnum, const gTowersType &twrs, gTowersType &gBlkSum, int gBlockthreshold, int A, int B, int &eth, int &ets, int &etw) const
virtual std::vector< std::unique_ptr< gFEXJwoJTOB > > jwojAlgo(const gTowersType &Atwr, int pucA_JWJ, const gTowersType &Btwr, int pucB_JWJ, const gTowersType &Ctwr, int pucC_JWJ, std::array< int32_t, 4 > &outTOB) const override
virtual StatusCode initialize() override
standard Athena-Algorithm method
float cosLUT(unsigned int phiIDX, unsigned int aw) const
void metFPGA_rho(int FPGAnum, const gTowersType &twrs, int puc_jwj, const gTowersType &gBlkSum, int gBlockthreshold, int aFPGA, int bFPGA, int &MHT_x, int &MHT_y, int &MST_x, int &MST_y, int &MET_x, int &MET_y) const
void metTotal(int A_MET_x, int A_MET_y, int B_MET_x, int B_MET_y, int C_MET_x, int C_MET_y, int &MET_x, int &MET_y) const
float sinLUT(unsigned int phiIDX, unsigned int aw) const
void gBlockAB(const gTowersType &twrs, gTowersType &gBlkSum, gTowersType &hasSeed, int seedThreshold) const
gFEXJwoJAlgo(const std::string &type, const std::string &name, const IInterface *parent)
Constructors.
void etTotal(int A_ET, int B_ET, int C_ET, int &ET) const
virtual void setAlgoConstant(int aFPGA_A, int bFPGA_A, int aFPGA_B, int bFPGA_B, int aFPGA_C, int bFPGA_C, int gXE_seedThrA, int gXE_seedThrB, int gXE_seedThrC) override
eFexTowerBuilder creates xAOD::eFexTowerContainer from supercells (LATOME) and triggerTowers (TREX) i...
std::array< std::array< int, 12 >, 32 > gTowersType
Definition IgFEXFPGA.h:25
Definition MET.py:1
hold the test vectors and ease the comparison