ATLAS Offline Software
Loading...
Searching...
No Matches
SagittaRadiusEstimate.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 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
17TrigL2MuonSA::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
27void TrigL2MuonSA::SagittaRadiusEstimate::setMCFlag(bool 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
37StatusCode 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 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 else if ( i_station == 1 ) chamberID = xAOD::L2MuonParameters::Chamber::BarrelMiddle;
64 else if ( i_station == 2 ) chamberID = xAOD::L2MuonParameters::Chamber::BarrelOuter;
65 else 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
78 if ( count==2 ) {
79 y0 = 4230.; // radius of calorimeter.
80
81 if (superPoints[0]->R < ZERO_LIMIT) {
82 x2 = superPoints[1]->Z;
83 y2 = superPoints[1]->R;
84 x3 = superPoints[2]->Z;
85 y3 = superPoints[2]->R;
86 } else if (superPoints[1]->R < ZERO_LIMIT) {
87 x2 = superPoints[0]->Z;
88 y2 = superPoints[0]->R;
89 x3 = superPoints[2]->Z;
90 y3 = superPoints[2]->R;
91 } else if (superPoints[2]->R < ZERO_LIMIT) {
92 x2 = superPoints[0]->Z;
93 y2 = superPoints[0]->R;
94 x3 = superPoints[1]->Z;
95 y3 = superPoints[1]->R;
96 }
97
98 dx = x3 - x2;
99 dy = y3 - y2;
100
101 x0 = y0*x2/y2;
102
103 c3 = dy;
104 c2 = -y0*dx + 2.*(y2*x3-y3*x2);
105 c1 = -dy*(y2*y3-y0*y0)+ y3*x2*x2 - y2*x3*x3;
106 c0 = y0*x2*x3*dx + y0*x2*(y3-y0)*(y3-y0) - y0*x3*(y2-y0)*(y2-y0);
107 c22 = 2.*c2;
108 c33 = 3.*c3;
109
110 nit = 1;
111 while((nit++)<=nitmx&&std::abs(x0-xn)>=eps) {
112 xn = x0 - f(x0,c0,c1,c2,c3)/fp(x0,c33,c22,c1);
113 x0 = xn;
114 }
115 if (std::abs(xn)<ZERO_LIMIT) xn = ZERO_LIMIT;//To avoid divergence
116
117 x1 = xn;
118 y1 = y0;
119
120 if (superPoints[0]->R > ZERO_LIMIT ) {
121 theta = std::atan2(superPoints[0]->R,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 theta = std::atan2(y1,std::abs(x1));
125 signZ = (std::abs(x1) > ZERO_LIMIT)? x1/std::abs(x1): 1.;
126 }
127
128 trackPattern.etaMap = (-std::log(std::tan(theta/2.)))*signZ;
129 if (rpcFitResult.isSuccess ) {
130 one = (std::cos(rpcFitResult.phi)>0)? 1: -1;
131 } else {
132 one = (std::cos(p_roids->phi())>0)? 1: -1;
133 }
134 phi = std::atan2(trackPattern.phiMSDir*one,one);
135
136 if(phim>=M_PI+0.1) phim = phim - 2*M_PI;
137
138 if(phim>=0) trackPattern.phiMap = (phi>=0.)? phi - phim : phim -std::abs(phi);
139 else trackPattern.phiMap = phi - phim;
140
141 trackPattern.phiMS = phi;
142
143 c2 = x2 - x1;
144 c3 = x3 - x1;
145 e2 = y2 - y1;
146 e3 = y3 - y1;
147 c2q = c2*c2 + e2*e2;
148 c3q = c3*c3 + e3*e3;
149 d = c2*e3 - c3*e2;
150 da = -c2q*e3 + c3q*e2;
151 db = -c2*c3q + c3*c2q;
152 a = da/d;
153 b = db/d;
154
155 x0 = a/2.;
156 y0 = b/2.;
157 trackPattern.barrelRadius = std::sqrt(x0*x0 + y0*y0);
158 trackPattern.charge = -1;
159 if(a<=0.) trackPattern.charge = 1;
160
161 } else if (count==3) {
162
163 theta = std::atan2(superPoints[0]->R,std::abs(superPoints[0]->Z));
164 signZ = (std::abs(superPoints[0]->Z) > ZERO_LIMIT)? superPoints[0]->Z/std::abs(superPoints[0]->Z): 1.;
165
166 trackPattern.etaMap = (-std::log(std::tan(theta/2.)))*signZ;
167
168 if (rpcFitResult.isSuccess ) {
169 one = (std::cos(rpcFitResult.phi)>0)? 1: -1;
170 } else {
171 one = (std::cos(p_roids->phi())>0)? 1: -1;
172 }
173 phi = std::atan2(trackPattern.phiMSDir*one,one);
174 if(phim>=M_PI+0.1) phim = phim - 2*M_PI;
175
176 if(phim>=0) trackPattern.phiMap = (phi>=0.)? phi - phim : phim -std::abs(phi);
177 else trackPattern.phiMap = phi - phim;
178
179 trackPattern.phiMS = phi;
180
181 // Alignment correation to LargeSpecial
182 if ( trackPattern.s_address==1) {
183
184 if ( !m_alignmentBarrelLUT ) {
185 ATH_MSG_ERROR("Alignment correction service is not prepared");
186 return StatusCode::FAILURE;
187 }
188
189 double dZ = (*m_alignmentBarrelLUT)->GetDeltaZ(trackPattern.s_address,
190 trackPattern.etaMap,
191 trackPattern.phiMap,
192 trackPattern.phiMS,
193 superPoints[0]->R);
194 superPoints[1]->Z += 10*dZ;
195 }
196
197 m = ( superPoints[2]->Z - superPoints[0]->Z ) / ( superPoints[2]->R - superPoints[0]->R );
198
199 trackPattern.barrelSagitta = superPoints[1]->Z - superPoints[1]->R*m - superPoints[0]->Z + superPoints[0]->R*m;
200
201 cost = std::cos(std::atan(m));
202 x2 = superPoints[1]->R - superPoints[0]->R;
203 y2 = superPoints[1]->Z - superPoints[0]->Z;
204 x3 = superPoints[2]->R - superPoints[0]->R;
205 y3 = superPoints[2]->Z - superPoints[0]->Z;
206
207 tm = x2;
208 x2 = ( x2 + y2*m)*cost;
209 y2 = (-tm*m + y2 )*cost;
210
211 tm = x3;
212 x3 = ( x3 + y3*m)*cost;
213 y3 = (-tm*m + y3 )*cost;
214
215 x0 = x3/2.;
216 y0 = (y2*y2 + x2*x2 -x2*x3)/(2*y2);
217
218 trackPattern.barrelRadius = std::sqrt(x0*x0 + y0*y0);
219 trackPattern.charge = (trackPattern.barrelSagitta!=0)? -1.*trackPattern.barrelSagitta/std::abs(trackPattern.barrelSagitta): 0;
220
221 }
222 else if ( m_use_endcapInner == true && count == 1 && superPoints[3]->R > ZERO_LIMIT ) {
223 y0 = 4230.; // radius of calorimeter.
224
225 if (superPoints[0]->R > ZERO_LIMIT) {
226 x2 = superPoints[0]->Z; //BI
227 y2 = superPoints[0]->R;
228 x3 = superPoints[3]->Z; //EI
229 y3 = superPoints[3]->R;
230 } else if (superPoints[1]->R > ZERO_LIMIT) {
231 x2 = superPoints[3]->Z; //EI
232 y2 = superPoints[3]->R;
233 x3 = superPoints[1]->Z; //BM
234 y3 = superPoints[1]->R;
235 } else if (superPoints[2]->R > ZERO_LIMIT) {
236 x2 = superPoints[3]->Z; //EI
237 y2 = superPoints[3]->R;
238 x3 = superPoints[2]->Z; //BO
239 y3 = superPoints[2]->R;
240 }
241
242 dx = x3 - x2;
243 dy = y3 - y2;
244
245 x0 = y0*x2/y2;
246
247 c3 = dy;
248 c2 = -y0*dx + 2.*(y2*x3-y3*x2);
249 c1 = -dy*(y2*y3-y0*y0)+ y3*x2*x2 - y2*x3*x3;
250 c0 = y0*x2*x3*dx + y0*x2*(y3-y0)*(y3-y0) - y0*x3*(y2-y0)*(y2-y0);
251 c22 = 2.*c2;
252 c33 = 3.*c3;
253
254 nit = 1;
255 while((nit++)<=nitmx&&std::abs(x0-xn)>=eps) {
256 xn = x0 - f(x0,c0,c1,c2,c3)/fp(x0,c33,c22,c1);
257 x0 = xn;
258 }
259 if (std::abs(xn)<ZERO_LIMIT) xn = ZERO_LIMIT;//To avoid divergence
260
261 x1 = xn;
262 y1 = y0;
263
264 if (superPoints[0]->R > ZERO_LIMIT ) {
265 rad = superPoints[0]->R;
266 theta = std::atan2(rad,std::abs(superPoints[0]->Z));
267 signZ = (std::abs(superPoints[0]->Z) > ZERO_LIMIT)? superPoints[0]->Z/std::abs(superPoints[0]->Z): 1.;
268 } else {
269 rad = y1;
270 theta = std::atan2(rad,std::abs(x1));
271 signZ = (std::abs(x1) > ZERO_LIMIT)? x1/std::abs(x1): 1.;
272 }
273
274 trackPattern.etaMap = (-std::log(std::tan(theta/2.)))*signZ;
275 if (rpcFitResult.isSuccess ) {
276 // one = (std::abs(rpcFitResult.rpc1[0]) > 0.)? 1. * rpcFitResult.rpc1[0] / std::abs(rpcFitResult.rpc1[0]): 1.;
277 one = (std::cos(rpcFitResult.phi)>0)? 1: -1;
278 } else {
279 one = (std::cos(p_roids->phi())>0)? 1: -1;
280 }
281 phi = std::atan2(trackPattern.phiMSDir*one,one);
282
283 if(phim>=M_PI+0.1) phim = phim - 2*M_PI;
284
285 if(phim>=0) trackPattern.phiMap = (phi>=0.)? phi - phim : phim -std::abs(phi);
286 else trackPattern.phiMap = phi - phim;
287
288 trackPattern.phiMS = phi;
289
290 c2 = x2 - x1;
291 c3 = x3 - x1;
292 e2 = y2 - y1;
293 e3 = y3 - y1;
294 c2q = c2*c2 + e2*e2;
295 c3q = c3*c3 + e3*e3;
296 d = c2*e3 - c3*e2;
297 da = -c2q*e3 + c3q*e2;
298 db = -c2*c3q + c3*c2q;
299 a = da/d;
300 b = db/d;
301
302 x0 = -a/2. + x1;
303 y0 = -b/2. + y1;
304 double barrelRadius = std::sqrt(x0*x0 + y0*y0);
305 trackPattern.barrelRadius = barrelRadius;
306 trackPattern.charge = -1;
307 if(a<=0.) trackPattern.charge = 1;
308 }
309
310 ATH_MSG_DEBUG("... count/trackPattern.barrelSagitta/barrelRadius/charge/s_address/phi="
311 << count << " / " << trackPattern.barrelSagitta << "/" << trackPattern.barrelRadius << "/"
312 << trackPattern.charge << "/" << trackPattern.s_address << "/"
313 << trackPattern.phiMS);
314
315 return StatusCode::SUCCESS;
316}
317
318// --------------------------------------------------------------------------------
319// --------------------------------------------------------------------------------
320
M_PI
#define M_PI
Definition ActiveFraction.h:11
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
theta
Scalar theta() const
theta method
Definition AmgMatrixBasePlugin.h:75
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
ZERO_LIMIT
const float ZERO_LIMIT
Definition VP1TriggerHandleL2.cxx:37
AthAlgTool::AthAlgTool
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
Definition AthAlgTool.cxx:16
RoiDescriptor::phi
virtual double phi() const override final
Methods to retrieve data members.
Definition RoiDescriptor.h:100
TrigL2MuonSA::AlignmentBarrelLUTSvc::alignmentBarrelLUT
const ToolHandle< AlignmentBarrelLUT > * alignmentBarrelLUT(void) const
Definition AlignmentBarrelLUTSvc.h:35
TrigL2MuonSA::SagittaRadiusEstimate::m_alignmentBarrelLUT
const ToolHandle< AlignmentBarrelLUT > * m_alignmentBarrelLUT
Definition SagittaRadiusEstimate.h:47
TrigL2MuonSA::SagittaRadiusEstimate::setSagittaRadius
StatusCode setSagittaRadius(const TrigRoiDescriptor *p_roids, TrigL2MuonSA::RpcFitResult &rpcFitResult, TrigL2MuonSA::TrackPattern &trackPattern) const
Definition SagittaRadiusEstimate.cxx:37
TrigL2MuonSA::SagittaRadiusEstimate::f
float f(float x, float c0, float c1, float c2, float c3) const
Definition SagittaRadiusEstimate.h:54
TrigL2MuonSA::SagittaRadiusEstimate::SagittaRadiusEstimate
SagittaRadiusEstimate(const std::string &type, const std::string &name, const IInterface *parent)
Definition SagittaRadiusEstimate.cxx:17
TrigL2MuonSA::SagittaRadiusEstimate::setMCFlag
void setMCFlag(bool use_mcLUT, const AlignmentBarrelLUTSvc *alignmentBarrelLUTSvc)
Definition SagittaRadiusEstimate.cxx:27
TrigL2MuonSA::SagittaRadiusEstimate::fp
float fp(float x, float c33, float c22, float c1) const
Definition SagittaRadiusEstimate.h:59
TrigL2MuonSA::TrackPattern::superPoints
TrigL2MuonSA::SuperPoint superPoints[s_NCHAMBER]
Definition TrackData.h:60
TrigRoiDescriptor
nope - should be used for standalone also, perhaps need to protect the class def bits ifndef XAOD_ANA...
Definition TrigRoiDescriptor.h:56
cost
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)
Definition hcg.cxx:922
count
int count(std::string s, const std::string &regx)
count how many occurances of a regx are in a string
Definition hcg.cxx:146
xAOD::L2MuonParameters::BarrelInner
@ BarrelInner
Inner station in the barrel spectrometer.
Definition TrigMuonDefs.h:16
xAOD::L2MuonParameters::BarrelMiddle
@ BarrelMiddle
Middle station in the barrel spectrometer.
Definition TrigMuonDefs.h:17
xAOD::L2MuonParameters::BarrelOuter
@ BarrelOuter
Outer station in the barrel spectrometer.
Definition TrigMuonDefs.h:18
xAOD::L2MuonParameters::EndcapInner
@ EndcapInner
Inner station in the endcap spectrometer.
Definition TrigMuonDefs.h:19
type
Generated on for ATLAS Offline Software by doxygen 1.14.0