ATLAS Offline Software
gFEXJetAlgo.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 //***************************************************************************
5 // gFEXJetAlgo - JetFinder algorithm for gFEX
6 // -------------------
7 // begin : 01 04 2021
8 // email : cecilia.tosciri@cern.ch
9 //***************************************************************************
10 
11 #include <vector>
12 
13 #include "L1CaloFEXSim/gFEXJetAlgo.h"
14 #include "L1CaloFEXSim/gFEXJetTOB.h"
15 #include "L1CaloFEXSim/gTowerContainer.h"
16 #include "L1CaloFEXSim/gTower.h"
17 
18 namespace LVL1 {
19 
20  // default constructor for persistency
21 gFEXJetAlgo::gFEXJetAlgo(const std::string& type, const std::string& name, const IInterface* parent):
22  AthAlgTool(type, name, parent)
23  {
24  declareInterface<IgFEXJetAlgo>(this);
25  }
26 
27 
28 StatusCode gFEXJetAlgo::initialize(){
29 
30  ATH_CHECK(m_gFEXJetAlgo_gTowerContainerKey.initialize());
31 
32  return StatusCode::SUCCESS;
33 
34 }
35 
36 
37 
38 std::vector<std::unique_ptr<gFEXJetTOB>> gFEXJetAlgo::largeRfinder(
39  const gTowersType& Atwr,
40  const gTowersType& Btwr,
41  const gTowersType& Ctwr,
42  const gTowersType& Asat,
43  const gTowersType& Bsat,
44  const gTowersType& Csat,
45  int pucA, int pucB, int pucC, int gLJ_seedThrA, int gLJ_seedThrB, int gLJ_seedThrC,
46  int gJ_ptMinToTopoCounts1, int gJ_ptMinToTopoCounts2,
47  int jetThreshold, int gLJ_ptMinToTopoCounts1, int gLJ_ptMinToTopoCounts2,
48  std::array<uint32_t, 7> & ATOB1_dat, std::array<uint32_t, 7> & ATOB2_dat,
49  std::array<uint32_t, 7> & BTOB1_dat, std::array<uint32_t, 7> & BTOB2_dat,
50  std::array<uint32_t, 7> & CTOB1_dat, std::array<uint32_t, 7> & CTOB2_dat) const {
51 
52  // Arrays for gJets
53  // s = status (1 if above threshold)
54  // v = value 200 MeV LSB
55  // eta and phi are bin numbers
56  std::array<int, FEXAlgoSpaceDefs::gJetTOBfib> gTOBsat = {{0,0,0,0,0,0}};
57  std::array<int, FEXAlgoSpaceDefs::gJetTOBfib> gJetTOBs = {{0,0,0,0,0,0}};
58  std::array<int, FEXAlgoSpaceDefs::gJetTOBfib> gJetTOBv = {{0,0,0,0,0,0}};
59  std::array<int, FEXAlgoSpaceDefs::gJetTOBfib> gJetTOBeta = {{0,0,0,0,0,0}};
60  std::array<int, FEXAlgoSpaceDefs::gJetTOBfib> gJetTOBphi = {{0,0,0,0,0,0}};
61 
62  // Arrays for gBlocks (so far only 2 gBlocks per fiber)
63  // first index is for the fiber number (A 0 and 1, B 0 and 1, C P and C N )
64  // second index is for leading, sub-leading, and sub-sub-leading (not used)
65  std::array<std::array<int, 3>, FEXAlgoSpaceDefs::BTOBFIB> gBlockTOBv = {{ {{0,0,0}}, {{0,0,0}}, {{0,0,0}}, {{0,0,0}}, {{0,0,0}}, {{0,0,0}} }};
66  std::array<std::array<int, 3>, FEXAlgoSpaceDefs::BTOBFIB> gBlockTOBeta = {{ {{0,0,0}}, {{0,0,0}}, {{0,0,0}}, {{0,0,0}}, {{0,0,0}}, {{0,0,0}} }};
67  std::array<std::array<int, 3>, FEXAlgoSpaceDefs::BTOBFIB> gBlockTOBphi = {{ {{0,0,0}}, {{0,0,0}}, {{0,0,0}}, {{0,0,0}}, {{0,0,0}}, {{0,0,0}} }};
68 
69 
70  //Alternate format for FPGC C (split in two)
71 
72  gTowersJetEngine CNtwr;
73  gTowersJetEngine CPtwr;
74 
75  for(int irow=0; irow<FEXAlgoSpaceDefs::ABCrows; irow++){
76  for(int jcolumn=0;jcolumn<FEXAlgoSpaceDefs::ABCcolumnsEng;jcolumn++){
77  CNtwr[irow][jcolumn] = Ctwr[irow][jcolumn] ;
78  CPtwr[irow][jcolumn] = Ctwr[irow][jcolumn+FEXAlgoSpaceDefs::ABCcolumnsEng] ;
79  }
80 
81  for(unsigned int jcolumn=0;jcolumn<FEXAlgoSpaceDefs::ABcolumns;jcolumn++){
82 
83  if( jcolumn < 6 ) {
84  // TOB1
85  if( Asat[irow][jcolumn] == 1 ) gTOBsat[0] = 1;
86  if( Bsat[irow][jcolumn] == 1 ) gTOBsat[2] = 1;
87  if( Csat[irow][jcolumn] == 1 ) gTOBsat[4] = 1;
88  } else {
89  //TOB2
90  if( Asat[irow][jcolumn] == 1 ) gTOBsat[1] = 1;
91  if( Bsat[irow][jcolumn] == 1 ) gTOBsat[3] = 1;
92  if( Csat[irow][jcolumn] == 1 ) gTOBsat[5] = 1;
93  }
94  }
95  }
96 
97 
98 
99  gTowersType AjetsAlt;
100  singleHalf(Atwr,AjetsAlt);
101 
102  gTowersPartialSums AlpsAltOut;
103  gTowersPartialSums ArpsAltOut;
104 
105  InternalPartialAB(Atwr, AlpsAltOut, ArpsAltOut);
106 
107  gTowersType BjetsAlt;
108  singleHalf(Btwr,BjetsAlt);
109 
110  gTowersPartialSums BlpsAltOut;
111  gTowersPartialSums BrpsAltOut;
112 
113  InternalPartialAB(Btwr, BlpsAltOut, BrpsAltOut);
114 
115  gTowersType CjetsAlt;
116  singleHalf(Ctwr,CjetsAlt);
117 
118  // Add the left local partial sum (zero if right column)
119  addInternalLin(AjetsAlt, AlpsAltOut);
120  addInternalLin(BjetsAlt, BlpsAltOut);
121 
122  pileUpCorrectionAB(AjetsAlt,pucA);
123  pileUpCorrectionAB(BjetsAlt,pucB);
124  pileUpCorrectionAB(CjetsAlt,pucC);
125 
126  // If the local sum is less than zero set it to zero
127  ZeroNegative(AjetsAlt);
128  ZeroNegative(BjetsAlt);
129  ZeroNegative(CjetsAlt);
130 
131  // Add the right local partial sum (zero if right column)
132  addInternalRin(AjetsAlt, ArpsAltOut);
133  addInternalRin(BjetsAlt, BrpsAltOut);
134 
135  for(unsigned int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
136  for(unsigned int icolumn =0; icolumn<FEXAlgoSpaceDefs::ABcolumns; icolumn++){
137  // set 18 bits on
138  if (Asat[irow][icolumn] != 0){
139  }
140  }
141  }
142 
143  SaturateJets( AjetsAlt, Asat );
144  SaturateJets( BjetsAlt, Bsat );
145  SaturateJets( CjetsAlt, Csat );
146 
147 
148  gTowersType AjetsRestricted;
149  gTowersType BjetsRestricted;
150  gTowersType CjetsRestricted;
151 
152 
153  for(int irow=0; irow<FEXAlgoSpaceDefs::ABCrows; irow++){
154  for(int icolumn=0; icolumn<FEXAlgoSpaceDefs::ABcolumns; icolumn++){
155  AjetsRestricted[irow][icolumn] = AjetsAlt[irow][icolumn];
156  BjetsRestricted[irow][icolumn] = BjetsAlt[irow][icolumn];
157  CjetsRestricted[irow][icolumn] = CjetsAlt[irow][icolumn];
158  }
159  }
160 
161  // find gBlocks
162  gTowersType gBLKA;
163  gTowersType hasSeedA;
164  gBlockAB(Atwr, gBLKA, hasSeedA, gLJ_seedThrA);
165 
166  SaturateBlocks(gBLKA, Asat);
167 
168  gTowersType gBLKB;
169  gTowersType hasSeedB;
170  gBlockAB(Btwr, gBLKB, hasSeedB, gLJ_seedThrB);
171  SaturateBlocks(gBLKA, Bsat);
172 
173 
174  gTowersType gBLKC;
175  gTowersType hasSeedC;
176  gBlockAB(Ctwr, gBLKC, hasSeedC, gLJ_seedThrC);
177  SaturateBlocks(gBLKA, Csat);
178 
179  // sorting by jet engine -- not done in FPGA
180  std::array<int, 32> AgBlockOutL{};
181  std::array<int, 32> AgBlockEtaIndL{};
182  std::array<int, 32> AgBlockOutR{};
183  std::array<int, 32> AgBlockEtaIndR{};
184  blkOutAB(gBLKA, AgBlockOutL, AgBlockEtaIndL, AgBlockOutR, AgBlockEtaIndR);
185 
186  std::array<int, 32> BgBlockOutL{};
187  std::array<int, 32> BgBlockEtaIndL{};
188  std::array<int, 32> BgBlockOutR{};
189  std::array<int, 32> BgBlockEtaIndR{};
190  blkOutAB(gBLKB, BgBlockOutL, BgBlockEtaIndL, BgBlockOutR, BgBlockEtaIndR);
191 
192  std::array<int, 32> CgBlockOutL{};
193  std::array<int, 32> CgBlockEtaIndL{};
194  std::array<int, 32> CgBlockOutR{};
195  std::array<int, 32> CgBlockEtaIndR{};
196  blkOutAB(gBLKC, CgBlockOutL, CgBlockEtaIndL, CgBlockOutR, CgBlockEtaIndR) ;
197 
198  // sorting by 192 blocks 0 = left, 1 = right
199  // second index is column number, the eta value is 2 + 6*column number + local eta
200 
201  // find the leading and subleading gBlock in the left column of FPGA A
202  gBlockMax2(gBLKA, 1, 0, gBlockTOBv[0], gBlockTOBeta[0], gBlockTOBphi[0]);
203  // find the leading and subleading gBlock in the right column of FPGA A
204  gBlockMax2(gBLKA, 2, 1, gBlockTOBv[1], gBlockTOBeta[1], gBlockTOBphi[1]);
205 
206  // find the leading and subleading gBlock in the left column of FPGA B
207  gBlockMax2(gBLKB, 3, 0, gBlockTOBv[2], gBlockTOBeta[2], gBlockTOBphi[2]);
208  // find the leading and subleading gBlock in the right column of FPGA B
209  gBlockMax2(gBLKB, 4, 1, gBlockTOBv[3], gBlockTOBeta[3], gBlockTOBphi[3]);
210 
211  // find the leading and subleading gBlock in the CN column of FPGA C
212  gBlockMax2( gBLKC, 0, 0, gBlockTOBv[4], gBlockTOBeta[4], gBlockTOBphi[4]);
213  // find the leading and subleading gBlock in the CP column of FPGA C
214  gBlockMax2( gBLKC, 5, 1, gBlockTOBv[5], gBlockTOBeta[5], gBlockTOBphi[5]);
215 
216 
217  // in hardware vetos happen before remote sums come in.
218  gBlockVetoAB(AjetsRestricted, hasSeedA);
219  gBlockVetoAB(BjetsRestricted, hasSeedB);
220  gBlockVetoAB(CjetsRestricted, hasSeedC);
221 
222  // calculate A & B remote partial sums first
223  gTowersPartialSums RAlps_out, RArps_out;
224  RemotePartialAB(Atwr, RAlps_out, RArps_out);
225 
226  gTowersPartialSums RBlps_out, RBrps_out;
227  RemotePartialAB(Btwr, RBlps_out, RBrps_out);
228 
229  gTowersPartialSums RCNrps_out, RCPlps_out;
230  RemotePartialCN(CNtwr, RCNrps_out);
231  RemotePartialCP(CPtwr, RCPlps_out);
232 
233 
234  if( FEXAlgoSpaceDefs::ENABLE_INTER_AB ) {
235 
236  // input partial sums from FPGA B to A (lps_out -> rps_in)
237  addRemoteRin(AjetsRestricted, RBlps_out, FEXAlgoSpaceDefs::PS_UPPER_AB, FEXAlgoSpaceDefs::PS_LOWER_AB, FEXAlgoSpaceDefs::PS_SHIFT_AB);
238 
239  // input partial sums from FPGA B to A (lps_out -> rps_in)
240  addRemoteLin(BjetsRestricted, RArps_out, FEXAlgoSpaceDefs::PS_UPPER_AB, FEXAlgoSpaceDefs::PS_LOWER_AB, FEXAlgoSpaceDefs::PS_SHIFT_AB);
241 
242  }
243 
244  if( FEXAlgoSpaceDefs::ENABLE_INTER_C ) {
245 
246  addRemoteLin(AjetsRestricted, RCNrps_out, FEXAlgoSpaceDefs::PS_UPPER_C, FEXAlgoSpaceDefs::PS_LOWER_C, FEXAlgoSpaceDefs::PS_SHIFT_C);
247 
248  addRemoteRin(BjetsRestricted, RCPlps_out, FEXAlgoSpaceDefs::PS_UPPER_C, FEXAlgoSpaceDefs::PS_LOWER_C, FEXAlgoSpaceDefs::PS_SHIFT_C);
249 
250  }
251 
252 // this are added into the jets you have to truncate to the number of bits on the interFPGA communication
253 
254  if( FEXAlgoSpaceDefs::ENABLE_INTER_ABC ) {
255 
256  addRemoteCNin(CjetsRestricted, RAlps_out, FEXAlgoSpaceDefs::PS_UPPER_C, FEXAlgoSpaceDefs::PS_LOWER_C, FEXAlgoSpaceDefs::PS_SHIFT_C );
257 
258  addRemoteCPin(CjetsRestricted, RBrps_out, FEXAlgoSpaceDefs::PS_UPPER_C, FEXAlgoSpaceDefs::PS_LOWER_C, FEXAlgoSpaceDefs::PS_SHIFT_C );
259 
260  }
261 
262 
263  //Emulate switch to unsigned values by zeroing everything below the jet threshold
264  //https://gitlab.cern.ch/atlas-l1calo/gfex/firmware/-/blob/devel/common/jet_finder/HDL/jet_eng.vhd#L550
265 
266  gJetVetoAB(AjetsRestricted, jetThreshold);
267  gJetVetoAB(BjetsRestricted, jetThreshold);
268  gJetVetoAB(CjetsRestricted, jetThreshold);
269 
270 
271  std::array<int, 32> AjetOutL;
272  std::array<int, 32> AetaIndL;
273  std::array<int, 32> AjetOutR;
274  std::array<int, 32> AetaIndR;
275 
276 
277  jetOutAB(AjetsRestricted, AjetOutL, AetaIndL, AjetOutR, AetaIndR);
278 
279  std::array<int, 32> BjetOutL;
280  std::array<int, 32> BetaIndL;
281  std::array<int, 32> BjetOutR;
282  std::array<int, 32> BetaIndR;
283 
284  jetOutAB(BjetsRestricted, BjetOutL, BetaIndL, BjetOutR, BetaIndR);
285 
286  std::array<int, 32> CNjetOut;
287  std::array<int, 32> CNetaInd;
288  std::array<int, 32> CPjetOut;
289  std::array<int, 32> CPetaInd;
290 
291  jetOutAB(CjetsRestricted, CNjetOut, CNetaInd, CPjetOut, CPetaInd);
292 
293 
294 
295  gJetTOBgen(AjetOutL, AetaIndL, 0, jetThreshold, gJetTOBs, gJetTOBv, gJetTOBeta, gJetTOBphi);
296  gJetTOBgen(AjetOutR, AetaIndR, 1, jetThreshold, gJetTOBs, gJetTOBv, gJetTOBeta, gJetTOBphi);
297 
298  gJetTOBgen(BjetOutL, BetaIndL, 2, jetThreshold, gJetTOBs, gJetTOBv, gJetTOBeta, gJetTOBphi);
299  gJetTOBgen(BjetOutR, BetaIndR, 3, jetThreshold, gJetTOBs, gJetTOBv, gJetTOBeta, gJetTOBphi);
300 
301  gJetTOBgen(CNjetOut, CNetaInd, 4, jetThreshold, gJetTOBs, gJetTOBv, gJetTOBeta, gJetTOBphi);
302  gJetTOBgen(CPjetOut, CPetaInd, 5, jetThreshold, gJetTOBs, gJetTOBv, gJetTOBeta, gJetTOBphi);
303 
305  std::vector<std::unique_ptr<gFEXJetTOB>> tobs_v;
306  tobs_v.resize(14);
307 
308  // fill in TOBs
309  ATOB1_dat[0] = 0;
310  BTOB1_dat[0] = 0;
311  CTOB1_dat[0] = 0;
312 
313  ATOB2_dat[0] = (( pucA & 0x0000FFFF ) << 8);
314  BTOB2_dat[0] = (( pucB & 0x0000FFFF ) << 8);
315  CTOB2_dat[0] = (( pucC & 0x0000FFFF ) << 8);
316 
317  // //First available TOBs are the gRho for each central FPGA
318  tobs_v[0] = std::make_unique<gFEXJetTOB>();
319  tobs_v[0]->setWord(ATOB2_dat[0]);
320  tobs_v[0]->setET(pucA);
321  tobs_v[0]->setEta(0);
322  tobs_v[0]->setPhi(0);
323  tobs_v[0]->setTobID(0);
324  tobs_v[0]->setStatus(1);//is 1 only if rho calculation is valid, but rho is not calculated at the moment
325 
326  tobs_v[1] = std::make_unique<gFEXJetTOB>();
327  tobs_v[1]->setWord(BTOB2_dat[0]);
328  tobs_v[1]->setET(pucB);
329  tobs_v[1]->setEta(0);
330  tobs_v[1]->setPhi(0);
331  tobs_v[1]->setTobID(0);
332  tobs_v[1]->setStatus(1);//is 1 only if rho calculation is valid, but rho is not calculated at the moment
333 
334  // leading gBlocks (available first and go in first TOB value)
335  // TOBs 2-5 are leading gBlocks
336  ATOB1_dat[1] = 0x00000001; //set the TOB ID in the corresponding slot (LSB)
337  if( gBlockTOBv[0][0] > gJ_ptMinToTopoCounts2 ) ATOB1_dat[1] = ATOB1_dat[1] | 0x00000080;//status
338  ATOB1_dat[1] = ATOB1_dat[1] | ( ( gBlockTOBv[0][0] & 0x00000FFF ) << 8);
339  ATOB1_dat[1] = ATOB1_dat[1] | ( ( gBlockTOBeta[0][0] & 0x0000003F ) <<20);
340  ATOB1_dat[1] = ATOB1_dat[1] | ( ( gBlockTOBphi[0][0] & 0x0000001F ) <<26);
341 
342  tobs_v[2] = std::make_unique<gFEXJetTOB>();
343  tobs_v[2]->setWord(ATOB1_dat[1]);
344  tobs_v[2]->setET(gBlockTOBv[0][0]);
345  tobs_v[2]->setEta(gBlockTOBeta[0][0]);
346  tobs_v[2]->setPhi(gBlockTOBphi[0][0]);
347  tobs_v[2]->setTobID(1);
348  if( gBlockTOBv[0][0] > gJ_ptMinToTopoCounts2 ) tobs_v[2]->setStatus(1);
349  else tobs_v[2]->setStatus(0);
350 
351  ATOB2_dat[1] = 0x00000002;
352  if( gBlockTOBv[1][0] > gJ_ptMinToTopoCounts1 ) ATOB2_dat[1] = ATOB2_dat[1] | 0x00000080;
353  ATOB2_dat[1] = ATOB2_dat[1] | ( ( gBlockTOBv[1][0] & 0x00000FFF ) << 8);
354  ATOB2_dat[1] = ATOB2_dat[1] | ( ( gBlockTOBeta[1][0] & 0x0000003F ) <<20);
355  ATOB2_dat[1] = ATOB2_dat[1] | ( ( gBlockTOBphi[1][0] & 0x0000001F ) <<26);
356 
357  tobs_v[3] = std::make_unique<gFEXJetTOB>();
358  tobs_v[3]->setWord(ATOB2_dat[1]);
359  tobs_v[3]->setET(gBlockTOBv[1][0]);
360  tobs_v[3]->setEta(gBlockTOBeta[1][0]);
361  tobs_v[3]->setPhi(gBlockTOBphi[1][0]);
362  tobs_v[3]->setTobID(2);
363  if( gBlockTOBv[1][0] > gJ_ptMinToTopoCounts1 ) tobs_v[3]->setStatus(1);
364  else tobs_v[3]->setStatus(0);
365 
366 
367  BTOB1_dat[1] = 0x00000001;
368  if( gBlockTOBv[2][0] > gJ_ptMinToTopoCounts1 ) BTOB1_dat[1] = BTOB1_dat[1] | 0x00000080;
369  BTOB1_dat[1] = BTOB1_dat[1] | ( ( gBlockTOBv[2][0] & 0x00000FFF ) << 8);
370  BTOB1_dat[1] = BTOB1_dat[1] | ( ( gBlockTOBeta[2][0] & 0x0000003F ) <<20);
371  BTOB1_dat[1] = BTOB1_dat[1] | ( ( gBlockTOBphi[2][0] & 0x0000001F) <<26);
372 
373  tobs_v[4] = std::make_unique<gFEXJetTOB>();
374  tobs_v[4]->setWord(BTOB1_dat[1]);
375  tobs_v[4]->setET(gBlockTOBv[2][0]);
376  tobs_v[4]->setEta(gBlockTOBeta[2][0]);
377  tobs_v[4]->setPhi(gBlockTOBphi[2][0]);
378  tobs_v[4]->setTobID(1);
379  if( gBlockTOBv[2][0] > gJ_ptMinToTopoCounts1 ) tobs_v[4]->setStatus(1);
380  else tobs_v[4]->setStatus(0);
381 
382 
383  BTOB2_dat[1] = 0x00000002;
384  if( gBlockTOBv[3][0] > gJ_ptMinToTopoCounts2 ) BTOB2_dat[1] = BTOB2_dat[1] | 0x00000080;
385  BTOB2_dat[1] = BTOB2_dat[1] | ( ( gBlockTOBv[3][0] & 0x00000FFF ) << 8);
386  BTOB2_dat[1] = BTOB2_dat[1] | ( ( gBlockTOBeta[3][0] & 0x0000003F ) <<20);
387  BTOB2_dat[1] = BTOB2_dat[1] | ( ( gBlockTOBphi[3][0] & 0x0000001F ) <<26);
388 
389  tobs_v[5] = std::make_unique<gFEXJetTOB>();
390  tobs_v[5]->setWord(BTOB2_dat[1]);
391  tobs_v[5]->setET(gBlockTOBv[3][0]);
392  tobs_v[5]->setEta(gBlockTOBeta[3][0]);
393  tobs_v[5]->setPhi(gBlockTOBphi[3][0]);
394  tobs_v[5]->setTobID(2);
395  if( gBlockTOBv[3][0] > gJ_ptMinToTopoCounts2 ) tobs_v[5]->setStatus(1);
396  else tobs_v[5]->setStatus(0);
397 
398 
399  CTOB1_dat[1] = 0x00000001;
400  if( gBlockTOBv[4][0] > gJ_ptMinToTopoCounts1 ) CTOB1_dat[1] = CTOB1_dat[1] | 0x00000080;
401  CTOB1_dat[1] = CTOB1_dat[1] | ( ( gBlockTOBv[4][0] & 0x00000FFF ) << 8);
402  CTOB1_dat[1] = CTOB1_dat[1] | ( ( gBlockTOBeta[4][0] & 0x0000003F ) <<20);
403  CTOB1_dat[1] = CTOB1_dat[1] | ( ( gBlockTOBphi[4][0] & 0x0000001F) <<26);
404 
405 
406  CTOB2_dat[1] = 0x00000002;
407  if( gBlockTOBv[5][0] > gJ_ptMinToTopoCounts2) CTOB2_dat[1] = CTOB2_dat[1] | 0x00000080;
408  CTOB2_dat[1] = CTOB2_dat[1] | ( ( gBlockTOBv[5][0] & 0x00000FFF ) << 8);
409  CTOB2_dat[1] = CTOB2_dat[1] | ( ( gBlockTOBeta[5][0] & 0x0000003F ) <<20);
410  CTOB2_dat[1] = CTOB2_dat[1] | ( ( gBlockTOBphi[5][0] & 0x0000001F ) <<26);
411 
412 
413  // subleading gBlocks
414  // TOBs 6-9 are subleading gBlocks
415  ATOB1_dat[2] = 0x00000003;
416  if( gBlockTOBv[0][1] > gJ_ptMinToTopoCounts2 ) ATOB1_dat[2] = ATOB1_dat[2] | 0x00000080;
417  ATOB1_dat[2] = ATOB1_dat[2] | ( ( gBlockTOBv[0][1] & 0x00000FFF ) << 8);
418  ATOB1_dat[2] = ATOB1_dat[2] | ( ( gBlockTOBeta[0][1] & 0x0000003F ) <<20);
419  ATOB1_dat[2] = ATOB1_dat[2] | ( ( gBlockTOBphi[0][1] & 0x0000001F ) <<26);
420 
421  tobs_v[6] = std::make_unique<gFEXJetTOB>();
422  tobs_v[6]->setWord(ATOB1_dat[2]);
423  tobs_v[6]->setET(gBlockTOBv[0][1]);
424  tobs_v[6]->setEta(gBlockTOBeta[0][1]);
425  tobs_v[6]->setPhi(gBlockTOBphi[0][1]);
426  tobs_v[6]->setTobID(3);
427  if( gBlockTOBv[0][1] > gJ_ptMinToTopoCounts2 ) tobs_v[6]->setStatus(1);
428  else tobs_v[6]->setStatus(0);
429 
430 
431  ATOB2_dat[2] = 0x00000004;
432  if( gBlockTOBv[1][1] > gJ_ptMinToTopoCounts1 ) ATOB2_dat[2] = ATOB2_dat[2] | 0x00000080;
433  ATOB2_dat[2] = ATOB2_dat[2] | ( ( gBlockTOBv[1][1] & 0x00000FFF ) << 8);
434  ATOB2_dat[2] = ATOB2_dat[2] | ( ( gBlockTOBeta[1][1] & 0x0000003F ) <<20);
435  ATOB2_dat[2] = ATOB2_dat[2] | ( ( gBlockTOBphi[1][1] & 0x0000001F ) <<26);
436 
437  tobs_v[7] = std::make_unique<gFEXJetTOB>();
438  tobs_v[7]->setWord(ATOB2_dat[2]);
439  tobs_v[7]->setET(gBlockTOBv[1][1]);
440  tobs_v[7]->setEta(gBlockTOBeta[1][1]);
441  tobs_v[7]->setPhi(gBlockTOBphi[1][1]);
442  tobs_v[7]->setTobID(4);
443  if( gBlockTOBv[1][1] > gJ_ptMinToTopoCounts1 ) tobs_v[7]->setStatus(1);
444  else tobs_v[7]->setStatus(0);
445 
446 
447  BTOB1_dat[2] = 0x00000003;
448  if( gBlockTOBv[2][1] > gJ_ptMinToTopoCounts1 ) BTOB1_dat[2] = BTOB1_dat[2] | 0x00000080;
449  BTOB1_dat[2] = BTOB1_dat[2] | ( ( gBlockTOBv[2][1] & 0x00000FFF ) << 8);
450  BTOB1_dat[2] = BTOB1_dat[2] | ( ( gBlockTOBeta[2][1] & 0x0000003F ) <<20);
451  BTOB1_dat[2] = BTOB1_dat[2] | ( ( gBlockTOBphi[2][1] & 0x0000001F ) <<26);
452 
453  tobs_v[8] = std::make_unique<gFEXJetTOB>();
454  tobs_v[8]->setWord(BTOB1_dat[2]);
455  tobs_v[8]->setET(gBlockTOBv[2][1]);
456  tobs_v[8]->setEta(gBlockTOBeta[2][1]);
457  tobs_v[8]->setPhi(gBlockTOBphi[2][1]);
458  tobs_v[8]->setTobID(3);
459  if( gBlockTOBv[2][1] > gJ_ptMinToTopoCounts1 ) tobs_v[8]->setStatus(1);
460  else tobs_v[8]->setStatus(0);
461 
462 
463  BTOB2_dat[2] = 0x00000004;
464  if( gBlockTOBv[3][1] > gJ_ptMinToTopoCounts2 ) BTOB2_dat[2] = BTOB2_dat[2] | 0x00000080;
465  BTOB2_dat[2] = BTOB2_dat[2] | ( ( gBlockTOBv[3][1] & 0x00000FFF ) << 8);
466  BTOB2_dat[2] = BTOB2_dat[2] | ( ( gBlockTOBeta[3][1] & 0x0000003F ) <<20);
467  BTOB2_dat[2] = BTOB2_dat[2] | ( ( gBlockTOBphi[3][1] & 0x0000001F ) <<26);
468 
469  tobs_v[9] = std::make_unique<gFEXJetTOB>();
470  tobs_v[9]->setWord(BTOB2_dat[2]);
471  tobs_v[9]->setET(gBlockTOBv[3][1]);
472  tobs_v[9]->setEta(gBlockTOBeta[3][1]);
473  tobs_v[9]->setPhi(gBlockTOBphi[3][1]);
474  tobs_v[9]->setTobID(4);
475  if( gBlockTOBv[3][1] > gJ_ptMinToTopoCounts2 ) tobs_v[9]->setStatus(1);
476  else tobs_v[9]->setStatus(0);
477 
478 
479  CTOB1_dat[2] = 0x00000003;
480  if( gBlockTOBv[4][1] > gJ_ptMinToTopoCounts1 ) CTOB1_dat[2] = CTOB1_dat[2] | 0x00000080;
481  CTOB1_dat[2] = CTOB1_dat[2] | ( ( gBlockTOBv[4][1] & 0x00000FFF ) << 8);
482  CTOB1_dat[2] = CTOB1_dat[2] | ( ( gBlockTOBeta[4][1] & 0x0000003F ) <<20);
483  CTOB1_dat[2] = CTOB1_dat[2] | ( ( gBlockTOBphi[4][1] & 0x0000001F ) <<26);
484 
485  CTOB2_dat[2] = 0x00000004;
486  if( gBlockTOBv[5][1] > gJ_ptMinToTopoCounts2 ) CTOB2_dat[2] = CTOB2_dat[2] | 0x00000080;
487  CTOB2_dat[2] = CTOB2_dat[2] | ( ( gBlockTOBv[5][1] & 0x00000FFF ) << 8);
488  CTOB2_dat[2] = CTOB2_dat[2] | ( ( gBlockTOBeta[5][1] & 0x0000003F ) <<20);
489  CTOB2_dat[2] = CTOB2_dat[2] | ( ( gBlockTOBphi[5][1] & 0x0000001F ) <<26);
490 
491 
492  // finally the main event -- lead gJET
493  // according the specification https://docs.google.com/spreadsheets/d/15YVVtGofhXMtV7jXRFzWO0FVUtUAjS-X-aQjh3FKE_w/edit#gid=523371660
494  // we should have an lsb of 3.2 GeV -- so ignore lower 4 bits.
495 
496  // TOBs 10-13 are leading gJets
497  // shift is done before sorting as in firmware!
498  int tobvShift = 8;
499  int tobvMask = 0x00000FFF;
500 
501 
502  ATOB1_dat[3] = 0x00000005;
503  if( gJetTOBv[0] > gLJ_ptMinToTopoCounts2 ) ATOB1_dat[3] = ATOB1_dat[3] | 0x00000080;
504  ATOB1_dat[3] = ATOB1_dat[3] | ( ( gJetTOBv[0] & tobvMask) << tobvShift);
505  ATOB1_dat[3] = ATOB1_dat[3] | ( ( gJetTOBeta[0] & 0x0000003F ) <<20);
506  ATOB1_dat[3] = ATOB1_dat[3] | ( ( gJetTOBphi[0] & 0x0000001F ) <<26);
507 
508  tobs_v[10] = std::make_unique<gFEXJetTOB>();
509  tobs_v[10]->setWord(ATOB1_dat[3]);
510  tobs_v[10]->setET(gJetTOBv[0]);
511  tobs_v[10]->setEta(gJetTOBeta[0]);
512  tobs_v[10]->setPhi(gJetTOBphi[0]);
513  tobs_v[10]->setTobID(5);
514  if(gJetTOBv[0] > gLJ_ptMinToTopoCounts2 ) tobs_v[10]->setStatus(1);
515  else tobs_v[10]->setStatus(0);
516 
517 
518  ATOB2_dat[3] = 0x00000006;
519  if( gJetTOBv[1] > gLJ_ptMinToTopoCounts1 ) ATOB2_dat[3] = ATOB2_dat[3] | 0x00000080;
520  ATOB2_dat[3] = ATOB2_dat[3] | ( ( gJetTOBv[1] & tobvMask ) << tobvShift);
521  ATOB2_dat[3] = ATOB2_dat[3] | ( ( gJetTOBeta[1] & 0x0000003F ) <<20);
522  ATOB2_dat[3] = ATOB2_dat[3] | ( ( gJetTOBphi[1] & 0x0000001F ) <<26);
523 
524  tobs_v[11] = std::make_unique<gFEXJetTOB>();
525  tobs_v[11]->setWord(ATOB2_dat[3]);
526  tobs_v[11]->setET(gJetTOBv[1]);
527  tobs_v[11]->setEta(gJetTOBeta[1]);
528  tobs_v[11]->setPhi(gJetTOBphi[1]);
529  tobs_v[11]->setTobID(6);
530  if( gJetTOBv[1] > gLJ_ptMinToTopoCounts1 ) tobs_v[11]->setStatus(1);
531  else tobs_v[11]->setStatus(0);
532 
533 
534  BTOB1_dat[3] = 0x00000005;
535  if( gJetTOBv[2] > gLJ_ptMinToTopoCounts1 ) BTOB1_dat[3] = BTOB1_dat[3] | 0x00000080;
536  BTOB1_dat[3] = BTOB1_dat[3] | ( ( gJetTOBv[2] & tobvMask) << tobvShift);
537  BTOB1_dat[3] = BTOB1_dat[3] | ( ( gJetTOBeta[2] & 0x0000003F ) <<20);
538  BTOB1_dat[3] = BTOB1_dat[3] | ( ( gJetTOBphi[2] & 0x0000001F ) <<26);
539 
540  tobs_v[12] = std::make_unique<gFEXJetTOB>();
541  tobs_v[12]->setWord(BTOB1_dat[3]);
542  tobs_v[12]->setET(gJetTOBv[2]);
543  tobs_v[12]->setEta(gJetTOBeta[2]);
544  tobs_v[12]->setPhi(gJetTOBphi[2]);
545  tobs_v[12]->setTobID(5);
546  if( gJetTOBv[2] > gLJ_ptMinToTopoCounts1 ) tobs_v[12]->setStatus(1);
547  else tobs_v[12]->setStatus(0);
548 
549 
550  BTOB2_dat[3] = 0x00000006;
551  if( gJetTOBv[3] > gLJ_ptMinToTopoCounts2 ) BTOB2_dat[3] = BTOB2_dat[3] | 0x00000080;
552  BTOB2_dat[3] = BTOB2_dat[3] | ( ( gJetTOBv[3] & tobvMask ) << tobvShift);
553  BTOB2_dat[3] = BTOB2_dat[3] | ( ( gJetTOBeta[3] & 0x0000003F ) <<20);
554  BTOB2_dat[3] = BTOB2_dat[3] | ( ( gJetTOBphi[3] & 0x0000001F ) <<26);
555 
556  tobs_v[13] = std::make_unique<gFEXJetTOB>();
557  tobs_v[13]->setWord(BTOB2_dat[3]);
558  tobs_v[13]->setET(gJetTOBv[3]);
559  tobs_v[13]->setEta(gJetTOBeta[3]);
560  tobs_v[13]->setPhi(gJetTOBphi[3]);
561  tobs_v[13]->setTobID(6);
562  if( gJetTOBv[3] > gLJ_ptMinToTopoCounts2 ) tobs_v[13]->setStatus(1);
563  else tobs_v[13]->setStatus(0);
564 
565  // // CECILIA - TO BE CHECKED
566  CTOB1_dat[3] = 0x00000005;
567  if( gJetTOBv[4] > gLJ_ptMinToTopoCounts1 ) CTOB1_dat[3] = CTOB1_dat[3] | 0x00000080;
568  CTOB1_dat[3] = CTOB1_dat[3] | ( ( gJetTOBv[4] & tobvMask) << tobvShift);
569  CTOB1_dat[3] = CTOB1_dat[3] | ( ( gJetTOBeta[4] & 0x0000003F ) <<20);
570  CTOB1_dat[3] = CTOB1_dat[3] | ( ( gJetTOBphi[4] & 0x0000001F ) <<26);
571 
572  CTOB2_dat[3] = 0x00000006;
573  if( gJetTOBv[5] > gLJ_ptMinToTopoCounts2 ) CTOB2_dat[3] = CTOB2_dat[3] | 0x00000080;
574  CTOB2_dat[3] = CTOB2_dat[3] | ( ( gJetTOBv[5] & tobvMask ) << tobvShift);
575  CTOB2_dat[3] = CTOB2_dat[3] | ( ( gJetTOBeta[5] & 0x0000003F ) <<20);
576  CTOB2_dat[3] = CTOB2_dat[3] | ( ( gJetTOBphi[5] & 0x0000001F ) <<26);
577 
578 
579 // saturation currently not set in firmware
580 //
581 //
582  bool setSat = 0;
583  if ( setSat) {
584  if( gTOBsat[0] ) {
585  ATOB1_dat[1] = ( ATOB1_dat[1] | 0x80000000 );
586  ATOB1_dat[3] = ( ATOB1_dat[3] | 0x80000000 );
587  }
588 
589  if( gTOBsat[1] ) {
590  ATOB2_dat[1] = ( ATOB2_dat[1] | 0x80000000 );
591  ATOB2_dat[3] = ( ATOB2_dat[3] | 0x80000000 );
592  }
593 
594  if( gTOBsat[2] ) {
595  BTOB1_dat[1] = ( BTOB1_dat[1] | 0x80000000 );
596  BTOB1_dat[3] = ( BTOB1_dat[3] | 0x80000000 );
597  }
598 
599  if( gTOBsat[3] ) {
600  BTOB2_dat[1] = ( BTOB2_dat[1] | 0x80000000 );
601  BTOB2_dat[3] = ( BTOB2_dat[3] | 0x80000000 );
602  }
603 
604  if( gTOBsat[4] ) {
605  CTOB1_dat[1] = ( CTOB1_dat[1] | 0x80000000 );
606  CTOB1_dat[3] = ( CTOB1_dat[3] | 0x80000000 );
607  }
608 
609  if( gTOBsat[5] ) {
610  CTOB2_dat[1] = ( CTOB2_dat[1] | 0x80000000 );
611  CTOB2_dat[3] = ( CTOB2_dat[3] | 0x80000000 );
612  }
613  }
614 
615  // zero tob 4 word
616 
617  ATOB1_dat[4] = 0 ;
618  ATOB2_dat[4] = 0 ;
619  BTOB1_dat[4] = 0 ;
620  BTOB2_dat[4] = 0 ;
621  CTOB1_dat[4] = 0 ;
622  CTOB2_dat[4] = 0 ;
623 
624 
625  // add seven bits as in the hardware of BCID in 6th TOB word ( index 5)
626  int BCID = 0;
627 
628  ATOB1_dat[5] = ( (BCID&0x0000007F)<<8 ) ;
629  ATOB2_dat[5] = ( (BCID&0x0000007F)<<8 ) ;
630  BTOB1_dat[5] = ( (BCID&0x0000007F)<<8 ) ;
631  BTOB2_dat[5] = ( (BCID&0x0000007F)<<8 ) ;
632  CTOB1_dat[5] = ( (BCID&0x0000007F)<<8 ) ;
633  CTOB2_dat[5] = ( (BCID&0x0000007F)<<8 ) ;
634 
635  // 3 is gFEX FEX number, set CRC to zero for now
636  int CRC = 0;
637 
638 
639  ATOB1_dat[6] = 0x000000BC | ( (BCID&0x0000000F)<<8 ) | (3<<12) | (CRC<<23);
640  ATOB2_dat[6] = 0x000000BC | ( (BCID&0x0000000F)<<8 ) | (3<<12) | (CRC<<23);
641  BTOB1_dat[6] = 0x000000BC | ( (BCID&0x0000000F)<<8 ) | (3<<12) | (CRC<<23);
642  BTOB2_dat[6] = 0x000000BC | ( (BCID&0x0000000F)<<8 ) | (3<<12) | (CRC<<23);
643  CTOB1_dat[6] = 0x000000BC | ( (BCID&0x0000000F)<<8 ) | (3<<12) | (CRC<<23);
644  CTOB2_dat[6] = 0x000000BC | ( (BCID&0x0000000F)<<8 ) | (3<<12) | (CRC<<23);
645 
646 
647  return tobs_v;
648 
649 }
650 
651 
652 void gFEXJetAlgo::RemotePartialAB(const gTowersType& twrs, gTowersPartialSums & lps, gTowersPartialSums & rps) const {
653 
654  // Computes partial sums for FPGA A or B
655  // twrs are the 32 x 12 = 284 gTowers in FPGA A or B
656 
657  // number of rows above/below for right partial sum
658  // when sending to the right send values for largest partial sum first, i.e. for column 0
659  typedef std::array<std::array<int, 4>, 4> gTowersNeighbours;
660  gTowersNeighbours NUpDwnR = {{ {{2,3,4,4}}, {{0,2,3,4}}, {{0,0,2,3}}, {{0,0,0,2}} }};
661 
662  // number of rows above or below for left partial sum
663  // when sending to the left send values for smallest partial sum first (i.e. for column 8 )
664  gTowersNeighbours NUpDwnL = {{ {{2,0,0,0}}, {{3,2,0,0}}, {{4,3,2,0}}, {{4,4,3,2}} }};
665 
666  // Do partial sum for output to right FPGA first
667  for( int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
668  for(int rcolumn = 0; rcolumn<FEXAlgoSpaceDefs::n_partial; rcolumn++){
669  // start calculating right partial sums for remote column rcolumn
670  rps[irow][rcolumn] = 0;
671  for(int lcolumn = 0; lcolumn<FEXAlgoSpaceDefs::n_partial; lcolumn++){
672  // add in any energy from towers in this row
673  // no need of modular arithmetic
674  if ( NUpDwnR[rcolumn][lcolumn] > 0 ) {
675  rps[irow][rcolumn] = rps[irow][rcolumn] + twrs[irow][lcolumn+8];
676  }
677  // now add rup1, rup2, rup3, rup4, ldn1, ldn2, ln3, ln4 -- use a loop instead of enumeratin in firmware
678  for( int rowOff = 1 ; rowOff < NUpDwnR[rcolumn][lcolumn]+1; rowOff++){
679  int rowModUp = (irow + rowOff)%32;
680  int rowModDn = (irow - rowOff + 32 )%32;
681  // printf("row %d remote column %d local column %d rowOff %d rowModUp %d rowModDn %d \n",
682  rps[irow][rcolumn] = rps[irow][rcolumn] + twrs[rowModUp][lcolumn+8] + twrs[rowModDn][lcolumn+8];
683  }
684  }
685  }
686  }
687 
688  // do partial sum for output to left FPGA second
689  for( int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
690  for(int rcolumn = 0; rcolumn<FEXAlgoSpaceDefs::n_partial; rcolumn++){
691  // start calculating right partial sums for remote column rcolumn
692  lps[irow][rcolumn] = 0;
693  for(int lcolumn = 0; lcolumn<FEXAlgoSpaceDefs::n_partial; lcolumn++){
694  // add in any energy from towers in this row
695  // no need of modular arithmetic
696  if ( NUpDwnL[rcolumn][lcolumn] > 0 ) {
697  lps[irow][rcolumn] = lps[irow][rcolumn] + twrs[irow][lcolumn];
698  }
699  // now add rup1, rup2, rup3, rup4, ldn1, ldn2, ln3, ln4 -- use a loop instead of enumeratin in firmware
700  for( int rowOff = 1 ; rowOff < NUpDwnL[rcolumn][lcolumn]+1; rowOff++){
701  int rowModUp = (irow + rowOff)%32;
702  int rowModDn = (irow - rowOff + 32 )%32;
703  lps[irow][rcolumn] = lps[irow][rcolumn] + twrs[rowModUp][lcolumn] + twrs[rowModDn][lcolumn];
704  }
705  }
706  }
707  }
708 
709 }
710 
711 
712 void gFEXJetAlgo::RemotePartialCN(const gTowersJetEngine& twrs, gTowersPartialSums & rps ) const {
713 
714  typedef std::array<std::array<int, 4>, 4> gTowersNeighbours;
715  gTowersNeighbours NUpDwnR = {{ {{2,3,4,4}}, {{0,2,3,4}}, {{0,0,2,3}}, {{0,0,0,2}} }};
716 
717  // same as AB for right partial sum
718  // starts from the righ most partial sum -- largest partial sum
719  for( int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
720  for(int rcolumn = 0; rcolumn<FEXAlgoSpaceDefs::n_partial; rcolumn++){
721  rps[irow][rcolumn] = 0;
722  for(int lcolumn = 0; lcolumn<FEXAlgoSpaceDefs::n_partial; lcolumn++){
723  if ( NUpDwnR[rcolumn][lcolumn] > 0 ) {
724  rps[irow][rcolumn] = rps[irow][rcolumn] + twrs[irow][lcolumn+2];
725  }
726  for( int rowOff = 1 ; rowOff < NUpDwnR[rcolumn][lcolumn]+1; rowOff++){
727  int rowModUp = (irow + rowOff)%32;
728  int rowModDn = (irow - rowOff + 32 )%32;
729  rps[irow][rcolumn] = rps[irow][rcolumn] + twrs[rowModUp][lcolumn+2] + twrs[rowModDn][lcolumn+2];
730  }
731  }
732  }
733  }
734 }
735 
736 
737 void gFEXJetAlgo::RemotePartialCP(const gTowersJetEngine& twrs, gTowersPartialSums & lps ) const {
738 
739  typedef std::array<std::array<int, 4>, 4> gTowersNeighbours;
740  gTowersNeighbours NUpDwnL = {{ {{2,0,0,0}}, {{3,2,0,0}}, {{4,3,2,0}}, {{4,4,3,2}} }};
741  // do partial sum for output to left FPGA second
742  for( int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
743  for(int rcolumn = 0; rcolumn<FEXAlgoSpaceDefs::n_partial; rcolumn++){
744  // start calculating right partial sums for remote column rcolumn
745  lps[irow][rcolumn] = 0;
746  for(int lcolumn = 0; lcolumn<FEXAlgoSpaceDefs::n_partial; lcolumn++){
747  // add in any energy from towers in this row
748  // no need of modular arithmetic
749  if ( NUpDwnL[rcolumn][lcolumn] > 0 ) {
750  lps[irow][rcolumn] = lps[irow][rcolumn] + twrs[irow][lcolumn];
751  }
752  for( int rowOff = 1 ; rowOff < NUpDwnL[rcolumn][lcolumn]+1; rowOff++){
753  int rowModUp = (irow + rowOff)%32;
754  int rowModDn = (irow - rowOff + 32 )%32;
755  lps[irow][rcolumn] = lps[irow][rcolumn] + twrs[rowModUp][lcolumn] + twrs[rowModDn][lcolumn];
756  }
757  }
758  }
759  }
760 }
761 
762 
763 void gFEXJetAlgo::singleHalf(const gTowersType& twrs, gTowersType & FPGAsum) const {
764  // Finds jets in a single FPGA
765  // This version only gives sum in the right or left half of the FPGA
766 
767  // Number of up and down FPGAs to add
768  std::array<int, 9> NUpDwn = {{2,3,4,4,4,4,4,3,2}};
769 
770  // Evaluate jet centered on each possible eta and phi gTower
771  for( int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
772  for(int jcolumn = 0; jcolumn<FEXAlgoSpaceDefs::ABcolumns; jcolumn++){
773  // zero jet sum here
774  FPGAsum[irow][jcolumn] = 0;
775 
776  // local column goes from 0 to 8
777  for(int localColumn = 0; localColumn<9; localColumn++){
778  int sumColumn = jcolumn + localColumn - 4;
779 
780  // calculate min and max columns for FPGA half
781  int lmin = 0;
782  int lmax = 5;
783 
784  if( jcolumn > 5 ) {
785  lmin = 6;
786  lmax = 11;
787  }
788 
789  // check if tower to be summed is actually in FPGA
790  if( (sumColumn >= lmin) && (sumColumn <= lmax) ) {
791 
792  // sum tower in the same row as jet center
793  FPGAsum[irow][jcolumn] = FPGAsum[irow][jcolumn] + twrs[irow][sumColumn];
794 
795  // add rows above and below according NUpDwn
796  // localRow goes from 1 to 2, 1 to 3 or 1 to 4
797  for(int localRow = 1; localRow<=NUpDwn[localColumn]; localRow++){
798  int krowUp = (irow + localRow);
799  if( krowUp > 31 ) krowUp = krowUp - 32;
800  int krowDn = (irow - localRow);
801  if( krowDn < 0 ) krowDn = krowDn + 32;
802  FPGAsum[irow][jcolumn] = FPGAsum[irow][jcolumn] +
803  twrs[krowUp][sumColumn] +
804  twrs[krowDn][sumColumn];
805  }
806  }
807  }
808  }
809  }
810 }
811 
812 
813 void gFEXJetAlgo::InternalPartialAB(const gTowersType & twrs, gTowersPartialSums & lps, gTowersPartialSums & rps ) const{
814 
815  // Computes partial sums for FPGA A or B
816  // twrs are the 32 x 12 = 284 gTowers in FPGA A or B
817  //
818  // lps: partial sum for left hand side of FPGA
819  // rps: partial sum for right half of the FPGA
820 
821  // NOTE rps is available first
822 
823  // number of rows above/below for right partial sum
824  // when sending to the right send values for largest partial sum first, i.e. for column 6
825  unsigned int NUpDwnR[4][4] = { {2,3,4,4}, {0,2,3,4}, {0,0,2,3}, {0,0,0,2} };
826 
827  // number of rows above or below for left partial sum
828  // when sending to the left send values for smallest partial sumn first (i.e. for column 2 )
829  unsigned int NUpDwnL[4][4] = { {2,0,0,0}, {3,2,0,0}, {4,3,2,0}, {4,4,3,2} };
830 
831 
832  // Do partial sum for output to right side FPGA first
833  // for the right partial sum this starts from the right-most value first
834  for(unsigned int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
835  for(unsigned int rcolumn = 0; rcolumn<FEXAlgoSpaceDefs::n_partial; rcolumn++){
836  // start calculating right partial sums for remote column rcolumn
837  rps[irow][rcolumn] = 0;
838  for(unsigned int lcolumn = 0; lcolumn<FEXAlgoSpaceDefs::n_partial; lcolumn++){
839  // add in any energy from towers in this row
840  // no need of modular arithmetic
841  if ( NUpDwnR[rcolumn][lcolumn] > 0 ) {
842  // this is partial sum for the right half of the FPGA -- columns 2,3,4,5
843  rps[irow][rcolumn] = rps[irow][rcolumn] + twrs[irow][lcolumn+2];
844  }
845  // now add rup1, rup2, rup3, rup4, ldn1, ldn2, ln3, ln4 -- use a loop instead of enumeratin in firmware
846  for(unsigned int rowOff = 1 ; rowOff < NUpDwnR[rcolumn][lcolumn]+1; rowOff++){
847  int rowModUp = (irow + rowOff)%32;
848  int rowModDn = (irow - rowOff + 32 )%32;
849  // this is partial sum for the right half of the FPGA -- columns 2,3,4,5
850  rps[irow][rcolumn] = rps[irow][rcolumn] + twrs[rowModUp][lcolumn+2] + twrs[rowModDn][lcolumn+2];
851  }
852  }
853  }
854  }
855  // do partial sum for output to left half of the FPGA second
856  // for the left partial sum this starts also with the right most column -- which is the smallest sum
857  for(unsigned int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
858  for(unsigned int rcolumn = 0; rcolumn<FEXAlgoSpaceDefs::n_partial; rcolumn++){
859  // start calculating right partial sums for remote column rcolumn
860  lps[irow][rcolumn] = 0;
861  for(unsigned int lcolumn = 0; lcolumn<FEXAlgoSpaceDefs::n_partial; lcolumn++){
862  // add in any energy from towers in this row
863  // no need of modular arithmetic
864  if ( NUpDwnL[rcolumn][lcolumn] > 0 ) {
865  lps[irow][rcolumn] = lps[irow][rcolumn] + twrs[irow][lcolumn+6];
866  }
867  // now add rup1, rup2, rup3, rup4, ldn1, ldn2, ln3, ln4 -- use a loop instead of enumeratin in firmware
868  for(unsigned int rowOff = 1 ; rowOff < NUpDwnL[rcolumn][lcolumn]+1; rowOff++){
869  int rowModUp = (irow + rowOff)%32;
870  int rowModDn = (irow - rowOff + 32 )%32;
871  // this is partial sum for the left half of the FPGA -- columns 6,7,8,9
872  lps[irow][rcolumn] = lps[irow][rcolumn] + twrs[rowModUp][lcolumn+6] + twrs[rowModDn][lcolumn+6];
873  }
874  }
875  }
876  }
877 }
878 
879 
880 void gFEXJetAlgo::addInternalLin(gTowersType & jets, gTowersPartialSums & partial) const{
881  // add parial sums for left side of FPGA
882  for(unsigned int irow=0;irow<FEXAlgoSpaceDefs::ABCrows;irow++){
883  for(unsigned int ipartial= 0; ipartial <FEXAlgoSpaceDefs::n_partial; ipartial++){
884  jets[irow][2+ipartial] = jets[irow][2+ipartial] + partial[irow][ipartial];
885  }
886  }
887 }
888 
889 
890 void gFEXJetAlgo::addInternalRin(gTowersType & jets, gTowersPartialSums & partial) const{
891  // add parial sums for right side of FPGA
892  for(unsigned int irow=0;irow<FEXAlgoSpaceDefs::ABCrows;irow++){
893  for(unsigned int ipartial= 0; ipartial <FEXAlgoSpaceDefs::n_partial; ipartial++){
894  jets[irow][6+ipartial] = jets[irow][6+ipartial] + partial[irow][ipartial];
895  }
896  }
897 }
898 
899 
900 void gFEXJetAlgo::pileUpCorrectionAB(gTowersType &jets, int puc) const {
901  int rows = jets.size();
902  int cols = jets[0].size();
903  for(int irow=0;irow<rows;irow++){
904  for(int icolumn=0;icolumn<cols;icolumn++){
905  jets[irow][icolumn] = jets[irow][icolumn] - puc;
906  }
907  }
908 }
909 
910 
911 void gFEXJetAlgo::ZeroNegative(gTowersType & jets) const{
912  for(unsigned int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
913  for(unsigned int icolumn =0; icolumn<FEXAlgoSpaceDefs::ABcolumns; icolumn++){
914  if(jets[irow][icolumn] < 0) jets[irow][icolumn] = 0;
915  }
916  }
917 }
918 
919 
920 // https://gitlab.cern.ch/atlas-l1calo/gfex/firmware/-/blob/devel/common/jet_finder/HDL/jet_eng.vhd#L538
921 void gFEXJetAlgo::SaturateJets( gTowersType & jets, const gTowersType & sat ) const {
922  for(unsigned int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
923  for(unsigned int icolumn =0; icolumn<FEXAlgoSpaceDefs::ABcolumns; icolumn++){
924  if(static_cast<unsigned>(sat[irow][icolumn])) {
925  // this should perhaps be 0xfff -- firmware turncates to 12 bits before sorting.
926  jets[irow][icolumn] = 0x0003ffff;
927  }
928  }
929  }
930 }
931 
932 void gFEXJetAlgo::SaturateBlocks( gTowersType & gBlkSum, const gTowersType & sat ) const {
933  for(unsigned int irow = 0; irow < FEXAlgoSpaceDefs::ABCrows; irow++ ){
934  for(unsigned int icolumn =0; icolumn<FEXAlgoSpaceDefs::ABcolumns; icolumn++){
935  if(static_cast<unsigned>(sat[irow][icolumn])) {
936  gBlkSum[irow][icolumn] = 0x00000fff;
937  }
938  }
939  }
940 }
941 
942 
943 
944 void gFEXJetAlgo::gBlockAB(const gTowersType & twrs, gTowersType & gBlkSum, gTowersType & hasSeed, int seedThreshold) const {
945 
946  int rows = twrs.size();
947  int cols = twrs[0].size();
948  for( int irow = 0; irow < rows; irow++ ){
949  for(int jcolumn = 0; jcolumn<cols; jcolumn++){
950  // zero jet sum here
951  gBlkSum[irow][jcolumn] = 0;
952  int krowUp = (irow + 1)%32;
953  int krowDn = (irow - 1 +32)%32;
954  if( (jcolumn == 0) || (jcolumn == 6) ) {
955  //left edge case
956  gBlkSum[irow][jcolumn] =
957  twrs[irow][jcolumn] + twrs[krowUp][jcolumn] + twrs[krowDn][jcolumn] +
958  twrs[irow][jcolumn+1] + twrs[krowUp][jcolumn+1] + twrs[krowDn][jcolumn+1];
959  } else if( (jcolumn == 5) || (jcolumn == 11)) {
960  // right edge case
961  gBlkSum[irow][jcolumn] =
962  twrs[irow][jcolumn] + twrs[krowUp][jcolumn] + twrs[krowDn][jcolumn] +
963  twrs[irow][jcolumn-1] + twrs[krowUp][jcolumn-1] + twrs[krowDn][jcolumn-1];
964  } else{
965  // normal case
966  gBlkSum[irow][jcolumn] =
967  twrs[irow][jcolumn] + twrs[krowUp][jcolumn] + twrs[krowDn][jcolumn] +
968  twrs[irow][jcolumn-1] + twrs[krowUp][jcolumn-1] + twrs[krowDn][jcolumn-1] +
969  twrs[irow][jcolumn+1] + twrs[krowUp][jcolumn+1] + twrs[krowDn][jcolumn+1];
970  }
971 
972  if( gBlkSum[irow][jcolumn] > seedThreshold) {
973  hasSeed[irow][jcolumn] = 1;
974  } else {
975  hasSeed[irow][jcolumn] = 0;
976  }
977 
978  if ( gBlkSum[irow][jcolumn] < 0 )
979  gBlkSum[irow][jcolumn] = 0;
980 
981  // was bits 11+3 downto 3, now is 11 downto 0
982  if ( gBlkSum[irow][jcolumn] > FEXAlgoSpaceDefs::gBlockMax ) {
983  gBlkSum[irow][jcolumn] = FEXAlgoSpaceDefs::gBlockMax;
984  }
985  }
986  }
987 }
988 
989 
990 
991 void gFEXJetAlgo::blkOutAB(gTowersType & blocks, std::array<int, 32> & jetOutL, std::array<int, 32> & etaIndL, std::array<int, 32> & jetOutR, std::array<int, 32> & etaIndR) const{
992 
993  // find maximum in each jet engine for gBlocks (not done in hardware)
994 
995  //loop over left engines
996  for(unsigned int ieng=0; ieng<FEXAlgoSpaceDefs::ABCrows; ieng++){
997  for(unsigned int localEta = 0; localEta<6; localEta++){
998  if( blocks[ieng][localEta] > jetOutL[ieng] ){
999  jetOutL[ieng] = blocks[ieng][localEta];
1000  etaIndL[ieng] = localEta;
1001  }
1002  }
1003  }
1004  // loop over right engines
1005  for(unsigned int ieng=0; ieng<FEXAlgoSpaceDefs::ABCrows; ieng++){
1006  for(unsigned int localEta = 0; localEta<6; localEta++){
1007  if( blocks[ieng][localEta+6] > jetOutR[ieng] ) {
1008  jetOutR[ieng] = blocks[ieng][localEta+6];
1009  etaIndR[ieng] = localEta;
1010  }
1011  }
1012  }
1013 }
1014 
1015 
1016 void gFEXJetAlgo::gBlockMax2(const gTowersType & gBlkSum, int BjetColumn, int localColumn, std::array<int, 3> & gBlockV, std::array<int, 3> & gBlockEta, std::array<int, 3> & gBlockPhi) const {
1017  // gBlkSum are the 9 or 6 gTower sums
1018  // BjetColumn is the Block Column -- 0 for CN, 1, 2 for A 3, 4 for B and 5 for CP
1019  // gBlockV is the array of values currently 2
1020  // gBlockEta is the eta in global
1021 
1022  gTowersJetEngine gBlkSumC;
1023 
1024  // copy the correct column
1025  for( int irow = 0; irow<FEXAlgoSpaceDefs::ABCrows; irow++){
1026  for( int icolumn =0; icolumn<FEXAlgoSpaceDefs::ABCcolumnsEng; icolumn++){
1027  gBlkSumC[irow][icolumn] = gBlkSum[irow][icolumn + localColumn*FEXAlgoSpaceDefs::ABCcolumnsEng];
1028  }
1029  }
1030 
1031  gBlockMax192(gBlkSumC, gBlockV, gBlockEta, gBlockPhi, 0);
1032 
1033  for(int i = -1; i<2; i++){
1034  int iGlobal = i + gBlockPhi[0];
1035  // map row number in to 0-31
1036  if ( iGlobal < 0 ) iGlobal = iGlobal + 32;
1037  if ( iGlobal > 31 ) iGlobal = iGlobal - 32;
1038  // don't do anything outside of the six columsn
1039  for( int j = -1; j<2; j++){
1040  int jGlobal = j + gBlockEta[0];
1041  if( (jGlobal > -1) && (jGlobal < 6) ) {
1042  gBlkSumC[iGlobal][jGlobal] = 0;
1043  }
1044  }
1045  }
1046 
1047  gBlockMax192(gBlkSumC, gBlockV, gBlockEta, gBlockPhi, 1);
1048 
1049  // need to wait until the end for this!
1050  gBlockEta[0] = gBlockEta[0] + 2 + 6*BjetColumn;
1051  gBlockEta[1] = gBlockEta[1] + 2 + 6*BjetColumn;
1052 }
1053 
1054 
1055 void gFEXJetAlgo::gBlockMax192(const gTowersJetEngine& gBlkSum,
1056  std::array<int, 3> & gBlockVp,
1057  std::array<int, 3> & gBlockEtap,
1058  std::array<int, 3> & gBlockPhip,
1059  int index) const {
1060  // gBLKSum are the sums
1061  // gBlockVp is the returned value for the max block
1062  // gBlockEtap is the eta in the local corrdinate system
1063  // gBlockPhip is the phi (global and local are the same)
1064 
1065  int inpv[192]{};
1066  int maxv1[96]{};
1067  int maxv2[48]{};
1068  int maxv3[24]{};
1069  int maxv4[12]{};
1070  int maxv5[6]{};
1071  int maxv6[3]{};
1072 
1073  int inpi[192]{};
1074  int maxi1[96]{};
1075  int maxi2[48]{};
1076  int maxi3[24]{};
1077  int maxi4[12]{};
1078  int maxi5[6]{};
1079  int maxi6[3]{};
1080 
1081  int maxvall;
1082  int maxiall;
1083 
1084 
1085  // ABCcolumnsEng is 6 in hardware
1086  int maxv = 0;
1087 
1088  for(unsigned int icolumn = 0; icolumn<FEXAlgoSpaceDefs::ABCcolumnsEng; icolumn++){
1089  for(unsigned int irow = 0; irow<FEXAlgoSpaceDefs::ABCrows; irow++){
1090  inpv[irow + icolumn*FEXAlgoSpaceDefs::ABCrows] = gBlkSum[irow][icolumn];
1091  inpi[irow + icolumn*FEXAlgoSpaceDefs::ABCrows] = irow + icolumn*FEXAlgoSpaceDefs::ABCrows;
1092 
1093  if( gBlkSum[irow][icolumn] > maxv){
1094  maxv = gBlkSum[irow][icolumn];
1095  }
1096  }
1097  }
1098 
1099  for(int i = 0; i<96; i++){
1100  if( inpv[2*i+1] > inpv[2*i] ) {
1101  maxv1[i] = inpv[2*i+1];
1102  maxi1[i] = inpi[2*i+1];
1103  } else {
1104  maxv1[i] = inpv[2*i];
1105  maxi1[i] = inpi[2*i];
1106  }
1107  }
1108 
1109  for(int i = 0; i<48; i++){
1110  if( maxv1[2*i+1] > maxv1[2*i] ) {
1111  maxv2[i] = maxv1[2*i+1];
1112  maxi2[i] = maxi1[2*i+1];
1113  } else {
1114  maxv2[i] = maxv1[2*i];
1115  maxi2[i] = maxi1[2*i];
1116  }
1117  }
1118 
1119  for(int i = 0; i<24; i++){
1120  if( maxv2[2*i+1] > maxv2[2*i] ) {
1121  maxv3[i] = maxv2[2*i+1];
1122  maxi3[i] = maxi2[2*i+1];
1123  } else {
1124  maxv3[i] = maxv2[2*i];
1125  maxi3[i] = maxi2[2*i];
1126  }
1127  }
1128 
1129  for(int i = 0; i<12; i++){
1130  if( maxv3[2*i+1] > maxv3[2*i] ) {
1131  maxv4[i] = maxv3[2*i+1];
1132  maxi4[i] = maxi3[2*i+1];
1133  } else {
1134  maxv4[i] = maxv3[2*i];
1135  maxi4[i] = maxi3[2*i];
1136  }
1137  }
1138 
1139  for(int i = 0; i<6; i++){
1140  if( maxv4[2*i+1] > maxv4[2*i] ) {
1141  maxv5[i] = maxv4[2*i+1];
1142  maxi5[i] = maxi4[2*i+1];
1143  } else {
1144  maxv5[i] = maxv4[2*i];
1145  maxi5[i] = maxi4[2*i];
1146  }
1147  }
1148 
1149  for(int i = 0; i<3; i++){
1150  if( maxv5[2*i+1] > maxv5[2*i] ) {
1151  maxv6[i] = maxv5[2*i+1];
1152  maxi6[i] = maxi5[2*i+1];
1153  } else {
1154  maxv6[i] = maxv5[2*i];
1155  maxi6[i] = maxi5[2*i];
1156  }
1157  }
1158 
1159  if( ( maxv6[1] > maxv6[0] ) && ( maxv6[1] > maxv6[2] ) ) {
1160  maxvall = maxv6[1];
1161  maxiall = maxi6[1];
1162  } else {
1163  if( maxv6[0] > maxv6[2] ){
1164  maxvall = maxv6[0];
1165  maxiall = maxi6[0];
1166  } else {
1167  maxvall = maxv6[2];
1168  maxiall = maxi6[2];
1169  }
1170  }
1171 
1172 
1173  gBlockVp[index] = maxvall ;
1174  gBlockEtap[index] = maxiall/32;
1175  gBlockPhip[index] = maxiall%32;
1176 
1177 }
1178 
1179 
1180 void gFEXJetAlgo::gBlockVetoAB(gTowersType &jets, gTowersType& hasSeed) const {
1181  int rows = jets.size();
1182  int cols = jets[0].size();
1183  for( int irow = 0; irow < rows; irow++ ){
1184  for(int jcolumn = 0; jcolumn < cols; jcolumn++){
1185  if( hasSeed[irow][jcolumn] == 0 ) {
1186  // set to negative value as in hardware
1187  // https://gitlab.cern.ch/atlas-l1calo/gfex/firmware/-/blob/devel/common/jet_finder/HDL/jet_eng.vhd#L561
1188  jets[irow][jcolumn] = 0XFFFFF000;
1189  }
1190  }
1191  }
1192 }
1193 
1194 
1195 void gFEXJetAlgo::addRemoteRin(gTowersType &jets, const gTowersPartialSums &partial,
1196  int ps_upper, int ps_lower, int ps_shift) const {
1197 
1198  // See https://gitlab.cern.ch/atlas-l1calo/gfex/firmware/-/blob/kw-dev/jwj_verif/common/jet_finder/HDL/ab_encode.vhd
1199 
1200  int rows = partial.size();
1201  int cols = partial[0].size();
1202  // add partial sums
1203  for(int irow=0;irow<rows;irow++){
1204  for(int ipartial= 0; ipartial <cols; ipartial++){
1205  // current behavior
1206  int truncPart = -1;
1207  if( partial[irow][ipartial] > ps_upper ) {
1208  truncPart = ps_upper;
1209  } else if ( partial[irow][ipartial] < ps_lower ) {
1210  truncPart = ps_lower;
1211  } else {
1212  truncPart = partial[irow][ipartial];
1213  }
1214  // change LSB from 200 MeV to 1600 MeV and then back to 200 MeV.
1215  truncPart = (truncPart >> ps_shift );
1216  truncPart = (truncPart << ps_shift );
1217 
1218  jets[irow][8+ipartial] = jets[irow][8+ipartial] + truncPart;
1219  }
1220  }
1221 }
1222 
1223 void gFEXJetAlgo::addRemoteLin(gTowersType &jets, const gTowersPartialSums &partial,
1224  int ps_upper, int ps_lower, int ps_shift) const {
1225 
1226  int rows = partial.size();
1227  int cols = partial[0].size();
1228  // add partial sums
1229  for(int irow=0;irow<rows;irow++){
1230  for(int ipartial= 0; ipartial <cols; ipartial++){
1231  // current behavior
1232  int truncPart = -1;
1233  if( partial[irow][ipartial] > ps_upper ) {
1234  truncPart = ps_upper;
1235  } else if ( partial[irow][ipartial] < ps_lower ) {
1236  truncPart = ps_lower;
1237  } else {
1238  truncPart = partial[irow][ipartial];
1239  }
1240  // change LSB from 200 MeV to 1600 MeV and then back to 200 MeV.
1241  truncPart = (truncPart >> ps_shift );
1242  truncPart = (truncPart << ps_shift );
1243 
1244  jets[irow][ipartial] = jets[irow][ipartial] + truncPart;
1245  }
1246  }
1247 }
1248 
1249 
1250 void gFEXJetAlgo::addRemoteCNin(gTowersType & jets, const gTowersPartialSums & partial,
1251  int ps_upper, int ps_lower, int ps_shift ) const {
1252 
1253  int rows = partial.size();
1254  int cols = partial[0].size();
1255  // add partial sums
1256  for(int irow=0; irow < rows; irow++){
1257  for(int ipartial= 0; ipartial <cols; ipartial++){
1258  // current behavior
1259  int truncPart = -1;
1260  if( partial[irow][ipartial] > ps_upper ) {
1261  truncPart = ps_upper;
1262  } else if ( partial[irow][ipartial] < ps_lower ) {
1263  truncPart = ps_lower;
1264  } else {
1265  truncPart = partial[irow][ipartial];
1266  }
1267 
1268  truncPart = (truncPart >> ps_shift );
1269  truncPart = (truncPart << ps_shift );
1270 
1271  jets[irow][ipartial+2] = jets[irow][ipartial+2] + truncPart;
1272 
1273  }
1274  }
1275 }
1276 
1277 void gFEXJetAlgo::addRemoteCPin(gTowersType & jets, const gTowersPartialSums & partial,
1278  int ps_upper, int ps_lower, int ps_shift ) const {
1279 
1280  int rows = partial.size();
1281  int cols = partial[0].size();
1282  // add partial sums
1283  for(int irow=0;irow<rows;irow++){
1284  for(int ipartial= 0; ipartial <cols; ipartial++){
1285  int truncPart = -1;
1286  if( partial[irow][ipartial] > ps_upper ) {
1287  truncPart = ps_upper;
1288  } else if ( partial[irow][ipartial] < ps_lower ) {
1289  truncPart = ps_lower;
1290  } else {
1291  truncPart = partial[irow][ipartial];
1292  }
1293 
1294  truncPart = (truncPart >> ps_shift );
1295  truncPart = (truncPart << ps_shift );
1296 
1297  jets[irow][ipartial+6] = jets[irow][ipartial+6] + truncPart;
1298 
1299  }
1300  }
1301 }
1302 
1303 
1304 void gFEXJetAlgo::gJetVetoAB( gTowersType &twrs ,int jet_threshold ) const {
1305  int rows = twrs.size();
1306  int cols = twrs[0].size();
1307  for( int irow = 0; irow < rows; irow++ ){
1308  for(int jcolumn = 0; jcolumn<cols; jcolumn++){
1309  if( twrs[irow][jcolumn] < jet_threshold+1 ) {
1310  twrs[irow][jcolumn] = 0;
1311  }
1312  }
1313  }
1314 }
1315 
1316 
1317 void gFEXJetAlgo::jetOutAB(const gTowersType & jets,
1318  std::array<int, 32> & jetOutL, std::array<int, 32> & etaIndL,
1319  std::array<int, 32> & jetOutR, std::array<int, 32> & etaIndR ) const {
1320  // find maximum in each jet engine or either gJets and requires corresponding gBlock be above threhsold
1321  //loop over left engines
1322  for(int ieng=0; ieng<FEXAlgoSpaceDefs::ABCrows; ieng++){
1323  jetOutL[ieng] = 0;
1324  etaIndL[ieng] = 0;
1325  for(int localEta = 0; localEta<6; localEta++){
1326  if( jets[ieng][localEta] > jetOutL[ieng] ){
1327  jetOutL[ieng] = jets[ieng][localEta];
1328  etaIndL[ieng] = localEta;
1329  }
1330  }
1331  // truncation and mapping from 12 to 15 bits
1332  // https://gitlab.cern.ch/atlas-l1calo/gfex/firmware/-/blob/kw-dev/jet_finder_abc/common/jet_finder/HDL/jet_eng.vhd#L561
1333  // https://gitlab.cern.ch/atlas-l1calo/gfex/firmware/-/blob/kw-dev/jet_finder_abc/common/jet_finder/HDL/jet_eng.vhd#L561
1334  // Turncate to 15 bits as in VHDL (Corresponds to 13 TeV)
1335 
1336  if( jetOutL[ieng] > FEXAlgoSpaceDefs::gJetMax ) jetOutL[ieng] = FEXAlgoSpaceDefs::gJetMax;
1337  if( jetOutL[ieng] < 0 ) jetOutL[ieng] = 0;
1338  }
1339  // loop over right engines
1340  for(int ieng=0; ieng<FEXAlgoSpaceDefs::ABCrows; ieng++){
1341  jetOutR[ieng] = 0;
1342  etaIndR[ieng] = 0;
1343  for(int localEta = 0; localEta < 6; localEta++){
1344  if( jets[ieng][localEta+6] > jetOutR[ieng] ){
1345  jetOutR[ieng] = jets[ieng][localEta+6];
1346  etaIndR[ieng] = localEta;
1347  }
1348  }
1349  // Turncate to 15 bits as in VHDL (orresponds to 13 TeV)
1350  if( jetOutR[ieng] > FEXAlgoSpaceDefs::gJetMax ) jetOutR[ieng] = FEXAlgoSpaceDefs::gJetMax;
1351  if( jetOutR[ieng] < 0 ) jetOutR[ieng] = 0;
1352  // https://gitlab.cern.ch/atlas-l1calo/gfex/firmware/-/blob/kw-dev/jet_finder_abc/common/jet_finder/HDL/jet_finder_abc.vhd#L1103
1353  // https://gitlab.cern.ch/atlas-l1calo/gfex/firmware/-/blob/kw-dev/jet_finder_abc/common/jet_finder/HDL/jet_finder_abc.vhd#L1119
1354 
1355  }
1356 
1357 }
1358 
1359 
1360 void gFEXJetAlgo::pileUpCalculation(gTowersType &twrs, int rhoThreshold_Max, int inputScale, int &PUCp /*, int &PUChres*/) const {
1361  // input are 50 MeV "fine" scale towers (i.e. inputScale = 1)
1362  // to use 200 MeV towers use inputScale = 4
1363  // PUCp output is the pileup correction for 69 towers at 200 MeV energy scale
1364 
1365  int rows = twrs.size();
1366  int cols = twrs[0].size();
1367  int pucSum = 0;
1368  int nSum = 0;
1369  for(int irow=0; irow<rows; irow++){
1370  for( int icolumn=0; icolumn<cols; icolumn++){
1371  int fineGT = twrs[irow][icolumn]*inputScale;
1372  // set floor and ceiling here (currently 255 and -256 )
1373  if (fineGT > FEXAlgoSpaceDefs::fineCeiling ) fineGT = FEXAlgoSpaceDefs::fineCeiling;
1374  if (fineGT < FEXAlgoSpaceDefs::fineFloor ) fineGT = FEXAlgoSpaceDefs::fineFloor;
1375 
1376  if( fineGT < rhoThreshold_Max ) {
1377  pucSum = pucSum + fineGT;
1378  nSum = nSum + 1;
1379  }
1380  }
1381  }
1382  // in firmware this is done with a lookup table see https://gitlab.cern.ch/atlas-l1calo/gfex/firmware/-/blob/devel/common/jet_finder/HDL/inv_lut19.vhd
1383  // See also https://gitlab.cern.ch/atlas-l1calo/gfex/firmware/-/blob/devel/common/jet_finder/HDL/gt_build_all_AB.vhd#L1471
1384 
1385  // oneOverN is stored as a 32 bit number in inv_lut19
1386  unsigned int oneOverNTab = 69<<25;
1387  if( nSum > 0 && nSum < 385 ) {
1388  oneOverNTab = FEXAlgoSpaceDefs::inv19[nSum];
1389  } else {
1390  oneOverNTab = 0;
1391  }
1392 
1393  int oneOverN = oneOverNTab ;
1394 
1395  // largest value should be 255*2^14 ~ 2^22 -- should easily fit in int -- try expliciting putting in int here.
1396  pucSum = pucSum * oneOverN;
1397  // PUChres = pucSum;
1398 
1399  // current system 19 bits in table (69*4096 is max value) sign exten
1400  pucSum = ( pucSum >> 14);
1401  PUCp = pucSum;
1402 
1403 }
1404 
1405 void gFEXJetAlgo::gJetTOBgen(const std::array<int, FEXAlgoSpaceDefs::ABCrows>& jetOut,
1406  const std::array<int, FEXAlgoSpaceDefs::ABCrows>& etaInd,
1407  int TOBnum, int jetThreshold, std::array<int, FEXAlgoSpaceDefs::gJetTOBfib> & gJetTOBs,
1408  std::array<int, FEXAlgoSpaceDefs::gJetTOBfib> & gJetTOBv,
1409  std::array<int, FEXAlgoSpaceDefs::gJetTOBfib> & gJetTOBeta,
1410  std::array<int, FEXAlgoSpaceDefs::gJetTOBfib> & gJetTOBphi ) const {
1411 
1412  int jetOutZS[FEXAlgoSpaceDefs::ABCrows]{};
1413  // apply the tobthreshold to the values
1414  //note that jetThreshold is not a configurable parameter in firmware, it is used to check that jet values are positive
1415  for( int irow =0; irow<FEXAlgoSpaceDefs::ABCrows; irow++){
1416  if( jetOut[irow] > jetThreshold ) {
1417  jetOutZS[irow] = jetOut[irow];
1418  } else {
1419  jetOutZS[irow] = 0 ;
1420  }
1421  }
1422 
1423 
1424  // offset of TOBs according to official format
1425  int etaOff[FEXAlgoSpaceDefs::gJetTOBfib] = {8,14,20,26,2,32};
1426 
1427  // see tob_gen.vhd
1428  int l1mv[16]{};
1429  int l1met[16]{};
1430  int l1mphi[16]{};
1431 
1432  int l2mv[8]{};
1433  int l2met[8]{};
1434  int l2mphi[8]{};
1435 
1436  int l3mv[4]{};
1437  int l3met[4]{};
1438  int l3mphi[4]{};
1439 
1440  int l4mv[2]{};
1441  int l4met[2]{};
1442  int l4mphi[2]{};
1443 
1444  for(int i=0; i<16; i++){
1445  if( jetOut[2*i + 1] > jetOutZS[2*i] ){
1446  l1mv[i] = jetOutZS[2*i + 1];
1447  l1met[i] = etaInd[2*i + 1];
1448  l1mphi[i] = 2*i + 1;
1449  } else {
1450  l1mv[i] = jetOutZS[2*i ];
1451  l1met[i] = etaInd[2*i ];
1452  l1mphi[i] = 2*i ;
1453  }
1454  }
1455 
1456  for(int i=0; i<8; i++){
1457  if( l1mv[2*i + 1] > l1mv[2*i] ){
1458  l2mv[i] = l1mv[ 2*i + 1 ];
1459  l2met[i] = l1met[ 2*i + 1];
1460  l2mphi[i] = l1mphi[2*i + 1];
1461  } else {
1462  l2mv[i] = l1mv[ 2*i ];
1463  l2met[i] = l1met[ 2*i ];
1464  l2mphi[i] = l1mphi[2*i ];
1465  }
1466  }
1467 
1468  for(int i=0; i<4; i++){
1469  if( l2mv[2*i + 1] > l2mv[2*i] ){
1470  l3mv[i] = l2mv[ 2*i + 1];
1471  l3met[i] = l2met[ 2*i + 1];
1472  l3mphi[i] = l2mphi[2*i + 1];
1473  } else {
1474  l3mv[i] = l2mv[2*i ];
1475  l3met[i] = l2met[2*i ];
1476  l3mphi[i] = l2mphi[2*i ];
1477  }
1478  }
1479 
1480  for(int i=0; i<2; i++){
1481  if( l3mv[2*i + 1] > l3mv[2*i] ){
1482  l4mv[i] = l3mv[ 2*i + 1];
1483  l4met[i] = l3met[ 2*i + 1];
1484  l4mphi[i] = l3mphi[2*i + 1];
1485  } else {
1486  l4mv[i] = l3mv[2*i ];
1487  l4met[i] = l3met[2*i ];
1488  l4mphi[i] = l3mphi[2*i ];
1489  }
1490  }
1491 
1492  if( l4mv[1] > l4mv[0] ){
1493  if(l4mv[1] > jetThreshold){
1494  gJetTOBv[TOBnum] = l4mv[1];
1495  gJetTOBeta[TOBnum] = l4met[1] + etaOff[TOBnum];
1496  gJetTOBphi[TOBnum] = l4mphi[1];
1497  gJetTOBs[TOBnum] = 1;
1498  } else {
1499  gJetTOBv[TOBnum] = l4mv[1];
1500  gJetTOBeta[TOBnum] = l4met[1] + etaOff[TOBnum];
1501  gJetTOBphi[TOBnum] = l4mphi[1];
1502  gJetTOBs[TOBnum] = 0;
1503  }
1504  } else {
1505  if(l4mv[0] > jetThreshold){
1506  gJetTOBv[TOBnum] = l4mv[0];
1507  gJetTOBeta[TOBnum] = l4met[0] + etaOff[TOBnum];
1508  gJetTOBphi[TOBnum] = l4mphi[0];
1509  gJetTOBs[TOBnum] = 1;
1510  } else {
1511  gJetTOBv[TOBnum] = l4mv[0];
1512  gJetTOBeta[TOBnum] = l4met[0] + etaOff[TOBnum];
1513  gJetTOBphi[TOBnum] = l4mphi[0];
1514  gJetTOBs[TOBnum] = 0;
1515  }
1516  }
1517 
1518 }
1519 
1520 
1521 } // namespace LVL1
LVL1::gFEXJetAlgo::ZeroNegative
virtual void ZeroNegative(gTowersType &jets) const
Definition: gFEXJetAlgo.cxx:911
LVL1::FEXAlgoSpaceDefs::ENABLE_INTER_ABC
constexpr static bool ENABLE_INTER_ABC
Definition: FEXAlgoSpaceDefs.h:164
LVL1::FEXAlgoSpaceDefs::inv19
constexpr static std::array< unsigned int, 385 > inv19
Definition: FEXAlgoSpaceDefs.h:182
LVL1::gFEXJetAlgo::pileUpCorrectionAB
virtual void pileUpCorrectionAB(gTowersType &jets, int puc) const
Definition: gFEXJetAlgo.cxx:900
gFEXJetTOB.h
LVL1::FEXAlgoSpaceDefs::PS_LOWER_C
constexpr static int PS_LOWER_C
Definition: FEXAlgoSpaceDefs.h:176
LVL1::gFEXJetAlgo::InternalPartialAB
virtual void InternalPartialAB(const gTowersType &twrs, gTowersPartialSums &lps, gTowersPartialSums &rps) const
Definition: gFEXJetAlgo.cxx:813
index
Definition: index.py:1
LVL1::gFEXJetAlgo::singleHalf
virtual void singleHalf(const gTowersType &twrs, gTowersType &FPGAsum) const
Definition: gFEXJetAlgo.cxx:763
LVL1::FEXAlgoSpaceDefs::ENABLE_INTER_AB
constexpr static bool ENABLE_INTER_AB
Definition: FEXAlgoSpaceDefs.h:162
gTowerContainer.h
LVL1::gFEXJetAlgo::gBlockMax2
virtual void gBlockMax2(const gTowersType &gBlkSum, int BjetColumn, int localColumn, std::array< int, 3 > &gBlockV, std::array< int, 3 > &gBlockEta, std::array< int, 3 > &gBlockPhi) const
Definition: gFEXJetAlgo.cxx:1016
LVL1::FEXAlgoSpaceDefs::PS_UPPER_AB
constexpr static int PS_UPPER_AB
Definition: FEXAlgoSpaceDefs.h:172
LVL1::gFEXJetAlgo::SaturateJets
virtual void SaturateJets(gTowersType &jets, const gTowersType &sat) const
Definition: gFEXJetAlgo.cxx:921
LVL1::FEXAlgoSpaceDefs::ABCcolumnsEng
constexpr static int ABCcolumnsEng
Definition: FEXAlgoSpaceDefs.h:157
LVL1::gFEXJetAlgo::RemotePartialAB
virtual void RemotePartialAB(const gTowersType &twrs, gTowersPartialSums &lps, gTowersPartialSums &rps) const
Definition: gFEXJetAlgo.cxx:652
LVL1::gFEXJetAlgo::RemotePartialCP
virtual void RemotePartialCP(const gTowersJetEngine &twrs, gTowersPartialSums &lps) const
Definition: gFEXJetAlgo.cxx:737
LVL1::gFEXJetAlgo::addRemoteLin
virtual void addRemoteLin(gTowersType &jets, const gTowersPartialSums &partial, int ps_upper, int ps_lower, int ps_shift) const
Definition: gFEXJetAlgo.cxx:1223
LVL1::gFEXJetAlgo::addRemoteCNin
virtual void addRemoteCNin(gTowersType &jets, const gTowersPartialSums &partial, int ps_upper, int ps_lower, int ps_shift) const
Definition: gFEXJetAlgo.cxx:1250
LVL1
eFexTowerBuilder creates xAOD::eFexTowerContainer from supercells (LATOME) and triggerTowers (TREX) i...
Definition: ICMMCPHitsCnvTool.h:18
LVL1::gFEXJetAlgo::gJetVetoAB
virtual void gJetVetoAB(gTowersType &twrs, int jet_threshold) const
Definition: gFEXJetAlgo.cxx:1304
LVL1::FEXAlgoSpaceDefs::gJetMax
constexpr static int gJetMax
Definition: FEXAlgoSpaceDefs.h:158
LVL1::gFEXJetAlgo::pileUpCalculation
virtual void pileUpCalculation(gTowersType &twrs, int rhoThreshold_Max, int inputScale, int &PUCp) const override
Definition: gFEXJetAlgo.cxx:1360
PixelByteStreamErrors::BCID
@ BCID
Definition: PixelByteStreamErrors.h:13
LVL1::gFEXJetAlgo::addInternalRin
virtual void addInternalRin(gTowersType &jets, gTowersPartialSums &partial) const
Definition: gFEXJetAlgo.cxx:890
LVL1::FEXAlgoSpaceDefs::ENABLE_INTER_C
constexpr static bool ENABLE_INTER_C
Definition: FEXAlgoSpaceDefs.h:163
gFEXJetAlgo.h
beamspotnt.cols
list cols
Definition: bin/beamspotnt.py:1114
LVL1::gTowersPartialSums
std::array< std::array< int, 4 >, 32 > gTowersPartialSums
Definition: gFEXJetAlgo.h:29
LVL1::FEXAlgoSpaceDefs::PS_UPPER_C
constexpr static int PS_UPPER_C
Definition: FEXAlgoSpaceDefs.h:175
LVL1::gFEXJetAlgo::initialize
virtual StatusCode initialize() override
standard Athena-Algorithm method
Definition: gFEXJetAlgo.cxx:28
LVL1::FEXAlgoSpaceDefs::gJetTOBfib
constexpr static int gJetTOBfib
Definition: FEXAlgoSpaceDefs.h:165
lumiFormat.i
int i
Definition: lumiFormat.py:85
LVL1::gFEXJetAlgo::m_gFEXJetAlgo_gTowerContainerKey
SG::ReadHandleKey< LVL1::gTowerContainer > m_gFEXJetAlgo_gTowerContainerKey
Definition: gFEXJetAlgo.h:112
LVL1::gFEXJetAlgo::gFEXJetAlgo
gFEXJetAlgo(const std::string &type, const std::string &name, const IInterface *parent)
Constructors.
Definition: gFEXJetAlgo.cxx:21
LVL1::gFEXJetAlgo::gJetTOBgen
virtual void gJetTOBgen(const std::array< int, FEXAlgoSpaceDefs::ABCrows > &jetOut, const std::array< int, FEXAlgoSpaceDefs::ABCrows > &etaInd, int TOBnum, int jetThreshold, std::array< int, FEXAlgoSpaceDefs::gJetTOBfib > &gJetTOBs, std::array< int, FEXAlgoSpaceDefs::gJetTOBfib > &gJetTOBv, std::array< int, FEXAlgoSpaceDefs::gJetTOBfib > &gJetTOBeta, std::array< int, FEXAlgoSpaceDefs::gJetTOBfib > &gJetTOBphi) const
Definition: gFEXJetAlgo.cxx:1405
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
LVL1::FEXAlgoSpaceDefs::BTOBFIB
constexpr static int BTOBFIB
Definition: FEXAlgoSpaceDefs.h:166
test_pyathena.parent
parent
Definition: test_pyathena.py:15
LVL1::gFEXJetAlgo::blkOutAB
virtual void blkOutAB(gTowersType &blocks, std::array< int, 32 > &jetOutL, std::array< int, 32 > &etaIndL, std::array< int, 32 > &jetOutR, std::array< int, 32 > &etaIndR) const
Definition: gFEXJetAlgo.cxx:991
LVL1::gFEXJetAlgo::jetOutAB
virtual void jetOutAB(const gTowersType &jets, std::array< int, 32 > &jetOutL, std::array< int, 32 > &etaIndL, std::array< int, 32 > &jetOutR, std::array< int, 32 > &etaIndR) const
Definition: gFEXJetAlgo.cxx:1317
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
beamspotnt.rows
list rows
Definition: bin/beamspotnt.py:1112
LVL1::gFEXJetAlgo::gBlockAB
virtual void gBlockAB(const gTowersType &twrs, gTowersType &gBlkSum, gTowersType &hasSeed, int seedThreshold) const
Definition: gFEXJetAlgo.cxx:944
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
TrigConf::name
Definition: HLTChainList.h:35
LVL1::FEXAlgoSpaceDefs::fineCeiling
constexpr static int fineCeiling
Definition: FEXAlgoSpaceDefs.h:179
LVL1::FEXAlgoSpaceDefs::ABcolumns
constexpr static int ABcolumns
Definition: FEXAlgoSpaceDefs.h:156
LVL1::FEXAlgoSpaceDefs::PS_SHIFT_C
constexpr static int PS_SHIFT_C
Definition: FEXAlgoSpaceDefs.h:177
LVL1::FEXAlgoSpaceDefs::gBlockMax
constexpr static int gBlockMax
Definition: FEXAlgoSpaceDefs.h:159
LVL1::FEXAlgoSpaceDefs::ABCrows
constexpr static int ABCrows
Definition: FEXAlgoSpaceDefs.h:155
LVL1::gFEXJetAlgo::gBlockVetoAB
virtual void gBlockVetoAB(gTowersType &jets, gTowersType &hasSeed) const
Definition: gFEXJetAlgo.cxx:1180
LVL1::gFEXJetAlgo::largeRfinder
virtual std::vector< std::unique_ptr< gFEXJetTOB > > largeRfinder(const gTowersType &Atwr, const gTowersType &Btwr, const gTowersType &CNtwr, const gTowersType &Asat, const gTowersType &Bsat, const gTowersType &CNsat, int pucA, int pucB, int pucC, int gLJ_seedThrA, int gLJ_seedThrB, int gLJ_seedThrC, int gJ_ptMinToTopoCounts1, int gJ_ptMinToTopoCounts2, int jetThreshold, int gLJ_ptMinToTopoCounts1, int gLJ_ptMinToTopoCounts2, std::array< uint32_t, 7 > &ATOB1_dat, std::array< uint32_t, 7 > &ATOB2_dat, std::array< uint32_t, 7 > &BTOB1_dat, std::array< uint32_t, 7 > &BTOB2_dat, std::array< uint32_t, 7 > &CTOB1_dat, std::array< uint32_t, 7 > &CTOB2_dat) const override
Definition: gFEXJetAlgo.cxx:38
LVL1::gFEXJetAlgo::SaturateBlocks
virtual void SaturateBlocks(gTowersType &gBlkSum, const gTowersType &sat) const
Definition: gFEXJetAlgo.cxx:932
LVL1::FEXAlgoSpaceDefs::PS_SHIFT_AB
constexpr static int PS_SHIFT_AB
Definition: FEXAlgoSpaceDefs.h:174
DeMoScan.index
string index
Definition: DeMoScan.py:364
LVL1::FEXAlgoSpaceDefs::PS_LOWER_AB
constexpr static int PS_LOWER_AB
Definition: FEXAlgoSpaceDefs.h:173
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
LVL1::gFEXJetAlgo::addRemoteRin
virtual void addRemoteRin(gTowersType &jets, const gTowersPartialSums &partial, int ps_upper, int ps_lower, int ps_shift) const
Definition: gFEXJetAlgo.cxx:1195
LVL1::gTowersJetEngine
std::array< std::array< int, 6 >, 32 > gTowersJetEngine
Definition: gFEXJetAlgo.h:28
defineDB.jets
list jets
Definition: JetTagCalibration/share/defineDB.py:24
LVL1::gTowersType
std::array< std::array< int, 12 >, 32 > gTowersType
Definition: IgFEXFPGA.h:25
LVL1::gFEXJetAlgo::RemotePartialCN
virtual void RemotePartialCN(const gTowersJetEngine &twrs, gTowersPartialSums &rps) const
Definition: gFEXJetAlgo.cxx:712
LVL1::FEXAlgoSpaceDefs::fineFloor
constexpr static int fineFloor
Definition: FEXAlgoSpaceDefs.h:180
AthAlgTool
Definition: AthAlgTool.h:26
LVL1::gFEXJetAlgo::addInternalLin
virtual void addInternalLin(gTowersType &jets, gTowersPartialSums &partial) const
Definition: gFEXJetAlgo.cxx:880
gTower.h
LVL1::gFEXJetAlgo::addRemoteCPin
virtual void addRemoteCPin(gTowersType &jets, const gTowersPartialSums &partial, int ps_upper, int ps_lower, int ps_shift) const
Definition: gFEXJetAlgo.cxx:1277
LVL1::FEXAlgoSpaceDefs::n_partial
constexpr static int n_partial
Definition: FEXAlgoSpaceDefs.h:154
LVL1::gFEXJetAlgo::gBlockMax192
virtual void gBlockMax192(const gTowersJetEngine &gBlkSum, std::array< int, 3 > &gBlockVp, std::array< int, 3 > &gBlockEtap, std::array< int, 3 > &gBlockPhip, int index) const
Definition: gFEXJetAlgo.cxx:1055
Generated on Tue Dec 24 2024 21:10:55 for ATLAS Offline Software by   doxygen 1.8.18