ATLAS Offline Software
DiscBounds.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // DiscBounds.cxx, (c) ATLAS Detector Software
8 
9 // Trk
10 #include "TrkSurfaces/DiscBounds.h"
11 // Gaudi
12 #include "GaudiKernel/MsgStream.h"
13 // STD
14 #include <iomanip>
15 #include <iostream>
16 
17 Trk::DiscBounds::DiscBounds()
18  : m_boundValues(DiscBounds::bv_length, 0.)
19 {}
20 
21 Trk::DiscBounds::DiscBounds(double minrad, double maxrad, double hphisec)
22  : m_boundValues(DiscBounds::bv_length, 0.)
23 {
24  m_boundValues[DiscBounds::bv_rMin] = minrad;
25  m_boundValues[DiscBounds::bv_rMax] = maxrad;
26  m_boundValues[DiscBounds::bv_averagePhi] = 0.;
27  m_boundValues[DiscBounds::bv_halfPhiSector] = hphisec;
28  if (m_boundValues[DiscBounds::bv_rMin] > m_boundValues[DiscBounds::bv_rMax])
29  swap(m_boundValues[DiscBounds::bv_rMin], m_boundValues[DiscBounds::bv_rMax]);
30 }
31 
32 Trk::DiscBounds::DiscBounds(double minrad, double maxrad, double avephi, double hphisec)
33  : m_boundValues(DiscBounds::bv_length, 0.)
34 {
35  m_boundValues[DiscBounds::bv_rMin] = minrad;
36  m_boundValues[DiscBounds::bv_rMax] = maxrad;
37  m_boundValues[DiscBounds::bv_averagePhi] = avephi;
38  m_boundValues[DiscBounds::bv_halfPhiSector] = hphisec;
39  if (m_boundValues[DiscBounds::bv_rMin] > m_boundValues[DiscBounds::bv_rMax])
40  swap(m_boundValues[DiscBounds::bv_rMin], m_boundValues[DiscBounds::bv_rMax]);
41 }
42 
43 
44 bool
45 Trk::DiscBounds::operator==(const Trk::SurfaceBounds& sbo) const
46 {
47  // check the type first not to compare apples with oranges
48  const Trk::DiscBounds* discbo = dynamic_cast<const Trk::DiscBounds*>(&sbo);
49  if (!discbo)
50  return false;
51  return (m_boundValues == discbo->m_boundValues);
52 }
53 
54 bool
55 Trk::DiscBounds::inside(const Amg::Vector2D& locpo, double tol1, double tol2) const
56 {
57  double alpha = std::abs(locpo[locPhi] - m_boundValues[DiscBounds::bv_averagePhi]);
58  if (alpha > M_PI)
59  alpha = 2 * M_PI - alpha;
60  bool insidePhi =
61  (alpha <= (m_boundValues[DiscBounds::bv_halfPhiSector] + tol2));
62  return (locpo[locR] > (m_boundValues[DiscBounds::bv_rMin] - tol1) &&
63  locpo[locR] < (m_boundValues[DiscBounds::bv_rMax] + tol1) &&
64  insidePhi);
65 }
66 
67 bool
68 Trk::DiscBounds::inside(const Amg::Vector2D& locpo, const BoundaryCheck& bchk) const
69 {
70  if (bchk.bcType == 0 || bchk.nSigmas == 0 ||
71  m_boundValues[DiscBounds::bv_rMin] != 0 ||
72  m_boundValues[DiscBounds::bv_halfPhiSector] != M_PI)
73  return DiscBounds::inside(locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
74 
75  // a fast FALSE
76  sincosCache scResult = bchk.FastSinCos(locpo(1, 0));
77  double dx = bchk.nSigmas * sqrt(bchk.lCovariance(0, 0));
78  double dy =
79  bchk.nSigmas * sqrt(scResult.sinC * scResult.sinC * bchk.lCovariance(0, 0) +
80  locpo(0, 0) * locpo(0, 0) * scResult.cosC *
81  scResult.cosC * bchk.lCovariance(1, 1) +
82  2 * scResult.cosC * scResult.sinC * locpo(0, 0) *
83  bchk.lCovariance(0, 1));
84  double max_ell = dx > dy ? dx : dy;
85  if (locpo(0, 0) > (m_boundValues[DiscBounds::bv_rMax] + max_ell))
86  return false;
87  // a fast TRUE
88  double min_ell = dx < dy ? dx : dy;
89  if (locpo(0, 0) < (m_boundValues[DiscBounds::bv_rMax] + min_ell))
90  return true;
91 
92  // we are not using the KDOP approach here but rather a highly optimized one
93  class EllipseCollisionTest
94  {
95  private:
96  int m_maxIterations;
97  bool iterate(double x,
98  double y,
99  double c0x,
100  double c0y,
101  double c2x,
102  double c2y,
103  double rr) const
104  {
105  std::vector<double> innerPolygonCoef(m_maxIterations + 1);
106  std::vector<double> outerPolygonCoef(m_maxIterations + 1);
107  /*
108  t2______t4
109  --_ \
110  --_ \ /¨¨¨ ¨¨\
111  t1 = (0, 0) ( t )
112  | \ \__ _ /
113  | \
114  | t3
115  | /
116  | /
117  t0
118  */
119  for (int t = 1; t <= m_maxIterations; t++) {
120  int numNodes = 4 << t;
121  innerPolygonCoef[t] = 0.5 / cos(2.0*M_PI / numNodes);
122  double c1x = (c0x + c2x) * innerPolygonCoef[t];
123  double c1y = (c0y + c2y) * innerPolygonCoef[t];
124  double tx = x - c1x; // t indicates a translated coordinate
125  double ty = y - c1y;
126  if (tx * tx + ty * ty <= rr) {
127  return true; // collision with t1
128  }
129  double t2x = c2x - c1x;
130  double t2y = c2y - c1y;
131  if (tx * t2x + ty * t2y >= 0 &&
132  tx * t2x + ty * t2y <= t2x * t2x + t2y * t2y &&
133  (ty * t2x - tx * t2y >= 0 ||
134  rr * (t2x * t2x + t2y * t2y) >=
135  (ty * t2x - tx * t2y) * (ty * t2x - tx * t2y))) {
136  return true; // collision with t1---t2
137  }
138  double t0x = c0x - c1x;
139  double t0y = c0y - c1y;
140  if (tx * t0x + ty * t0y >= 0 &&
141  tx * t0x + ty * t0y <= t0x * t0x + t0y * t0y &&
142  (ty * t0x - tx * t0y <= 0 ||
143  rr * (t0x * t0x + t0y * t0y) >=
144  (ty * t0x - tx * t0y) * (ty * t0x - tx * t0y))) {
145  return true; // collision with t1---t0
146  }
147  outerPolygonCoef[t] =
148  0.5 / (std::cos(M_PI / numNodes) * std::cos(M_PI / numNodes));
149  double c3x = (c0x + c1x) * outerPolygonCoef[t];
150  double c3y = (c0y + c1y) * outerPolygonCoef[t];
151  if ((c3x - x) * (c3x - x) + (c3y - y) * (c3y - y) < rr) {
152  c2x = c1x;
153  c2y = c1y;
154  continue; // t3 is inside circle
155  }
156  double c4x = c1x - c3x + c1x;
157  double c4y = c1y - c3y + c1y;
158  if ((c4x - x) * (c4x - x) + (c4y - y) * (c4y - y) < rr) {
159  c0x = c1x;
160  c0y = c1y;
161  continue; // t4 is inside circle
162  }
163  double t3x = c3x - c1x;
164  double t3y = c3y - c1y;
165  if (ty * t3x - tx * t3y <= 0 ||
166  rr * (t3x * t3x + t3y * t3y) >
167  (ty * t3x - tx * t3y) * (ty * t3x - tx * t3y)) {
168  if (tx * t3x + ty * t3y > 0) {
169  if (std::abs(tx * t3x + ty * t3y) <= t3x * t3x + t3y * t3y ||
170  (x - c3x) * (c0x - c3x) + (y - c3y) * (c0y - c3y) >= 0) {
171  c2x = c1x;
172  c2y = c1y;
173  continue; // circle center is inside t0---t1---t3
174  }
175  } else if (-(tx * t3x + ty * t3y) <= t3x * t3x + t3y * t3y ||
176  (x - c4x) * (c2x - c4x) + (y - c4y) * (c2y - c4y) >= 0) {
177  c0x = c1x;
178  c0y = c1y;
179  continue; // circle center is inside t1---t2---t4
180  }
181  }
182  return false; // no collision possible
183  }
184  return false; // out of iterations so it is unsure if there was a
185  // collision. But have to return something.
186  }
187 
188  public:
189  // test for collision between an ellipse of horizontal radius w and vertical
190  // radius h at (x0, y0) and a circle of radius r at (x1, y1)
191  bool collide(double x0,
192  double y0,
193  double w,
194  double h,
195  double x1,
196  double y1,
197  double r) const
198  {
199  double x = std::abs(x1 - x0);
200  double y = std::abs(y1 - y0);
201  if (x * x + (h - y) * (h - y) <= r * r ||
202  (w - x) * (w - x) + y * y <= r * r ||
203  x * h + y * w <= w * h // collision with (0, h)
204  || ((x * h + y * w - w * h) * (x * h + y * w - w * h) <=
205  r * r * (w * w + h * h) &&
206  x * w - y * h >= -h * h &&
207  x * w - y * h <= w * w)) { // collision with (0, h)---(w, 0)
208  return true;
209  } else {
210  if ((x - w) * (x - w) + (y - h) * (y - h) <= r * r ||
211  (x <= w && y - r <= h) || (y <= h && x - r <= w)) {
212  return iterate(
213  x, y, w, 0, 0, h, r * r); // collision within triangle (0, h) (w, h)
214  // (0, 0) is possible
215  }
216  return false;
217  }
218  }
219  explicit EllipseCollisionTest(int maxIterations) : m_maxIterations(maxIterations)
220  {
221 
222  }
223  };
224 
225  EllipseCollisionTest test(4);
226  // convert to cartesian coordinates
227  AmgMatrix(2, 2) covRotMatrix;
228  // cppcheck-suppress constStatement
229  covRotMatrix << scResult.cosC, -locpo(0, 0) * scResult.sinC, scResult.sinC,
230  locpo(0, 0) * scResult.cosC;
231  AmgMatrix(2, 2) lCovarianceCar =
232  covRotMatrix * bchk.lCovariance * covRotMatrix.transpose();
233  Amg::Vector2D locpoCar(covRotMatrix(1, 1), -covRotMatrix(0, 1));
234 
235  // ellipse is always at (0,0), surface is moved to ellipse position and then
236  // rotated
237  double w = bchk.nSigmas * sqrt(lCovarianceCar(0, 0));
238  double h = bchk.nSigmas * sqrt(lCovarianceCar(1, 1));
239  double x0 = 0;
240  double y0 = 0;
241  float theta =
242  (lCovarianceCar(1, 0) != 0 &&
243  (lCovarianceCar(1, 1) - lCovarianceCar(0, 0)) != 0)
244  ? .5 * bchk.FastArcTan(2 * lCovarianceCar(1, 0) /
245  (lCovarianceCar(1, 1) - lCovarianceCar(0, 0)))
246  : 0.;
247  scResult = bchk.FastSinCos(theta);
248  AmgMatrix(2, 2) rotMatrix;
249  rotMatrix << scResult.cosC, scResult.sinC, -scResult.sinC, scResult.cosC;
250  Amg::Vector2D tmp = rotMatrix * (-locpoCar);
251  double x1 = tmp(0, 0);
252  double y1 = tmp(1, 0);
253  double r = m_boundValues[DiscBounds::bv_rMax];
254  // check if ellipse and circle overlap and return result
255  return test.collide(x0, y0, w, h, x1, y1, r);
256 }
257 
258 double
259 Trk::DiscBounds::minDistance(const Amg::Vector2D& pos) const
260 {
261  const double pi2 = 2. * M_PI;
262 
263  double r = pos[locR];
264  if (r == 0.)
265  return m_boundValues[DiscBounds::bv_rMin];
266  double sf = 0.;
267  double dF = 0.;
268 
269  if (m_boundValues[DiscBounds::bv_halfPhiSector] < M_PI) {
270 
271  dF = std::abs(pos[locPhi] - m_boundValues[DiscBounds::bv_averagePhi]);
272  if (dF > M_PI)
273  dF = pi2 - dF;
274  dF -= m_boundValues[DiscBounds::bv_halfPhiSector];
275  sf = r * sin(dF);
276  if (dF > 0.)
277  r *= cos(dF);
278 
279  } else {
280  sf = -1.e+10;
281  }
282 
283  if (sf <= 0.) {
284 
285  double sr0 = m_boundValues[DiscBounds::bv_rMin] - r;
286  if (sr0 > 0.)
287  return sr0;
288  double sr1 = r - m_boundValues[DiscBounds::bv_rMax];
289  if (sr1 > 0.)
290  return sr1;
291  if (sf < sr0)
292  sf = sr0;
293  if (sf < sr1)
294  sf = sr1;
295  return sf;
296  }
297 
298  double sr0 = m_boundValues[DiscBounds::bv_rMin] - r;
299  if (sr0 > 0.)
300  return sqrt(sr0 * sr0 + sf * sf);
301  double sr1 = r - m_boundValues[DiscBounds::bv_rMax];
302  if (sr1 > 0.)
303  return sqrt(sr1 * sr1 + sf * sf);
304  return sf;
305 }
306 
307 // ostream operator overload
308 
309 MsgStream&
310 Trk::DiscBounds::dump(MsgStream& sl) const
311 {
312  sl << std::setiosflags(std::ios::fixed);
313  sl << std::setprecision(7);
314  sl << "Trk::DiscBounds: (innerRadius, outerRadius, averagePhi, hPhiSector) = ";
315  sl << "(" << this->rMin() << ", " << this->rMax() << ", " << this->averagePhi() << ", " << this->halfPhiSector()
316  << ")";
317  sl << std::setprecision(-1);
318  return sl;
319 }
320 
321 std::ostream&
322 Trk::DiscBounds::dump(std::ostream& sl) const
323 {
324  sl << std::setiosflags(std::ios::fixed);
325  sl << std::setprecision(7);
326  sl << "Trk::DiscBounds: (innerRadius, outerRadius, hPhiSector) = ";
327  sl << "(" << this->rMin() << ", " << this->rMax() << ", " << this->averagePhi() << ", " << this->halfPhiSector()
328  << ")";
329  sl << std::setprecision(-1);
330  return sl;
331 }
Trk::y
@ y
Definition: ParamDefs.h:56
beamspotman.r
def r
Definition: beamspotman.py:676
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:216
Trk::AmgMatrix
AmgMatrix(3, 3) NeutralParticleParameterCalculator
Definition: NeutralParticleParameterCalculator.cxx:233
DiscBounds.h
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition: GeoPrimitives.h:48
Trk::SurfaceBounds
Definition: SurfaceBounds.h:47
Trk::SurfaceBounds::swap
void swap(double &b1, double &b2)
Swap method to be called from DiscBounds or TrapezoidalBounds.
Definition: SurfaceBounds.h:133
python.PhysicalConstants.pi2
float pi2
Definition: PhysicalConstants.py:52
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Trk::DiscBounds::bv_rMin
@ bv_rMin
Definition: DiscBounds.h:50
Trk::locR
@ locR
Definition: ParamDefs.h:44
TrigInDetValidation_Base.test
test
Definition: TrigInDetValidation_Base.py:147
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
Trk::DiscBounds::bv_halfPhiSector
@ bv_halfPhiSector
Definition: DiscBounds.h:53
Trk::DiscBounds::inside
virtual bool inside(const Amg::Vector2D &locpo, double tol1=0., double tol2=0.) const override final
This method cheks if the radius given in the LocalPosition is inside [rMin,rMax] if only tol1 is give...
Definition: DiscBounds.cxx:55
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
Trk::DiscBounds::dump
virtual MsgStream & dump(MsgStream &sl) const override
Output Method for MsgStream.
Definition: DiscBounds.cxx:310
Trk::BoundaryCheck::toleranceLoc1
double toleranceLoc1
absolute tolerance in local 1 coordinate
Definition: BoundaryCheck.h:68
Trk::DiscBounds::m_boundValues
std::vector< TDD_real_t > m_boundValues
Internal members of the bounds (float/double)
Definition: DiscBounds.h:128
Trk::DiscBounds::bv_rMax
@ bv_rMax
Definition: DiscBounds.h:51
Trk::theta
@ theta
Definition: ParamDefs.h:66
Trk::DiscBounds::operator==
virtual bool operator==(const SurfaceBounds &sbo) const override
Equality operator.
Definition: DiscBounds.cxx:45
Trk::BoundaryCheck::bcType
BoundaryCheckType bcType
Definition: BoundaryCheck.h:70
Trk::sincosCache
Definition: BoundaryCheck.h:45
EllipseCollisionTest
Definition: AnnulusBounds.cxx:21
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
Trk::DiscBounds::minDistance
virtual double minDistance(const Amg::Vector2D &pos) const override final
Minimal distance to boundary ( > 0 if outside and <=0 if inside)
Definition: DiscBounds.cxx:259
Trk::DiscBounds::bv_averagePhi
@ bv_averagePhi
Definition: DiscBounds.h:52
Trk::sincosCache::sinC
double sinC
Definition: BoundaryCheck.h:46
Trk::BoundaryCheck::nSigmas
int nSigmas
allowed sigmas for chi2 boundary check
Definition: BoundaryCheck.h:67
Trk::locPhi
@ locPhi
local polar
Definition: ParamDefs.h:45
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
Trk::BoundaryCheck::FastArcTan
double FastArcTan(double x) const
makeTRTBarrelCans.dy
tuple dy
Definition: makeTRTBarrelCans.py:21
h
mapkey::sf
@ sf
Definition: TElectronEfficiencyCorrectionTool.cxx:38
Trk::BoundaryCheck
Definition: BoundaryCheck.h:51
Trk::sincosCache::cosC
double cosC
Definition: BoundaryCheck.h:47
Trk::DiscBounds::DiscBounds
DiscBounds()
Default Constructor.
Definition: DiscBounds.cxx:17
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
Trk::BoundaryCheck::FastSinCos
sincosCache FastSinCos(double x) const
Trk::BoundaryCheck::toleranceLoc2
double toleranceLoc2
absolute tolerance in local 2 coordinate
Definition: BoundaryCheck.h:69
rr
const boost::regex rr(r_r)
PrepareReferenceFile.iterate
def iterate(ROOT.TDirectory thisdir, ROOT.TDirectory targetdir, str prefix, typing.Pattern regex, bool excludeTrees)
Definition: PrepareReferenceFile.py:10
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
Trk::x
@ x
Definition: ParamDefs.h:55
TileDCSDataPlotter.tx
tx
Definition: TileDCSDataPlotter.py:878
Trk::DiscBounds
Definition: DiscBounds.h:44
Generated on Sun Dec 22 2024 21:09:26 for ATLAS Offline Software by   doxygen 1.8.18