ATLAS Offline Software
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 
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(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 
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 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 
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:215
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
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
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:142
CaloCondBlobAlgs_fillNoiseFromASCII.db
db
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:42
plotBeamSpotCompare.x2
x2
Definition: plotBeamSpotCompare.py:217
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:75
extractSporadic.c1
c1
Definition: extractSporadic.py:133
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:129
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
XMLtoHeader.count
count
Definition: XMLtoHeader.py:84
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
TrigL2MuonSA::SuperPoint::R
float R
Definition: SuperPointData.h:55
drawFromPickle.atan
atan
Definition: drawFromPickle.py:36
TrigL2MuonSA::TrackPattern
Definition: TrackData.h:16
hotSpotInTAG.c0
c0
Definition: hotSpotInTAG.py:191
TrigL2MuonSA::SuperPoint::Phim
float Phim
Definition: SuperPointData.h:57
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:140
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
python.SystemOfUnits.rad
float rad
Definition: SystemOfUnits.py:126
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:240
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
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
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
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:56
TrigL2MuonSA::SuperPoint
Definition: SuperPointData.h:38
AthAlgTool
Definition: AthAlgTool.h:26
xAOD::L2MuonParameters::BarrelMiddle
@ BarrelMiddle
Middle station in the barrel spectrometer.
Definition: TrigMuonDefs.h:17
TrigL2MuonSA::SagittaRadiusEstimate::setMCFlag
void setMCFlag(bool use_mcLUT, const AlignmentBarrelLUTSvc *alignmentBarrelLUTSvc)
Definition: SagittaRadiusEstimate.cxx:27
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
Generated on Mon Sep 29 2025 21:18:03 for ATLAS Offline Software by   doxygen 1.8.18