ATLAS Offline Software
SagittaRadiusEstimate.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 <cmath>
6 
7 #include "SagittaRadiusEstimate.h"
8 
9 #include "xAODTrigMuon/TrigMuonDefs.h"
10 
11 #include "AthenaBaseComps/AthMsgStreamMacros.h"
12 
13 
14 // --------------------------------------------------------------------------------
15 // --------------------------------------------------------------------------------
16 
17 TrigL2MuonSA::SagittaRadiusEstimate::SagittaRadiusEstimate(const std::string& type,
18  const std::string& name,
19  const IInterface* parent):
20  AthAlgTool(type, name, parent)
21 {
22 }
23 
24 // --------------------------------------------------------------------------------
25 // --------------------------------------------------------------------------------
26 
27 void TrigL2MuonSA::SagittaRadiusEstimate::setMCFlag(const BooleanProperty& use_mcLUT,
28  const AlignmentBarrelLUTSvc* alignmentBarrelLUTSvc)
29 {
30  m_use_mcLUT = use_mcLUT;
31  if ( alignmentBarrelLUTSvc ) m_alignmentBarrelLUT = alignmentBarrelLUTSvc->alignmentBarrelLUT();
32 }
33 
34 // --------------------------------------------------------------------------------
35 // --------------------------------------------------------------------------------
36 
37 StatusCode TrigL2MuonSA::SagittaRadiusEstimate::setSagittaRadius(const TrigRoiDescriptor* p_roids,
38  TrigL2MuonSA::RpcFitResult& rpcFitResult,
39  TrigL2MuonSA::TrackPattern& trackPattern) const
40 {
41  const int MAX_STATION = 4;
42  const float ZERO_LIMIT = 1e-5;
43 
44  int nit;
45  const int nitmx=10;
46  int count=0;
47  double a3,theta,rad,phi,one,phim=0,signZ;
48 
49  double c0,c1,c2,c3,c22,c33,e2,e3,c2q,c3q,d,da,db,a,b,dx,dy;
50  double m = 0.;
51  double cost = 0.;
52  double x0 = 0., y0 = 0., x1 = 0., y1 = 0., x2 = 0., y2 = 0., x3 = 0., y3 = 0.;
53  double tm = 0.;
54  double xn = 0.;
55  const double eps = 0.005;
56 
57  TrigL2MuonSA::SuperPoint* superPoints[4];
58 
59  for (int i_station=0; i_station<MAX_STATION; i_station++) {
60 
61  int chamberID = -1;
62  if ( i_station == 0 ) chamberID = xAOD::L2MuonParameters::Chamber::BarrelInner;
63  if ( i_station == 1 ) chamberID = xAOD::L2MuonParameters::Chamber::BarrelMiddle;
64  if ( i_station == 2 ) chamberID = xAOD::L2MuonParameters::Chamber::BarrelOuter;
65  if ( i_station == 3 ) chamberID = xAOD::L2MuonParameters::Chamber::EndcapInner;
66  superPoints[i_station] = &(trackPattern.superPoints[chamberID]);
67 
68  if (superPoints[i_station]->R > ZERO_LIMIT) {
69  count++;
70  if ( i_station != 3 ){
71  phim = superPoints[i_station]->Phim;
72  }
73  }
74  }
75 
76  if ( superPoints[3] -> R >ZERO_LIMIT ) count--; // Not use Endcap Inner
77  if ( count==2 ) {
78  y0 = 4230.; // radius of calorimeter.
79 
80  if (superPoints[0]->R < ZERO_LIMIT) {
81  x2 = superPoints[1]->Z;
82  y2 = superPoints[1]->R;
83  x3 = superPoints[2]->Z;
84  y3 = superPoints[2]->R;
85  } else if (superPoints[1]->R < ZERO_LIMIT) {
86  x2 = superPoints[0]->Z;
87  y2 = superPoints[0]->R;
88  x3 = superPoints[2]->Z;
89  y3 = superPoints[2]->R;
90  } else if (superPoints[2]->R < ZERO_LIMIT) {
91  x2 = superPoints[0]->Z;
92  y2 = superPoints[0]->R;
93  x3 = superPoints[1]->Z;
94  y3 = superPoints[1]->R;
95  }
96 
97  dx = x3 - x2;
98  dy = y3 - y2;
99 
100  x0 = y0*x2/y2;
101 
102  c3 = dy;
103  c2 = -y0*dx + 2.*(y2*x3-y3*x2);
104  c1 = -dy*(y2*y3-y0*y0)+ y3*x2*x2 - y2*x3*x3;
105  c0 = y0*x2*x3*dx + y0*x2*(y3-y0)*(y3-y0) - y0*x3*(y2-y0)*(y2-y0);
106  c22 = 2.*c2;
107  c33 = 3.*c3;
108 
109  nit = 1;
110  while((nit++)<=nitmx&&std::abs(x0-xn)>=eps) {
111  xn = x0 - f(x0,c0,c1,c2,c3)/fp(x0,c33,c22,c1);
112  x0 = xn;
113  }
114  if (std::abs(xn)<ZERO_LIMIT) xn = ZERO_LIMIT;//To avoid divergence
115 
116  x1 = xn;
117  y1 = y0;
118 
119  if (superPoints[0]->R > ZERO_LIMIT ) {
120  rad = superPoints[0]->R;
121  theta = std::atan2(rad,(double)std::abs(superPoints[0]->Z));
122  signZ = (std::abs(superPoints[0]->Z) > ZERO_LIMIT)? superPoints[0]->Z/std::abs(superPoints[0]->Z): 1.;
123  } else {
124  rad = y1;
125  theta = std::atan2(rad,(double)std::abs(x1));
126  signZ = (std::abs(x1) > ZERO_LIMIT)? x1/std::abs(x1): 1.;
127  }
128 
129  trackPattern.etaMap = (-std::log(std::tan(theta/2.)))*signZ;
130  if (rpcFitResult.isSuccess ) {
131  // one = (std::abs(rpcFitResult.rpc1[0]) > 0.)? 1. * rpcFitResult.rpc1[0] / std::abs(rpcFitResult.rpc1[0]): 1.;
132  one = (std::cos(rpcFitResult.phi)>0)? 1: -1;
133  } else {
134  one = (std::cos(p_roids->phi())>0)? 1: -1;
135  }
136  phi = std::atan2(trackPattern.phiMSDir*one,one);
137 
138  if(phim>=M_PI+0.1) phim = phim - 2*M_PI;
139 
140  if(phim>=0) trackPattern.phiMap = (phi>=0.)? phi - phim : phim -std::abs(phi);
141  else trackPattern.phiMap = phi - phim;
142 
143  trackPattern.phiMS = phi;
144 
145  c2 = x2 - x1;
146  c3 = x3 - x1;
147  e2 = y2 - y1;
148  e3 = y3 - y1;
149  c2q = c2*c2 + e2*e2;
150  c3q = c3*c3 + e3*e3;
151  d = c2*e3 - c3*e2;
152  da = -c2q*e3 + c3q*e2;
153  db = -c2*c3q + c3*c2q;
154  a = da/d;
155  b = db/d;
156 
157  x0 = -a/2. + x1;
158  y0 = -b/2. + y1;
159  trackPattern.barrelRadius = std::sqrt(x0*x0 + y0*y0);
160  trackPattern.charge = -1;
161  if(a<=0.) trackPattern.charge = 1;
162 
163  } else if (count==3) {
164 
165  rad = superPoints[0]->R;
166  theta = std::atan2(rad,(double)std::abs(superPoints[0]->Z));
167  signZ = (std::abs(superPoints[0]->Z) > ZERO_LIMIT)? superPoints[0]->Z/std::abs(superPoints[0]->Z): 1.;
168 
169  trackPattern.etaMap = (-std::log(std::tan(theta/2.)))*signZ;
170 
171  if (rpcFitResult.isSuccess ) {
172  // one = (std::abs(rpcFitResult.rpc1[0]) > 0.)? 1. * rpcFitResult.rpc1[0] / std::abs(rpcFitResult.rpc1[0]): 1;
173  one = (std::cos(rpcFitResult.phi)>0)? 1: -1;
174  } else {
175  one = (std::cos(p_roids->phi())>0)? 1: -1;
176  }
177  phi = std::atan2(trackPattern.phiMSDir*one,one);
178  if(phim>=M_PI+0.1) phim = phim - 2*M_PI;
179 
180  if(phim>=0) trackPattern.phiMap = (phi>=0.)? phi - phim : phim -std::abs(phi);
181  else trackPattern.phiMap = phi - phim;
182 
183  trackPattern.phiMS = phi;
184 
185  // Alignment correation to LargeSpecial
186  if ( trackPattern.s_address==1) {
187 
188  if ( !m_alignmentBarrelLUT ) {
189  ATH_MSG_ERROR("Alignment correction service is not prepared");
190  return StatusCode::FAILURE;
191  }
192 
193  double dZ = (*m_alignmentBarrelLUT)->GetDeltaZ(trackPattern.s_address,
194  trackPattern.etaMap,
195  trackPattern.phiMap,
196  trackPattern.phiMS,
197  superPoints[0]->R);
198  superPoints[1]->Z += 10*dZ;
199  }
200 
201  a3 = ( superPoints[2]->Z - superPoints[0]->Z ) / ( superPoints[2]->R - superPoints[0]->R );
202 
203  trackPattern.barrelSagitta = superPoints[1]->Z - superPoints[1]->R*a3 - superPoints[0]->Z + superPoints[0]->R*a3;
204 
205  m = a3;
206  cost = std::cos(std::atan(m));
207  x2 = superPoints[1]->R - superPoints[0]->R;
208  y2 = superPoints[1]->Z - superPoints[0]->Z;
209  x3 = superPoints[2]->R - superPoints[0]->R;
210  y3 = superPoints[2]->Z - superPoints[0]->Z;
211 
212  tm = x2;
213  x2 = ( x2 + y2*m)*cost;
214  y2 = (-tm*m + y2 )*cost;
215 
216  tm = x3;
217  x3 = ( x3 + y3*m)*cost;
218  y3 = (-tm*m + y3 )*cost;
219 
220  x0 = x3/2.;
221  y0 = (y2*y2 + x2*x2 -x2*x3)/(2*y2);
222 
223  trackPattern.barrelRadius = std::sqrt(x0*x0 + y0*y0);
224  trackPattern.charge = (trackPattern.barrelSagitta!=0)? -1.*trackPattern.barrelSagitta/std::abs(trackPattern.barrelSagitta): 0;
225 
226  }
227 
228  if ( m_use_endcapInner == true && count == 1 && superPoints[3]->R > ZERO_LIMIT ) {
229  y0 = 4230.; // radius of calorimeter.
230 
231  if (superPoints[0]->R > ZERO_LIMIT) {
232  x2 = superPoints[0]->Z; //BI
233  y2 = superPoints[0]->R;
234  x3 = superPoints[3]->Z; //EI
235  y3 = superPoints[3]->R;
236  } else if (superPoints[1]->R > ZERO_LIMIT) {
237  x2 = superPoints[3]->Z; //EI
238  y2 = superPoints[3]->R;
239  x3 = superPoints[1]->Z; //BM
240  y3 = superPoints[1]->R;
241  } else if (superPoints[2]->R > ZERO_LIMIT) {
242  x2 = superPoints[3]->Z; //EI
243  y2 = superPoints[3]->R;
244  x3 = superPoints[2]->Z; //BO
245  y3 = superPoints[2]->R;
246  }
247 
248  dx = x3 - x2;
249  dy = y3 - y2;
250 
251  x0 = y0*x2/y2;
252 
253  c3 = dy;
254  c2 = -y0*dx + 2.*(y2*x3-y3*x2);
255  c1 = -dy*(y2*y3-y0*y0)+ y3*x2*x2 - y2*x3*x3;
256  c0 = y0*x2*x3*dx + y0*x2*(y3-y0)*(y3-y0) - y0*x3*(y2-y0)*(y2-y0);
257  c22 = 2.*c2;
258  c33 = 3.*c3;
259 
260  nit = 1;
261  while((nit++)<=nitmx&&std::abs(x0-xn)>=eps) {
262  xn = x0 - f(x0,c0,c1,c2,c3)/fp(x0,c33,c22,c1);
263  x0 = xn;
264  }
265  if (std::abs(xn)<ZERO_LIMIT) xn = ZERO_LIMIT;//To avoid divergence
266 
267  x1 = xn;
268  y1 = y0;
269 
270  if (superPoints[0]->R > ZERO_LIMIT ) {
271  rad = superPoints[0]->R;
272  theta = std::atan2(rad,(double)std::abs(superPoints[0]->Z));
273  signZ = (std::abs(superPoints[0]->Z) > ZERO_LIMIT)? superPoints[0]->Z/std::abs(superPoints[0]->Z): 1.;
274  } else {
275  rad = y1;
276  theta = std::atan2(rad,(double)std::abs(x1));
277  signZ = (std::abs(x1) > ZERO_LIMIT)? x1/std::abs(x1): 1.;
278  }
279 
280  trackPattern.etaMap = (-std::log(std::tan(theta/2.)))*signZ;
281  if (rpcFitResult.isSuccess ) {
282  // one = (std::abs(rpcFitResult.rpc1[0]) > 0.)? 1. * rpcFitResult.rpc1[0] / std::abs(rpcFitResult.rpc1[0]): 1.;
283  one = (std::cos(rpcFitResult.phi)>0)? 1: -1;
284  } else {
285  one = (std::cos(p_roids->phi())>0)? 1: -1;
286  }
287  phi = std::atan2(trackPattern.phiMSDir*one,one);
288 
289  if(phim>=M_PI+0.1) phim = phim - 2*M_PI;
290 
291  if(phim>=0) trackPattern.phiMap = (phi>=0.)? phi - phim : phim -std::abs(phi);
292  else trackPattern.phiMap = phi - phim;
293 
294  trackPattern.phiMS = phi;
295 
296  c2 = x2 - x1;
297  c3 = x3 - x1;
298  e2 = y2 - y1;
299  e3 = y3 - y1;
300  c2q = c2*c2 + e2*e2;
301  c3q = c3*c3 + e3*e3;
302  d = c2*e3 - c3*e2;
303  da = -c2q*e3 + c3q*e2;
304  db = -c2*c3q + c3*c2q;
305  a = da/d;
306  b = db/d;
307 
308  x0 = -a/2. + x1;
309  y0 = -b/2. + y1;
310  double barrelRadius = std::sqrt(x0*x0 + y0*y0);
311  trackPattern.barrelRadius = barrelRadius;
312  trackPattern.charge = -1;
313  if(a<=0.) trackPattern.charge = 1;
314  }
315 
316  ATH_MSG_DEBUG("... count/trackPattern.barrelSagitta/barrelRadius/charge/s_address/phi="
317  << count << " / " << trackPattern.barrelSagitta << "/" << trackPattern.barrelRadius << "/"
318  << trackPattern.charge << "/" << trackPattern.s_address << "/"
319  << trackPattern.phiMS);
320 
321  return StatusCode::SUCCESS;
322 }
323 
324 // --------------------------------------------------------------------------------
325 // --------------------------------------------------------------------------------
326 
TrigL2MuonSA::TrackPattern::phiMap
double phiMap
Definition: TrackData.h:78
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
cost
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)
Definition: hcg.cxx:921
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:216
TrigL2MuonSA::SagittaRadiusEstimate::SagittaRadiusEstimate
SagittaRadiusEstimate(const std::string &type, const std::string &name, const IInterface *parent)
Definition: SagittaRadiusEstimate.cxx:17
TrigL2MuonSA::TrackPattern::superPoints
TrigL2MuonSA::SuperPoint superPoints[s_NCHAMBER]
Definition: TrackData.h:60
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
AthMsgStreamMacros.h
SagittaRadiusEstimate.h
TrigL2MuonSA::RpcFitResult::phi
double phi
Definition: RpcFitResult.h:44
TrigL2MuonSA::TrackPattern::barrelRadius
double barrelRadius
Definition: TrackData.h:84
TrigL2MuonSA::RpcFitResult
Definition: RpcFitResult.h:14
hist_file_dump.d
d
Definition: hist_file_dump.py:137
CaloCondBlobAlgs_fillNoiseFromASCII.db
db
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:43
plotBeamSpotCompare.x2
x2
Definition: plotBeamSpotCompare.py:218
extractSporadic.c1
c1
Definition: extractSporadic.py:134
TrigL2MuonSA::TrackPattern::s_address
int s_address
Definition: TrackData.h:72
M_PI
#define M_PI
Definition: ActiveFraction.h:11
xAOD::L2MuonParameters::BarrelInner
@ BarrelInner
Inner station in the barrel spectrometer.
Definition: TrigMuonDefs.h:16
Trk::one
@ one
Definition: TrkDetDescr/TrkSurfaces/TrkSurfaces/RealQuadraticEquation.h:22
TrigRoiDescriptor
nope - should be used for standalone also, perhaps need to protect the class def bits #ifndef XAOD_AN...
Definition: TrigRoiDescriptor.h:56
xAOD::barrelRadius
setSAddress setEtaMS setDirPhiMS setDirZMS barrelRadius
Definition: L2StandAloneMuon_v1.cxx:128
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
XMLtoHeader.count
count
Definition: XMLtoHeader.py:85
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
TrigL2MuonSA::SuperPoint::R
float R
Definition: SuperPointData.h:102
drawFromPickle.atan
atan
Definition: drawFromPickle.py:36
TrigL2MuonSA::TrackPattern
Definition: TrackData.h:16
hotSpotInTAG.c0
c0
Definition: hotSpotInTAG.py:192
TrigL2MuonSA::SuperPoint::Phim
float Phim
Definition: SuperPointData.h:104
TrigL2MuonSA::AlignmentBarrelLUTSvc
Definition: AlignmentBarrelLUTSvc.h:18
python.DataFormatRates.c3
c3
Definition: DataFormatRates.py:127
TrigL2MuonSA::TrackPattern::phiMS
double phiMS
Definition: TrackData.h:74
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
CheckAppliedSFs.e3
e3
Definition: CheckAppliedSFs.py:264
TrigL2MuonSA::SagittaRadiusEstimate::setSagittaRadius
StatusCode setSagittaRadius(const TrigRoiDescriptor *p_roids, TrigL2MuonSA::RpcFitResult &rpcFitResult, TrigL2MuonSA::TrackPattern &trackPattern) const
Definition: SagittaRadiusEstimate.cxx:37
trigmenu_modify_prescale_json.fp
fp
Definition: trigmenu_modify_prescale_json.py:53
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TrigL2MuonSA::TrackPattern::barrelSagitta
double barrelSagitta
Definition: TrackData.h:85
makeTRTBarrelCans.y2
tuple y2
Definition: makeTRTBarrelCans.py:18
test_pyathena.parent
parent
Definition: test_pyathena.py:15
hist_file_dump.f
f
Definition: hist_file_dump.py:135
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
ZERO_LIMIT
const float ZERO_LIMIT
Definition: VP1TriggerHandleL2.cxx:37
TrigL2MuonSA::TrackPattern::charge
double charge
Definition: TrackData.h:63
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
TrigMuonDefs.h
makeTRTBarrelCans.dy
tuple dy
Definition: makeTRTBarrelCans.py:21
TrigL2MuonSA::RpcFitResult::isSuccess
bool isSuccess
Definition: RpcFitResult.h:38
python.DataFormatRates.c2
c2
Definition: DataFormatRates.py:123
TrigL2MuonSA::SagittaRadiusEstimate::setMCFlag
void setMCFlag(const BooleanProperty &use_mcLUT, const AlignmentBarrelLUTSvc *alignmentBarrelLUTSvc)
Definition: SagittaRadiusEstimate.cxx:27
a
TList * a
Definition: liststreamerinfos.cxx:10
egammaEnergyPositionAllSamples::e2
double e2(const xAOD::CaloCluster &cluster)
return the uncorrected cluster energy in 2nd sampling
RoiDescriptor::phi
virtual double phi() const override final
Methods to retrieve data members.
Definition: RoiDescriptor.h:100
TrigL2MuonSA::TrackPattern::phiMSDir
double phiMSDir
Definition: TrackData.h:75
xAOD::L2MuonParameters::EndcapInner
@ EndcapInner
Inner station in the endcap spectrometer.
Definition: TrigMuonDefs.h:19
TrigL2MuonSA::TrackPattern::etaMap
double etaMap
Definition: TrackData.h:77
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
xAOD::L2MuonParameters::BarrelOuter
@ BarrelOuter
Outer station in the barrel spectrometer.
Definition: TrigMuonDefs.h:18
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
TrigL2MuonSA::AlignmentBarrelLUTSvc::alignmentBarrelLUT
const ToolHandle< AlignmentBarrelLUT > * alignmentBarrelLUT(void) const
Definition: AlignmentBarrelLUTSvc.h:35
TrigL2MuonSA::SuperPoint::Z
float Z
Definition: SuperPointData.h:103
TrigL2MuonSA::SuperPoint
Definition: SuperPointData.h:74
AthAlgTool
Definition: AthAlgTool.h:26
xAOD::L2MuonParameters::BarrelMiddle
@ BarrelMiddle
Middle station in the barrel spectrometer.
Definition: TrigMuonDefs.h:17
python.SystemOfUnits.rad
int rad
Definition: SystemOfUnits.py:111
Generated on Mon Dec 23 2024 21:17:40 for ATLAS Offline Software by   doxygen 1.8.18