ATLAS Offline Software
Loading...
Searching...
No Matches
SagittaRadiusEstimate.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5
6
7#include "SagittaRadiusEstimate.h"
8
9#include "xAODTrigMuon/TrigMuonDefs.h"
10
11#include "AthenaBaseComps/AthMsgStreamMacros.h"
12#include <cmath>
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 if (chamberID < 0)[[unlikely]]{
67 ATH_MSG_ERROR("Chamber ID is ill defined.");
68 return StatusCode::FAILURE;
69 }
70 superPoints[i_station] = &(trackPattern.superPoints[chamberID]);
71
72 if (superPoints[i_station]->R > ZERO_LIMIT) {
73 count++;
74 if ( i_station != 3 ){
75 phim = superPoints[i_station]->Phim;
76 }
77 }
78 }
79
80 if ( superPoints[3] -> R >ZERO_LIMIT ) count--; // Not use Endcap Inner
81
82 if ( count==2 ) {
83 y0 = 4230.; // radius of calorimeter.
84
85 if (superPoints[0]->R < ZERO_LIMIT) {
86 x2 = superPoints[1]->Z;
87 y2 = superPoints[1]->R;
88 x3 = superPoints[2]->Z;
89 y3 = superPoints[2]->R;
90 } else if (superPoints[1]->R < ZERO_LIMIT) {
91 x2 = superPoints[0]->Z;
92 y2 = superPoints[0]->R;
93 x3 = superPoints[2]->Z;
94 y3 = superPoints[2]->R;
95 } else if (superPoints[2]->R < ZERO_LIMIT) {
96 x2 = superPoints[0]->Z;
97 y2 = superPoints[0]->R;
98 x3 = superPoints[1]->Z;
99 y3 = superPoints[1]->R;
100 }
101
102 dx = x3 - x2;
103 dy = y3 - y2;
104 if (y2 == 0.)[[unlikely]]{
105 throw std::runtime_error("y2 is zero in setSagittaRadius");
106 }
107 x0 = y0*x2/y2;
108
109 c3 = dy;
110 c2 = -y0*dx + 2.*(y2*x3-y3*x2);
111 c1 = -dy*(y2*y3-y0*y0)+ y3*x2*x2 - y2*x3*x3;
112 c0 = y0*x2*x3*dx + y0*x2*(y3-y0)*(y3-y0) - y0*x3*(y2-y0)*(y2-y0);
113 c22 = 2.*c2;
114 c33 = 3.*c3;
115
116 nit = 1;
117 while((nit++)<=nitmx&&std::abs(x0-xn)>=eps) {
118 xn = x0 - f(x0,c0,c1,c2,c3)/fp(x0,c33,c22,c1);
119 x0 = xn;
120 }
121 if (std::abs(xn)<ZERO_LIMIT) xn = ZERO_LIMIT;//To avoid divergence
122
123 x1 = xn;
124 y1 = y0;
125
126 if (superPoints[0]->R > ZERO_LIMIT ) {
127 theta = std::atan2(superPoints[0]->R,std::abs(superPoints[0]->Z));
128 signZ = (std::abs(superPoints[0]->Z) > ZERO_LIMIT)? superPoints[0]->Z/std::abs(superPoints[0]->Z): 1.;
129 } else {
130 theta = std::atan2(y1,std::abs(x1));
131 signZ = (std::abs(x1) > ZERO_LIMIT)? x1/std::abs(x1): 1.;
132 }
133
134 trackPattern.etaMap = (-std::log(std::tan(theta/2.)))*signZ;
135 if (rpcFitResult.isSuccess ) {
136 one = (std::cos(rpcFitResult.phi)>0)? 1: -1;
137 } else {
138 one = (std::cos(p_roids->phi())>0)? 1: -1;
139 }
140 phi = std::atan2(trackPattern.phiMSDir*one,one);
141
142 if(phim>=M_PI+0.1) phim = phim - 2*M_PI;
143
144 if(phim>=0) trackPattern.phiMap = (phi>=0.)? phi - phim : phim -std::abs(phi);
145 else trackPattern.phiMap = phi - phim;
146
147 trackPattern.phiMS = phi;
148
149 c2 = x2 - x1;
150 c3 = x3 - x1;
151 e2 = y2 - y1;
152 e3 = y3 - y1;
153 c2q = c2*c2 + e2*e2;
154 c3q = c3*c3 + e3*e3;
155 d = c2*e3 - c3*e2;
156 da = -c2q*e3 + c3q*e2;
157 db = -c2*c3q + c3*c2q;
158 a = da/d;
159 b = db/d;
160
161 x0 = a/2.;
162 y0 = b/2.;
163 trackPattern.barrelRadius = std::sqrt(x0*x0 + y0*y0);
164 trackPattern.charge = -1;
165 if(a<=0.) trackPattern.charge = 1;
166
167 } else if (count==3) {
168
169 theta = std::atan2(superPoints[0]->R,std::abs(superPoints[0]->Z));
170 signZ = (std::abs(superPoints[0]->Z) > ZERO_LIMIT)? superPoints[0]->Z/std::abs(superPoints[0]->Z): 1.;
171
172 trackPattern.etaMap = (-std::log(std::tan(theta/2.)))*signZ;
173
174 if (rpcFitResult.isSuccess ) {
175 one = (std::cos(rpcFitResult.phi)>0)? 1: -1;
176 } else {
177 one = (std::cos(p_roids->phi())>0)? 1: -1;
178 }
179 phi = std::atan2(trackPattern.phiMSDir*one,one);
180 if(phim>=M_PI+0.1) phim = phim - 2*M_PI;
181
182 if(phim>=0) trackPattern.phiMap = (phi>=0.)? phi - phim : phim -std::abs(phi);
183 else trackPattern.phiMap = phi - phim;
184
185 trackPattern.phiMS = phi;
186
187 // Alignment correation to LargeSpecial
188 if ( trackPattern.s_address==1) {
189
190 if ( !m_alignmentBarrelLUT ) {
191 ATH_MSG_ERROR("Alignment correction service is not prepared");
192 return StatusCode::FAILURE;
193 }
194
195 double dZ = (*m_alignmentBarrelLUT)->GetDeltaZ(trackPattern.s_address,
196 trackPattern.etaMap,
197 trackPattern.phiMap,
198 trackPattern.phiMS,
199 superPoints[0]->R);
200 superPoints[1]->Z += 10*dZ;
201 }
202
203 m = ( superPoints[2]->Z - superPoints[0]->Z ) / ( superPoints[2]->R - superPoints[0]->R );
204
205 trackPattern.barrelSagitta = superPoints[1]->Z - superPoints[1]->R*m - superPoints[0]->Z + superPoints[0]->R*m;
206
207 cost = std::cos(std::atan(m));
208 x2 = superPoints[1]->R - superPoints[0]->R;
209 y2 = superPoints[1]->Z - superPoints[0]->Z;
210 x3 = superPoints[2]->R - superPoints[0]->R;
211 y3 = superPoints[2]->Z - superPoints[0]->Z;
212
213 tm = x2;
214 x2 = ( x2 + y2*m)*cost;
215 y2 = (-tm*m + y2 )*cost;
216
217 tm = x3;
218 x3 = ( x3 + y3*m)*cost;
219 y3 = (-tm*m + y3 )*cost;
220
221 x0 = x3/2.;
222 y0 = (y2*y2 + x2*x2 -x2*x3)/(2*y2);
223
224 trackPattern.barrelRadius = std::sqrt(x0*x0 + y0*y0);
225 trackPattern.charge = (trackPattern.barrelSagitta!=0)? -1.*trackPattern.barrelSagitta/std::abs(trackPattern.barrelSagitta): 0;
226
227 }
228 else 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,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,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
M_PI
#define M_PI
Definition ActiveFraction.h:14
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
theta
Scalar theta() const
theta method
Definition AmgMatrixBasePlugin.h:75
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x,...)
Definition AthMsgStreamMacros.h:43
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x,...)
Definition AthMsgStreamMacros.h:47
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
ZERO_LIMIT
static const double ZERO_LIMIT
Definition MuonMatchingTool.cxx:26
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:101
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:48
cost
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)
Definition hcg.cxx:926
count
int count(std::string s, const std::string &regx)
count how many occurances of a regx are in a string
Definition hcg.cxx:148
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
unlikely
#define unlikely(x)
Definition pythonic_coracool.cxx:9
type
Generated on for ATLAS Offline Software by doxygen 1.17.0