jFEXForwardJetsAlgo.cxx 17.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
/*
  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration  
*/
//***************************************************************************  
//		jFEXForwardJetsAlgo - Algorithm for forward Jets in jFEX
//                              -------------------
//     begin                : 07 06 2021
//     email                : varsiha.sothilingam@cern.ch
//***************************************************************************  
#include <iostream>
#include <vector>
#include <string>
#include <map>
#include "L1CaloFEXSim/jFEXForwardJetsAlgo.h"
#include "L1CaloFEXSim/jTower.h"
#include "L1CaloFEXSim/jTowerContainer.h"
#include "L1CaloFEXSim/jFEXForwardJetsInfo.h"
#include "L1CaloFEXSim/FEXAlgoSpaceDefs.h"
#include "CaloEvent/CaloCellContainer.h"
#include "CaloIdentifier/CaloIdManager.h"
#include "CaloIdentifier/CaloCell_SuperCell_ID.h"
#include "AthenaBaseComps/AthAlgorithm.h"
#include "StoreGate/StoreGateSvc.h"

namespace LVL1{

//Default Constructor
LVL1::jFEXForwardJetsAlgo::jFEXForwardJetsAlgo(const std::string& type, const std::string& name, const IInterface* parent):
   AthAlgTool(type, name, parent)   
  {
  declareInterface<IjFEXForwardJetsAlgo>(this);
  }

/** Destructor */
LVL1::jFEXForwardJetsAlgo::~jFEXForwardJetsAlgo()
{
}
StatusCode LVL1::jFEXForwardJetsAlgo::initialize()
{
   ATH_CHECK(m_jFEXForwardJetsAlgo_jTowerContainerKey.initialize());
   return StatusCode::SUCCESS;

}

//calls container for TT
46
StatusCode LVL1::jFEXForwardJetsAlgo::safetyTest() {
47

48
    SG::ReadHandle<jTowerContainer> jk_jFEXForwardJetsAlgo_jTowerContainer(m_jFEXForwardJetsAlgo_jTowerContainerKey);
49

50
51
    if(! jk_jFEXForwardJetsAlgo_jTowerContainer.isValid()) {
        ATH_MSG_FATAL("Could not retrieve  jk_jFEXForwardJetsAlgo_jTowerContainer " << m_jFEXForwardJetsAlgo_jTowerContainerKey.key());
52

53
54
55
56
        return StatusCode::FAILURE;
    }

    return StatusCode::SUCCESS;
57
58
}

59
60
61
62
void LVL1::jFEXForwardJetsAlgo::setup(int inputTable[FEXAlgoSpaceDefs::jFEX_algoSpace_height][FEXAlgoSpaceDefs::jFEX_wide_algoSpace_width], int jfex, int fpga) {
    std::copy(&inputTable[0][0], &inputTable[0][0] + (FEXAlgoSpaceDefs::jFEX_algoSpace_height*FEXAlgoSpaceDefs::jFEX_wide_algoSpace_width), &m_jFEXalgoTowerID[0][0]);
    m_jfex=jfex;
    m_fpga=fpga;
63
64
65
}

//Gets geometric global centre Phi coord of the TT
66
float LVL1::jFEXForwardJetsAlgo::globalPhi(int nphi, int neta, const jTowerContainer& towers) const {
67
68
69
70
    int TTID = m_jFEXalgoTowerID[nphi][neta];
    if(TTID == 0) {
        return 999;
    }
71

72
    float phi = towers.findTower(TTID)->centrePhi();
73
    return phi;
74
75
76
}

//Gets geometric global centre Eta coord of the TT
77
float LVL1::jFEXForwardJetsAlgo::globalEta(int nphi, int neta, const jTowerContainer& towers) const {
78
79
80
81
    int TTID = m_jFEXalgoTowerID[nphi][neta];
    if(TTID == 0) {
        return 999;
    }
82

83
    float eta = towers.findTower(TTID)->centreEta();
84
    return eta;
85
86
87
}

//Gets Phi of the TT
88
89
unsigned int LVL1::jFEXForwardJetsAlgo::localPhi(int nphi, int neta, const jTowerContainer& towers) const {
    unsigned int phi = towers.findTower(m_jFEXalgoTowerID[nphi][neta])->phi();
90
    return phi;
91
}
92

93
//Gets Eta of the TT
94
unsigned int LVL1::jFEXForwardJetsAlgo::localEta(int nphi, int neta, const jTowerContainer& towers) const {
95

96
    unsigned int eta  = towers.findTower(m_jFEXalgoTowerID[nphi][neta])->eta();
97
    return eta;
98
99
}

100
//Return ET of TT. Should be FCAL 0 + 1 + 2 //maybe check this
101
unsigned int LVL1::jFEXForwardJetsAlgo::getTTowerET(int nphi, int neta ) {
102
103
104
    unsigned int TTID = m_jFEXalgoTowerID[nphi][neta];
    if(TTID == 0) {
        return 0;
105
106
    }

107
108
109
    if(m_map_Etvalues.find(TTID) != m_map_Etvalues.end()) {
        return m_map_Etvalues[TTID][0];
    }
110

111
    //we shouldn't arrive here
112
    return 0;
113
114
}

115
std::unordered_map<int, jFEXForwardJetsInfo> LVL1::jFEXForwardJetsAlgo::FcalJetsTowerIDLists() {
116
    
117
    std::unordered_map<int, jFEXForwardJetsInfo> FCALJetTowerIDLists;
118

119
120
121
122
    std::vector<int> lower_centre_neta;
    std::vector<int> upper_centre_neta;
    m_lowerEM_eta = 0;
    m_upperEM_eta = 0;
123

124
125
126
127
128
129
130
131
132
133
134
135
    //STEP 1: check if we are in module 0 or 5 and assign corrrect eta FEXAlgoSpace parameters
    if(m_jfex == 0) {
        //Module 0
        lower_centre_neta.assign({FEXAlgoSpaceDefs::jFEX_algoSpace_C_EMB_start_eta, FEXAlgoSpaceDefs::jFEX_algoSpace_C_EMIE_start_eta, FEXAlgoSpaceDefs::jFEX_algoSpace_C_FCAL_start_eta});
        upper_centre_neta.assign({FEXAlgoSpaceDefs::jFEX_algoSpace_C_EMB_end_eta, FEXAlgoSpaceDefs::jFEX_algoSpace_C_EMIE_end_eta,FEXAlgoSpaceDefs::jFEX_algoSpace_C_FCAL_end_eta });
        m_lowerEM_eta = FEXAlgoSpaceDefs::jFEX_algoSpace_C_lowerEM_eta;
        m_upperEM_eta = FEXAlgoSpaceDefs::jFEX_algoSpace_C_upperEM_eta;
    }
    else {
        //Module 5
        lower_centre_neta.assign({FEXAlgoSpaceDefs::jFEX_algoSpace_A_EMB_eta, FEXAlgoSpaceDefs::jFEX_algoSpace_A_EMIE_eta, FEXAlgoSpaceDefs::jFEX_algoSpace_A_FCAL_start_eta});
        upper_centre_neta.assign({FEXAlgoSpaceDefs::jFEX_algoSpace_A_EMIE_eta, FEXAlgoSpaceDefs::jFEX_algoSpace_A_FCAL_start_eta, FEXAlgoSpaceDefs::jFEX_algoSpace_A_FCAL_end_eta});
136

137
138
        m_lowerEM_eta = FEXAlgoSpaceDefs::jFEX_algoSpace_A_lowerEM_eta;
        m_upperEM_eta = FEXAlgoSpaceDefs::jFEX_algoSpace_A_upperEM_eta;
139

140
    }
141

142
143
144
145
146
    //STEP 2: define phi FEXAlgoSpace parameters
    std::vector<int> lower_centre_nphi{FEXAlgoSpaceDefs::jFEX_algoSpace_EMB_start_phi, FEXAlgoSpaceDefs::jFEX_algoSpace_EMIE_start_phi,  FEXAlgoSpaceDefs::jFEX_algoSpace_FCAL_start_phi};
    std::vector<int> upper_centre_nphi{FEXAlgoSpaceDefs::jFEX_algoSpace_EMB_end_phi, FEXAlgoSpaceDefs::jFEX_algoSpace_EMIE_end_phi,  FEXAlgoSpaceDefs::jFEX_algoSpace_FCAL_end_phi};

    //STEP 3: loop over different EM/FCAL0 eta phi core fpga regions. These are potential central trigger towers for jets
147
    SG::ReadHandle<jTowerContainer> towers(m_jFEXForwardJetsAlgo_jTowerContainerKey/*,ctx*/);
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
    for(int myCounter = 0; myCounter<3; myCounter++) {
        for(int centre_nphi = lower_centre_nphi[myCounter]; centre_nphi < upper_centre_nphi[myCounter]; centre_nphi++) {
            for(int centre_neta = lower_centre_neta[myCounter]; centre_neta < upper_centre_neta[myCounter]; centre_neta++) {

                //STEP 4: define TTID which will be the key for class in map
                int myTTIDKey = m_jFEXalgoTowerID[centre_nphi][centre_neta];

                //STEP 5: ignore when tower ID is zero. Should not happend though
                if(myTTIDKey == 0) {
                    continue;
                }

                //STEP 6: define class
                jFEXForwardJetsInfo TriggerTowerInformation;

                //STEP 7: fill local eta phi coords
                TriggerTowerInformation.setCentreLocalTTPhi(centre_nphi);
                TriggerTowerInformation.setCentreLocalTTEta(centre_neta);

167
168
                float centreTT_phi = globalPhi(centre_nphi, centre_neta, *towers);
                float centreTT_eta = globalEta(centre_nphi, centre_neta, *towers);
169
170
171
172
173
174
175
176
177
178
179
180
181
182
                TriggerTowerInformation.setCentreTTPhi(centreTT_phi);
                TriggerTowerInformation.setCentreTTEta(centreTT_eta);
                TriggerTowerInformation.includeTTinSeed(m_jFEXalgoTowerID[centre_nphi][centre_neta]);
                TriggerTowerInformation.addToSeedET(getTTowerET(centreTT_phi, centreTT_eta));
                TriggerTowerInformation.includeTTinSearchWindow(m_jFEXalgoTowerID[centre_nphi][centre_neta]);
                TriggerTowerInformation.addToSearchWindowET(getTTowerET(centreTT_phi, centreTT_eta));
                
                //STEP 8: loop over all FCAL0 TTIDs and fill TT IDs for seed and calculate seed energy
                for(int nphi = 0; nphi < FEXAlgoSpaceDefs::jFEX_algoSpace_height; nphi++) {
                    for(int neta = m_lowerEM_eta; neta < m_upperEM_eta; neta++) {
                        
                        if(m_jFEXalgoTowerID[centre_nphi][centre_neta]==m_jFEXalgoTowerID[nphi][neta]  || m_jFEXalgoTowerID[nphi][neta]==0 ) {
                            continue;
                        }
183
                        float TT_phi = globalPhi(nphi, neta, *towers);
184
185
186
187
188
189
190
191
192
193
194
195
                        if(m_fpga==0 || m_fpga==3) { //This corrects the overlap of FPGA 0 with FPGA 3 and viceversa
                            if(m_fpga==0) {
                                if(TT_phi>M_PI){
                                   TT_phi = TT_phi-6.4; 
                                }
                            }
                            else {
                                if(TT_phi<M_PI){
                                   TT_phi = TT_phi+6.4; 
                                }
                            }
                        }
196
                        float TT_eta = globalEta(nphi, neta, *towers);
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
                        int DeltaR = std::round(std::sqrt(std::pow((centreTT_eta - TT_eta),2) + std::pow((centreTT_phi - TT_phi),2))*100); // cast float to int to avoid misbehaviours
                        if(DeltaR < 20 ) {
                            //STEP 9.0: fill TTID in seed
                            TriggerTowerInformation.includeTTinSeed(m_jFEXalgoTowerID[nphi][neta]);
                            //STEP 10.0: add ET value to seed
                            TriggerTowerInformation.addToSeedET(getTTowerET(nphi, neta));
                        }
                        else if(DeltaR < 40 ) {
                            TriggerTowerInformation.addToFirstEnergyRingET(getTTowerET(nphi, neta));
                            if(m_storeEnergyRingTTIDs) {
                                TriggerTowerInformation.includeTTIDinFirstER(m_jFEXalgoTowerID[nphi][neta]);
                            }
                        }
                        else if(DeltaR < 80){
                            TriggerTowerInformation.addToSecondEnergyRingET(getTTowerET(nphi, neta));
                            if(m_storeEnergyRingTTIDs) {
                                TriggerTowerInformation.includeTTIDinSecondER(m_jFEXalgoTowerID[nphi][neta]);
                            }    
                        }
                        
                        
                        if(DeltaR < 30 ) {
                            //STEP 9.1: fill TTID in search window
                            TriggerTowerInformation.includeTTinSearchWindow(m_jFEXalgoTowerID[nphi][neta]);
                            //STEP 10.1: add ET value to seed
                            TriggerTowerInformation.addToSearchWindowET(getTTowerET(nphi, neta));
                        }
                    }
                }
                //STEP 11: Store filled class into map, with central TT ID as Key
                FCALJetTowerIDLists[myTTIDKey] = TriggerTowerInformation;
            }//end of centre_neta loop
        }//end of centre_nphi loop
    }//end of myCounter loop
231
    
232
    return FCALJetTowerIDLists;
233
234
235
236
}



237
std::unordered_map<int, jFEXForwardJetsInfo> LVL1::jFEXForwardJetsAlgo::isSeedLocalMaxima() {
238
    //std::vector<int> localMaximas;
239
240
    std::unordered_map<int, jFEXForwardJetsInfo> localMaximaCandidates = FcalJetsTowerIDLists();
    std::unordered_map<int, jFEXForwardJetsInfo> localMaximaList ;
241
    SG::ReadHandle<jTowerContainer> my_jTowerContainer(m_jFEXForwardJetsAlgo_jTowerContainerKey/*,ctx*/);
242

243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
    size_t isLocalMaxima = 0;
    
    for (std::pair<int, jFEXForwardJetsInfo> element : localMaximaCandidates){
        
        int myTTKey = element.first;
        jFEXForwardJetsInfo myFCALJetInfoClass = element.second;
        
        //Local maxima check takes place here
        isLocalMaxima = 0;
        float centre_phi = myFCALJetInfoClass.getCentreTTPhi();
        float centre_eta = myFCALJetInfoClass.getCentreTTEta();
        int centre_energy = myFCALJetInfoClass.getSeedET();

        const std::vector<int> TTinSW = myFCALJetInfoClass.getTTinSearchWindow();
        for (const int iTTinSW : TTinSW) {
            if(iTTinSW == myTTKey) continue;
            const LVL1::jTower* seed_tower = my_jTowerContainer->findTower(iTTinSW);
            float seed_phi = seed_tower->centrePhi();
            float seed_eta = seed_tower->centreEta();

            int delta_phi = std::round((seed_phi - centre_phi)*100);
            int delta_eta = std::round((seed_eta - centre_eta)*100);
            
            int seed_energy = 0;
            
            //get the seed Et for iTTinSW, and check if it is a Local Maxima 
            if(localMaximaCandidates.count(iTTinSW)==1 ){
              seed_energy = localMaximaCandidates[iTTinSW].getSeedET();
            }
            else{
                for(int nphi = 0; nphi < FEXAlgoSpaceDefs::jFEX_algoSpace_height; nphi++) {
                    for(int neta = m_lowerEM_eta; neta < m_upperEM_eta; neta++) {
                        
                        if(m_jFEXalgoTowerID[nphi][neta]==0 ) {
                            continue;
                        }
279
                        float TT_phi = globalPhi(nphi, neta, *my_jTowerContainer);
280
281
282
283
284
285
286
287
288
289
290
291
                        if(m_fpga==0 || m_fpga==3) { //This corrects the overlap of FPGA 0 with FPGA 3 and viceversa
                            if(m_fpga==0) {
                                if(TT_phi>M_PI){
                                   TT_phi = TT_phi-6.4; 
                                }
                            }
                            else {
                                if(TT_phi<M_PI){
                                   TT_phi = TT_phi+6.4; 
                                }
                            }
                        }
292
                        float TT_eta = globalEta(nphi, neta, *my_jTowerContainer);
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
                        int DeltaR = std::round(std::sqrt(std::pow((seed_eta - TT_eta),2) + std::pow((seed_phi - TT_phi),2))*100); // cast float to int to avoid misbehaviours
                        if(DeltaR < 20 ) {
                            seed_energy+=getTTowerET(nphi, neta);
                        }
                    }
                }
            }
            
            if( (delta_eta + delta_phi) < 0 || ((delta_eta + delta_phi) == 0 && delta_eta < 0) ) {
                if(centre_energy >= seed_energy) {
                    isLocalMaxima++;
                }
            }
            else {   
                if(centre_energy > seed_energy) {
                    isLocalMaxima++;
                }
310
311
            }
        }
312
313
314
315
        //if it is a local maxima, here we save the TT ID
        if((isLocalMaxima == TTinSW.size()-1) && (isLocalMaxima !=0)) {
            localMaximaList[myTTKey] = myFCALJetInfoClass;
        }
316
    }
317
    return localMaximaList;
318
319
}

320

321
std::unordered_map<int, jFEXForwardJetsInfo> LVL1::jFEXForwardJetsAlgo::calculateJetETs() {
322

323
324
325
326
327
328
329
330
331
    // setting the lower/upper eta range for the FCAL 2 and 3 since they are not added in the seed information yet 
    int lowerFCAL_eta = FEXAlgoSpaceDefs::jFEX_algoSpace_C_lowerFCAL_eta;
    int upperFCAL_eta = FEXAlgoSpaceDefs::jFEX_algoSpace_C_upperFCAL_eta;
    if(m_jfex == 5) {
        //Module 5
        lowerFCAL_eta = FEXAlgoSpaceDefs::jFEX_algoSpace_A_lowerFCAL_eta;
        upperFCAL_eta = FEXAlgoSpaceDefs::jFEX_algoSpace_A_upperFCAL_eta;
    } 
    // Adding the FCAL 2 and 3 TT in the seed, 1st and 2nd energy rings
332
    SG::ReadHandle<jTowerContainer> towers(m_jFEXForwardJetsAlgo_jTowerContainerKey/*,ctx*/);
333
    std::unordered_map<int, jFEXForwardJetsInfo> localMaximas = isSeedLocalMaxima();
334
335
336
337
338
339
340
341
    for(std::pair<int, jFEXForwardJetsInfo> element : localMaximas) {
        jFEXForwardJetsInfo myFCALJetInfoClass = element.second;

        float centreTT_phi = myFCALJetInfoClass.getCentreTTPhi();
        float centreTT_eta = myFCALJetInfoClass.getCentreTTEta();
        for(int nphi = 0; nphi < 8; nphi++) {
            for(int neta = lowerFCAL_eta; neta < upperFCAL_eta; neta++) {
                
342
343
                float TT_phi = globalPhi(nphi, neta, *towers );
                float TT_eta = globalEta(nphi, neta, *towers);
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374

                if(m_fpga==0 || m_fpga==3) { //This corrects the overlap of FPGA 0 with FPGA 3 and viceversa
                    if(m_fpga==0) {
                        if(TT_phi>M_PI) {
                            TT_phi = TT_phi-6.4;
                        }
                    }
                    else {
                        if(TT_phi<M_PI) {
                            TT_phi = TT_phi+6.4;
                        }
                    }
                }
                int DeltaR = std::round(std::sqrt(std::pow((centreTT_eta - TT_eta),2) + std::pow((centreTT_phi - TT_phi),2))*100); // cast float to int to avoid misbehaviours
                if(DeltaR < 20 ) {
                    myFCALJetInfoClass.addToSeedET(getTTowerET(nphi, neta));
                    myFCALJetInfoClass.includeTTinSeed(m_jFEXalgoTowerID[nphi][neta]);
                }
                else if(DeltaR < 40 ) {
                    myFCALJetInfoClass.addToFirstEnergyRingET(getTTowerET(nphi, neta));
                    if(m_storeEnergyRingTTIDs) {
                        myFCALJetInfoClass.includeTTIDinFirstER(m_jFEXalgoTowerID[nphi][neta]);
                    }
                }
                else if(DeltaR < 80) {
                    myFCALJetInfoClass.addToSecondEnergyRingET(getTTowerET(nphi, neta));
                    if(m_storeEnergyRingTTIDs) {
                        myFCALJetInfoClass.includeTTIDinSecondER(m_jFEXalgoTowerID[nphi][neta]);
                    }
                }
            }
375
        }
376
        
377
    }
378
379
    //Returning the complite information 
    return localMaximas;
380
381
}

382

383
void LVL1::jFEXForwardJetsAlgo::setFPGAEnergy(std::unordered_map<int,std::vector<int> > et_map){
384
385
386
    m_map_Etvalues=et_map;
}

387
388
}// end of namespace LVL1