6 #include "CLHEP/Vector/TwoVector.h"
16 const std::string&
name,
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);
51 std::ostringstream
ss;
52 ss <<
"side=" << m_side <<
" charge=" << m_charge <<
" type=" << dt2s(
m_type);
61 std::ifstream ifs(lut_fileName.c_str());
64 return StatusCode::FAILURE;
77 if (
line.empty())
continue;
79 if (
line.substr(0, 5) ==
"side=") {
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;
89 ATH_MSG_ERROR(
"Invalid header line " << line_no <<
" in EndcapLUT file " << lut_fileName);
90 return StatusCode::FAILURE;
95 m_tables.insert(TableMap::value_type(
key,
table));
101 ATH_MSG_ERROR(
"Missing header line " << line_no <<
" in EndcapLUT file " << lut_fileName);
102 return StatusCode::FAILURE;
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;
117 return StatusCode::SUCCESS;
125 CLHEP::Hep2Vector mid1(z1, r1);
126 CLHEP::Hep2Vector mid2(z2, r2);
127 CLHEP::Hep2Vector
dir = mid2 - mid1;
129 if(std::abs(mid1.unit() *
dir.unit()) > 1.)
return 0;
131 double a = std::acos(mid1.unit() *
dir.unit());
141 float cr1 = 0.080/400;
145 if (std::abs(
x1)>=0.1) {
151 if(!(std::abs(
s2+999)<0.1)) {
157 double xm = (
x1+
x2)/2.;
158 double ym = (
y1+
y2)/2.;
160 double c2_1 = -
x1-A1*
y1;
161 double c2_2 = -
x2-A2*
y2;
164 double xR1 = -A1*yR1-c2_1;
165 double xR2 = -A2*yR2-c2_2;
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)));
183 <<
",type=" << dt2s(
type) <<
",iEta=" <<
iEta
184 <<
",iPhi=" <<
iPhi <<
",value=" <<
value <<
")");
189 if (iEta < 0 || iEta >= ETAS || iPhi < 0 || iPhi >= PHISEE) {
191 <<
", " <<
iEta <<
", " <<
iPhi <<
") Invalid indices");
197 if (
it == m_tables.end()) {
199 <<
", " <<
iEta <<
", " <<
iPhi <<
") Invalid key");
208 const double PT_MAX = 500;
227 if( ptinv < (1/PT_MAX) )
pt = PT_MAX;
228 else {
pt = 1/ptinv; }
239 double ptinv = (-
b +
pow(judge,0.5)) /
c / 2;
247 if( ptinv < (1/PT_MAX) ) {
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);
285 return (
"alphapol2");
break;
288 return (
"betapol2");
break;
291 return (
"tgcalphapol2");
break;
294 return (
"invradiuspol2");
break;
297 return (
"cscpol2");
break;
300 return (
"invalid");
break;
303 return (
"invalid");
break;
316 <<
"Alpha pT=" << ApT
323 if (iEta < 0 || iEta >= ETAS || iPhi < 0 || iPhi >= PHIS) {
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];
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];
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;
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;
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;
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.;
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;
393 if (!ifssigma.is_open()) {
394 ATH_MSG_ERROR(
"Cannot open EndcapLUT Sigma file " << lut_sigma);
395 return StatusCode::FAILURE;
400 for(
int ei=0; ei <
ETAS1; ei++){
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.;
413 while (!ifsmean.eof()) {
414 getline(ifsmean,
line);
415 if (
line.empty())
continue;
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;
424 m_meana[
iEta][
iPhi][iNP] = tmp_par1;
425 m_meanb[
iEta][
iPhi][iNP] = tmp_par2;
426 m_meanc[
iEta][
iPhi][iNP] = tmp_par3;
430 while (!ifssigma.eof()) {
431 getline(ifssigma,
line2);
432 if (
line2.empty())
continue;
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;
441 m_sigmaa[
iEta][
iPhi][iNP] = tmp_par1;
442 m_sigmab[
iEta][
iPhi][iNP] = tmp_par2;
443 m_sigmac[
iEta][
iPhi][iNP] = tmp_par3;
447 return StatusCode::SUCCESS;