ATLAS Offline Software
Loading...
Searching...
No Matches
PtEndcapLUT.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "PtEndcapLUT.h"
6#include "CLHEP/Vector/TwoVector.h"
7#include <fstream>
8#include <sstream>
9
10#include "AthenaBaseComps/AthMsgStreamMacros.h"
11
12// --------------------------------------------------------------------------------
13// --------------------------------------------------------------------------------
14
15TrigL2MuonSA::PtEndcapLUT::PtEndcapLUT(const std::string& type,
16 const std::string& name,
17 const IInterface* parent):
18 AthAlgTool(type, name, parent)
19{
20}
21
22// --------------------------------------------------------------------------------
23// --------------------------------------------------------------------------------
24
25TrigL2MuonSA::PtEndcapLUT::~PtEndcapLUT()
26{
27 for (TableMap::iterator it = m_tables.begin(); it != m_tables.end(); ++it) {
28 delete it->second;
29 }
30}
31
32// --------------------------------------------------------------------------------
33// --------------------------------------------------------------------------------
34
35bool TrigL2MuonSA::PtEndcapLUT::KeyType::operator<(const KeyType& other) const
36{
37 if (m_side < other.m_side) return (true);
38 if (m_side > other.m_side) return (false);
39 if (m_charge < other.m_charge) return (true);
40 if (m_charge > other.m_charge) return (false);
41 if (m_type < other.m_type) return (true);
42 if (m_type > other.m_type) return (false);
43 return (false);
44}
45
46// --------------------------------------------------------------------------------
47// --------------------------------------------------------------------------------
48
49std::string TrigL2MuonSA::PtEndcapLUT::KeyType::toString() const
50{
51 std::ostringstream ss;
52 ss << "side=" << m_side << " charge=" << m_charge << " type=" << dt2s(m_type);
53 return ss.str();
54}
55
56// --------------------------------------------------------------------------------
57// --------------------------------------------------------------------------------
58
59StatusCode TrigL2MuonSA::PtEndcapLUT::readLUT(const std::string& lut_fileName)
60{
61 std::ifstream ifs(lut_fileName.c_str());
62 if (!ifs.is_open()) {
63 ATH_MSG_ERROR("Cannot open EndcapLUT file " << lut_fileName);
64 return StatusCode::FAILURE;
65 }
66
67 m_tables.clear();
68 std::string line;
69 int line_no = 0;
70 TableType* table = NULL;
71
72 while (!ifs.eof()) {
73
74 getline(ifs, line);
75 line_no++;
76
77 if (line.empty()) continue;
78
79 if (line.substr(0, 5) == "side=") {
80 char side, charge, ctype[15];
81
82 if (sscanf(line.c_str(), "side=%c charge=%c %14s", &side, &charge, ctype) != 3) {
83 ATH_MSG_ERROR("Invalid header line " << line_no << " in EndcapLUT file " << lut_fileName);
84 return StatusCode::FAILURE;
85 }
86
87 DataType type = s2dt(ctype);
88 if ((side != '-' && side != '+') || (charge != '-' && charge != '+') || type == INVALID) {
89 ATH_MSG_ERROR("Invalid header line " << line_no << " in EndcapLUT file " << lut_fileName);
90 return StatusCode::FAILURE;
91 }
92
93 KeyType key(side == '+' ? 1 : 0, charge == '+' ? 1 : 0, type);
94 table = new TableType();
95 m_tables.insert(TableMap::value_type(key, table));
96 ATH_MSG_DEBUG("Created table for " << key.toString());
97
98 } else {
99
100 if (table == NULL) {
101 ATH_MSG_ERROR("Missing header line " << line_no << " in EndcapLUT file " << lut_fileName);
102 return StatusCode::FAILURE;
103 }
104
105 int iEta, iPhi;
106 double xcept, slope;
107 if (sscanf(line.c_str(), "%3d %3d %15lf %15lf", &iEta, &iPhi, &xcept, &slope) != 4) {
108 ATH_MSG_ERROR("Invalid data line " << line_no << " in EndcapLUT file " << lut_fileName);
109 return StatusCode::FAILURE;
110 }
111 table->m_xcepts[iEta-1][iPhi] = xcept;
112 table->m_slopes[iEta-1][iPhi] = slope;
113 }
114 }
115 ifs.close();
116
117 return StatusCode::SUCCESS;
118}
119
120// --------------------------------------------------------------------------------
121// --------------------------------------------------------------------------------
122
123double TrigL2MuonSA::PtEndcapLUT::alpha(double z1, double r1, double z2, double r2) const
124{
125 CLHEP::Hep2Vector mid1(z1, r1);
126 CLHEP::Hep2Vector mid2(z2, r2);
127 CLHEP::Hep2Vector dir = mid2 - mid1;
128
129 if(std::abs(mid1.unit() * dir.unit()) > 1.) return 0;
130
131 double a = std::acos(mid1.unit() * dir.unit());
132
133 return a;
134}
135
136// --------------------------------------------------------------------------------
137// --------------------------------------------------------------------------------
138
139double TrigL2MuonSA::PtEndcapLUT::radius(double z1, double r1, double s1, double z2, double r2, double s2, double deltar) const
140{
141 float cr1 = 0.080/400;
142 float cr2 = cr1;
143 double x1 = z1;
144
145 if (std::abs(x1)>=0.1) {
146 double x2 = z2;
147 double y1 = r1;
148 double y2 = r2;
149 double A1 = s1;
150 double A2 = s2;
151 if(!(std::abs(s2+999)<0.1)) {
152 A2 = s2;
153 cr2 = cr1/10;
154 }
155
156 // find centre of circonference
157 double xm = (x1+x2)/2.;
158 double ym = (y1+y2)/2.;
159 double c1 = (x2-x1)*xm+(y2-y1)*ym;
160 double c2_1 = -x1-A1*y1;
161 double c2_2 = -x2-A2*y2;
162 double yR1 = (-c1-c2_1*(x2-x1))/(A1*(x2-x1)-(y2-y1));
163 double yR2 = (-c1-c2_2*(x2-x1))/(A2*(x2-x1)-(y2-y1));
164 double xR1 = -A1*yR1-c2_1;
165 double xR2 = -A2*yR2-c2_2;
166
167 double xR = ((1./cr1)*xR1+(1./cr2)*xR2)/((1./cr1)+(1./cr2));
168 double yR = ((1./cr1)*yR1+(1./cr2)*yR2)/((1./cr1)+(1./cr2));
169 double sign = deltar / std::abs(deltar);
170 double radius = 0.5*(std::sqrt((xR-x1)*(xR-x1)+(yR-y1)*(yR-y1))+std::sqrt((xR-x2)*(xR-x2)+(yR-y2)*(yR-y2)));
171 return(sign * radius);
172 } else {
173 return 0.0;
174 }
175}
176
177// --------------------------------------------------------------------------------
178// --------------------------------------------------------------------------------
179
180double TrigL2MuonSA::PtEndcapLUT::lookup(int side, int charge, DataType type, int iEta, int iPhi, double value) const
181{
182 ATH_MSG_DEBUG("lookup(side=" << side << ",charge=" << charge
183 << ",type=" << dt2s(type) << ",iEta=" << iEta
184 << ",iPhi=" << iPhi << ",value=" << value << ")");
185
186 if (iEta == -1) iEta = 0;
187 if (iEta == 30) iEta = 29;
188
189 if (iEta < 0 || iEta >= ETAS || iPhi < 0 || iPhi >= PHISEE) {
190 ATH_MSG_WARNING("lookup(" << side << ", " << charge << ", " << dt2s(type)
191 << ", " << iEta << ", " << iPhi << ") Invalid indices");
192 return 0.0;
193 }
194
195 TableMap::const_iterator it = m_tables.find(KeyType(side, charge, type));
196
197 if (it == m_tables.end()) {
198 ATH_MSG_ERROR("lookup(" << side << ", " << charge << ", " << dt2s(type)
199 << ", " << iEta << ", " << iPhi << ") Invalid key");
200 return 0.0;
201 }
202
203 double pt = 0;
204 TableType* table = it->second;
205
206 double b = table->m_xcepts[iEta][iPhi];
207 double c = table->m_slopes[iEta][iPhi];
208 const double PT_MAX = 500; // 500 GeV upper limit
209 const double ZERO_LIMIT = 1e-5;
210
211 if( std::abs(c) < ZERO_LIMIT ) {
212
213 if( std::abs(b) < ZERO_LIMIT ) {
214
215 pt = PT_MAX;
216
217 } else {
218
219 double ptinv = value / b;
220
221 if( ptinv < 0 ) {
222
223 ATH_MSG_WARNING("Pol2: ptinv < 0, ptinv=" << ptinv);
224
225 } else {
226
227 if( ptinv < (1/PT_MAX) ) pt = PT_MAX;
228 else { pt = 1/ptinv; }
229
230 }
231 }
232
233 } else {
234
235 double judge = b*b + 4*c*value;
236
237 if( judge > 0 ) {
238
239 double ptinv = (-b + pow(judge,0.5)) / c / 2;
240
241 if( ptinv < 0 ) {
242
243 ATH_MSG_WARNING("Pol2: ptinv < 0, ptinv=" << ptinv);
244
245 } else {
246
247 if( ptinv < (1/PT_MAX) ) {
248 pt = PT_MAX;
249 } else {
250 pt = 1/ptinv;
251 }
252 }
253 }
254 }
255
256 ATH_MSG_DEBUG("Pol2: value=" << value);
257 ATH_MSG_DEBUG("Pol2: b=" << b << " c=" << c << " pt=" << pt);
258 pt *= 1000; // convert to MeV
259
260 return pt;
261}
262
263// --------------------------------------------------------------------------------
264// --------------------------------------------------------------------------------
265
266TrigL2MuonSA::PtEndcapLUT::DataType TrigL2MuonSA::PtEndcapLUT::s2dt(const char* type)
267{
268 std::string stype(type);
269 if (stype == "alphapol2") return (ALPHAPOL2);
270 if (stype == "betapol2") return (BETAPOL2);
271 if (stype == "tgcalphapol2") return (TGCALPHAPOL2);
272 if (stype == "invradiuspol2") return (INVRADIUSPOL2);
273 if (stype == "cscpol2") return (CSCPOL2);
274 return (INVALID);
275}
276
277// --------------------------------------------------------------------------------
278// --------------------------------------------------------------------------------
279
280const char* TrigL2MuonSA::PtEndcapLUT::dt2s(TrigL2MuonSA::PtEndcapLUT::DataType type)
281{
282 switch (type) {
283
284 case ALPHAPOL2:
285 return ("alphapol2"); break;
286
287 case BETAPOL2:
288 return ("betapol2"); break;
289
290 case TGCALPHAPOL2:
291 return ("tgcalphapol2"); break;
292
293 case INVRADIUSPOL2:
294 return ("invradiuspol2"); break;
295
296 case CSCPOL2:
297 return ("cscpol2"); break;
298
299 case INVALID:
300 return ("invalid"); break;
301
302 default:
303 return ("invalid"); break;
304 }
305}
306
307// --------------------------------------------------------------------------------
308// --------------------------------------------------------------------------------
309
310double TrigL2MuonSA::PtEndcapLUT::ptcombined(int iEta, int iPhi, double ApT, double BpT, double &CApT, \
311 double &CBpT) const
312{
313 ATH_MSG_DEBUG("pTcombined("
314 << "iEta=" << iEta
315 << "iPhi=" << iPhi
316 << "Alpha pT=" << ApT
317 << "Beta pT=" << BpT
318 << ")" );
319
320 if (iEta == -1) iEta = 0;
321 if (iEta == 30) iEta = 29;
322
323 if (iEta < 0 || iEta >= ETAS || iPhi < 0 || iPhi >= PHIS) {
324 ATH_MSG_WARNING("pTcombined("<< iEta << ", " << iPhi << ") Invalid indices");
325 return 0.0;
326 }
327
328 const float ZERO_LIMIT = 1e-5;
329 int iphibin=iPhi;
330 int ietabin=iEta/6;
331 if(iPhi==5||iPhi==6)iphibin=5;
332 if(iPhi==4||iPhi==7)iphibin=4;
333 if(iPhi==3||iPhi==8)iphibin=3;
334 if(iPhi==2||iPhi==9)iphibin=2;
335 if(iPhi==1||iPhi==10)iphibin=1;
336 if(iPhi==0||iPhi==11)iphibin=0;
337 CApT=0.;
338 CBpT=0.;
339 double Ameana = m_meana[ietabin][iphibin][0];
340 double Ameanb = m_meanb[ietabin][iphibin][0];
341 double Ameanc = m_meanc[ietabin][iphibin][0];
342 double Bmeana = m_meana[ietabin][iphibin][1];
343 double Bmeanb = m_meanb[ietabin][iphibin][1];
344 double Bmeanc = m_meanc[ietabin][iphibin][1];
345
346 double Asigmaa = m_sigmaa[ietabin][iphibin][0];
347 double Asigmab = m_sigmab[ietabin][iphibin][0];
348 double Asigmac = m_sigmac[ietabin][iphibin][0];
349 double Bsigmaa = m_sigmaa[ietabin][iphibin][1];
350 double Bsigmab = m_sigmab[ietabin][iphibin][1];
351 double Bsigmac = m_sigmac[ietabin][iphibin][1];
352
353
354 double MeanAP = ( std::abs(ApT) > ZERO_LIMIT)? (Ameana + Ameanb * std::exp( Ameanc / ApT)) : 1.0;
355 double MeanBP = ( std::abs(BpT) > ZERO_LIMIT)? (Bmeana + Bmeanb * std::exp( Bmeanc / BpT)) : 1.0;
356 double ApT_tmp = ( std::abs(1-MeanAP) > ZERO_LIMIT)? (std::abs(ApT) / (1-MeanAP)) : 1.0;
357 ApT_tmp = std::abs(ApT_tmp);
358 if(ApT_tmp >= 500) ApT_tmp = 500;
359 double BpT_tmp = ( std::abs(1-MeanBP) > ZERO_LIMIT)? (std::abs(BpT) / (1-MeanBP)) : 1.0;
360 BpT_tmp = std::abs(BpT_tmp);
361 if(BpT_tmp >= 500) BpT_tmp = 500;
362
363 CApT = ApT_tmp;
364 CBpT = BpT_tmp;
365
366 if(ApT == 0. ) CApT = 0.;
367 if(BpT == 0. ) CBpT = 0.;
368 double NSigmaA= Asigmaa * std::abs(ApT_tmp) * std::abs(ApT_tmp) + Asigmab * std::abs(ApT_tmp) + Asigmac;
369 double NSigmaB= Bsigmaa * std::abs(BpT_tmp) * std::abs(BpT_tmp) + Bsigmab * std::abs(BpT_tmp) + Bsigmac;
370
371 double NVsigpTA =(std::abs(ApT_tmp) > ZERO_LIMIT&& std::abs(NSigmaA) > ZERO_LIMIT)? (1/(NSigmaA * ApT_tmp) ): 1.0;
372 double NVsigpTB =(std::abs(BpT_tmp) > ZERO_LIMIT&& std::abs(NSigmaB) > ZERO_LIMIT)? (1/(NSigmaB * BpT_tmp) ): 1.0;
373 double NVsigAsq =(std::abs(NSigmaA) > ZERO_LIMIT)? (1/(NSigmaA * NSigmaA)) : 1.0;
374 double NVsigBsq =(std::abs(NSigmaB) > ZERO_LIMIT)? (1/(NSigmaB * NSigmaB)) : 1.0;
375
376 double NVsigpTAsq = NVsigpTA * NVsigpTA;
377 double NVsigpTBsq = NVsigpTB * NVsigpTB;
378 double pt = (std::abs(NVsigAsq + NVsigBsq) > ZERO_LIMIT)? (1/std::sqrt((NVsigpTAsq + NVsigpTBsq)/(NVsigAsq + NVsigBsq))) : 0.;
379 if(pt>500) pt = 500.;
380 return pt;
381}
382// --------------------------------------------------------------------------------
383// --------------------------------------------------------------------------------
384
385StatusCode TrigL2MuonSA::PtEndcapLUT::readLUTSigmaMean(const std::string& lut_mean, const std::string& lut_sigma)
386{
387 std::ifstream ifsmean(lut_mean.c_str());
388 std::ifstream ifssigma(lut_sigma.c_str());
389 if (!ifsmean.is_open()) {
390 ATH_MSG_ERROR("Cannot open EndcapLUT Mean file " << lut_mean);
391 return StatusCode::FAILURE;
392 }
393 if (!ifssigma.is_open()) {
394 ATH_MSG_ERROR("Cannot open EndcapLUT Sigma file " << lut_sigma);
395 return StatusCode::FAILURE;
396 }
397
398 std::string line;
399
400 for(int ei=0; ei < ETAS1; ei++){
401 for(int pi=0; pi < PHIS1; pi++){
402 for(int pti=0; pti < PTS1; pti++){
403 m_meana[ei][pi][pti] =0.;
404 m_meanb[ei][pi][pti] =0.;
405 m_meanc[ei][pi][pti] =0.;
406 m_sigmaa[ei][pi][pti]=0.;
407 m_sigmab[ei][pi][pti]=0.;
408 m_sigmac[ei][pi][pti]=0.;
409 }
410 }
411 }
412
413 while (!ifsmean.eof()) {
414 getline(ifsmean, line);
415 if (line.empty()) continue;
416
417 int iEta, iPhi, iNP;
418 double tmp_par1, tmp_par2, tmp_par3;
419 if (sscanf(line.c_str(), "%d %d %d %lf %lf %lf", &iEta, &iPhi, &iNP, &tmp_par1, &tmp_par2, &tmp_par3) != 6) {
420 ATH_MSG_ERROR(" Invalid data in mean EndcapLUT file " << lut_mean);
421 return StatusCode::FAILURE;
422 }
423
424 m_meana[iEta][iPhi][iNP] = tmp_par1;
425 m_meanb[iEta][iPhi][iNP] = tmp_par2;
426 m_meanc[iEta][iPhi][iNP] = tmp_par3;
427 }
428 ifsmean.close();
429 std::string line2;
430 while (!ifssigma.eof()) {
431 getline(ifssigma, line2);
432 if (line2.empty()) continue;
433
434 int iEta, iPhi, iNP;
435 double tmp_par1, tmp_par2, tmp_par3;
436 if (sscanf(line2.c_str(), "%d %d %d %lf %lf %lf", &iEta, &iPhi, &iNP, &tmp_par1, &tmp_par2, &tmp_par3) != 6) {
437 ATH_MSG_ERROR(" Invalid data in mean EndcapLUT file " << lut_mean);
438 return StatusCode::FAILURE;
439 }
440
441 m_sigmaa[iEta][iPhi][iNP] = tmp_par1;
442 m_sigmab[iEta][iPhi][iNP] = tmp_par2;
443 m_sigmac[iEta][iPhi][iNP] = tmp_par3;
444 }
445 ifssigma.close();
446
447 return StatusCode::SUCCESS;
448}
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
charge
double charge(const T &p)
Definition AtlasPID.h:997
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
ss
static Double_t ss
Definition LArPhysWaveHECTool.cxx:37
PTS1
#define PTS1
Definition PtEndcapLUT.h:17
PHIS1
#define PHIS1
Definition PtEndcapLUT.h:16
ETAS1
#define ETAS1
Definition PtEndcapLUT.h:15
sign
int sign(int a)
Definition TRT_StrawNeighbourSvc.h:108
pi
#define pi
Definition TileMuonFitter.cxx:65
ZERO_LIMIT
const float ZERO_LIMIT
Definition VP1TriggerHandleL2.cxx:37
pow
constexpr int pow(int base, int exp) noexcept
Definition ap_fixedTest.cxx:15
AthAlgTool::AthAlgTool
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
Definition AthAlgTool.cxx:16
TrigL2MuonSA::PtEndcapLUT::radius
double radius(double z1, double r1, double s1, double z2, double r2, double s2, double deltar) const
Definition PtEndcapLUT.cxx:139
TrigL2MuonSA::PtEndcapLUT::m_sigmaa
double m_sigmaa[ETAS1][PHIS1][PTS1]
Definition PtEndcapLUT.h:44
TrigL2MuonSA::PtEndcapLUT::ptcombined
double ptcombined(int iEta, int iPhi, double ApT, double BpT, double &CApT, double &CBpT) const
Definition PtEndcapLUT.cxx:310
TrigL2MuonSA::PtEndcapLUT::dt2s
static const char * dt2s(DataType type)
Definition PtEndcapLUT.cxx:280
TrigL2MuonSA::PtEndcapLUT::PtEndcapLUT
PtEndcapLUT(const std::string &type, const std::string &name, const IInterface *parent)
Definition PtEndcapLUT.cxx:15
TrigL2MuonSA::PtEndcapLUT::m_sigmab
double m_sigmab[ETAS1][PHIS1][PTS1]
Definition PtEndcapLUT.h:45
TrigL2MuonSA::PtEndcapLUT::lookup
double lookup(int side, int charge, DataType type, int iEta, int iPhi, double value) const
Definition PtEndcapLUT.cxx:180
TrigL2MuonSA::PtEndcapLUT::m_sigmac
double m_sigmac[ETAS1][PHIS1][PTS1]
Definition PtEndcapLUT.h:46
TrigL2MuonSA::PtEndcapLUT::alpha
double alpha(double z1, double r1, double z2, double r2) const
Definition PtEndcapLUT.cxx:123
TrigL2MuonSA::PtEndcapLUT::readLUT
StatusCode readLUT(const std::string &lut_fileName)
Definition PtEndcapLUT.cxx:59
TrigL2MuonSA::PtEndcapLUT::readLUTSigmaMean
StatusCode readLUTSigmaMean(const std::string &lut_mean, const std::string &lut_sigma)
Definition PtEndcapLUT.cxx:385
TrigL2MuonSA::PtEndcapLUT::s2dt
static DataType s2dt(const char *type)
Definition PtEndcapLUT.cxx:266
TrigL2MuonSA::PtEndcapLUT::m_meanc
double m_meanc[ETAS1][PHIS1][PTS1]
Definition PtEndcapLUT.h:43
TrigL2MuonSA::PtEndcapLUT::m_meana
double m_meana[ETAS1][PHIS1][PTS1]
Definition PtEndcapLUT.h:41
TrigL2MuonSA::PtEndcapLUT::m_meanb
double m_meanb[ETAS1][PHIS1][PTS1]
Definition PtEndcapLUT.h:42
type
TrigL2MuonSA::PtEndcapLUT::KeyType::KeyType
KeyType(int side, int charge, DataType type)
Definition PtEndcapLUT.h:55
TrigL2MuonSA::PtEndcapLUT::KeyType::operator<
bool operator<(const KeyType &other) const
Definition PtEndcapLUT.cxx:35
Generated on for ATLAS Offline Software by doxygen 1.14.0