ATLAS Offline Software
DiscTrapezoidalBounds.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 // DiscTrapezoidalBounds.cxx, (c) ATLAS Detector Software
8 
9 // Trk
11 // Gaudi
12 #include "GaudiKernel/MsgStream.h"
13 // STD
14 #include <iomanip>
15 #include <iostream>
16 
18  : m_boundValues(DiscTrapezoidalBounds::bv_length, 0.)
19 {}
20 
22  double maxhalfx,
23  double maxR,
24  double minR,
25  double avephi,
26  double stereo)
27  : m_boundValues(DiscTrapezoidalBounds::bv_length, 0.)
28 {
31 
36 
41 
44 
45  double hmax =
48 
49  double hmin =
52 
55 }
56 
58  : Trk::SurfaceBounds()
59  , m_boundValues(disctrbo.m_boundValues)
60 {}
61 
62 
65 {
66  if (this != &disctrbo)
67  m_boundValues = disctrbo.m_boundValues;
68  return *this;
69 }
70 
71 bool
73 {
74  // check the type first not to compare apples with oranges
75  const Trk::DiscTrapezoidalBounds* disctrbo = dynamic_cast<const Trk::DiscTrapezoidalBounds*>(&sbo);
76  if (!disctrbo)
77  return false;
78  return (m_boundValues == disctrbo->m_boundValues);
79 }
80 
81 bool
82 Trk::DiscTrapezoidalBounds::inside(const Amg::Vector2D& locpo, double, double) const
83 {
84  double rMedium = m_boundValues[DiscTrapezoidalBounds::bv_rCenter];
85  double phi = m_boundValues[DiscTrapezoidalBounds::bv_averagePhi];
86 
87  Amg::Vector2D cartCenter(rMedium * cos(phi), rMedium * sin(phi));
88 
89  Amg::Vector2D cartPos(locpo[Trk::locR] * cos(locpo[Trk::locPhi]),
90  locpo[Trk::locR] * sin(locpo[Trk::locPhi]));
91 
92  Amg::Vector2D Pos = cartPos - cartCenter;
93 
94  Amg::Vector2D DeltaPos(Pos[Trk::locX] * sin(phi) - Pos[Trk::locY] * cos(phi),
95  Pos[Trk::locY] * sin(phi) + Pos[Trk::locX] * cos(phi));
96 
97  bool insideX = (std::abs(DeltaPos[Trk::locX]) <
98  m_boundValues[DiscTrapezoidalBounds::bv_maxHalfX]);
99  bool insideY = (std::abs(DeltaPos[Trk::locY]) <
100  m_boundValues[DiscTrapezoidalBounds::bv_halfY]);
101 
102  if (!insideX || !insideY)
103  return false;
104 
105  if (std::abs(DeltaPos[Trk::locX]) <
106  m_boundValues[DiscTrapezoidalBounds::bv_minHalfX])
107  return true;
108 
109  double m = 2. * m_boundValues[DiscTrapezoidalBounds::bv_halfY] /
110  (m_boundValues[DiscTrapezoidalBounds::bv_maxHalfX] -
111  m_boundValues[DiscTrapezoidalBounds::bv_minHalfX]);
112 
113  double q = m_boundValues[DiscTrapezoidalBounds::bv_halfY] *
114  (m_boundValues[DiscTrapezoidalBounds::bv_maxHalfX] +
115  m_boundValues[DiscTrapezoidalBounds::bv_minHalfX]) /
116  (m_boundValues[DiscTrapezoidalBounds::bv_minHalfX] -
117  m_boundValues[DiscTrapezoidalBounds::bv_maxHalfX]);
118 
119  bool inside = (DeltaPos[Trk::locY] > (m * std::abs(DeltaPos[Trk::locX]) + q));
120 
121  return inside;
122 }
123 
124 bool
126  const BoundaryCheck& bchk) const
127 {
128  if (bchk.bcType == 0 || bchk.nSigmas == 0 ||
129  m_boundValues[DiscTrapezoidalBounds::bv_rMin] != 0)
131  locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
132  double alpha =
133  std::abs(locpo[locPhi] - m_boundValues[DiscTrapezoidalBounds::bv_averagePhi]);
134  if (alpha > M_PI)
135  alpha = 2 * M_PI - alpha;
136  // a fast FALSE
137  sincosCache scResult = bchk.FastSinCos(locpo(1, 0));
138  double dx = bchk.nSigmas * sqrt(bchk.lCovariance(0, 0));
139  double dy =
140  bchk.nSigmas * sqrt(scResult.sinC * scResult.sinC * bchk.lCovariance(0, 0) +
141  locpo(0, 0) * locpo(0, 0) * scResult.cosC *
142  scResult.cosC * bchk.lCovariance(1, 1) +
143  2 * scResult.cosC * scResult.sinC * locpo(0, 0) *
144  bchk.lCovariance(0, 1));
145  double max_ell = dx > dy ? dx : dy;
146  if (locpo(0, 0) >
147  (m_boundValues[DiscTrapezoidalBounds::bv_rMax] *
149  cos(alpha) +
150  max_ell))
151  return false;
152  // a fast TRUE
153  double min_ell = dx < dy ? dx : dy;
154  if (locpo(0, 0) <
155  (m_boundValues[DiscTrapezoidalBounds::bv_rMax] *
157  cos(alpha) +
158  min_ell))
159  return true;
160 
161  // we are not using the KDOP approach here but rather a highly optimized one
163  {
164  private:
165  int m_maxIterations;
166  bool iterate(double x,
167  double y,
168  double c0x,
169  double c0y,
170  double c2x,
171  double c2y,
172  double rr) const
173  {
174  std::vector<double> innerPolygonCoef(m_maxIterations + 1);
175  std::vector<double> outerPolygonCoef(m_maxIterations + 1);
176  /*
177  t2______t4
178  --_ \
179  --_ \ /¨¨¨ ¨¨\
180  t1 = (0, 0) ( t )
181  | \ \__ _ /
182  | \
183  | t3
184  | /
185  | /
186  t0
187  */
188  for (int t = 1; t <= m_maxIterations; t++) {
189  int numNodes = 4 << t;
190  innerPolygonCoef[t] = 0.5 / cos(2.0f*M_PI / numNodes);
191  double c1x = (c0x + c2x) * innerPolygonCoef[t];
192  double c1y = (c0y + c2y) * innerPolygonCoef[t];
193  double tx = x - c1x; // t indicates a translated coordinate
194  double ty = y - c1y;
195  if (tx * tx + ty * ty <= rr) {
196  return true; // collision with t1
197  }
198  double t2x = c2x - c1x;
199  double t2y = c2y - c1y;
200  if (tx * t2x + ty * t2y >= 0 &&
201  tx * t2x + ty * t2y <= t2x * t2x + t2y * t2y &&
202  (ty * t2x - tx * t2y >= 0 ||
203  rr * (t2x * t2x + t2y * t2y) >=
204  (ty * t2x - tx * t2y) * (ty * t2x - tx * t2y))) {
205  return true; // collision with t1---t2
206  }
207  double t0x = c0x - c1x;
208  double t0y = c0y - c1y;
209  if (tx * t0x + ty * t0y >= 0 &&
210  tx * t0x + ty * t0y <= t0x * t0x + t0y * t0y &&
211  (ty * t0x - tx * t0y <= 0 ||
212  rr * (t0x * t0x + t0y * t0y) >=
213  (ty * t0x - tx * t0y) * (ty * t0x - tx * t0y))) {
214  return true; // collision with t1---t0
215  }
216  outerPolygonCoef[t] =
217  0.5 / (std::cos(M_PI / numNodes) * std::cos(M_PI / numNodes));
218  double c3x = (c0x + c1x) * outerPolygonCoef[t];
219  double c3y = (c0y + c1y) * outerPolygonCoef[t];
220  if ((c3x - x) * (c3x - x) + (c3y - y) * (c3y - y) < rr) {
221  c2x = c1x;
222  c2y = c1y;
223  continue; // t3 is inside circle
224  }
225  double c4x = c1x - c3x + c1x;
226  double c4y = c1y - c3y + c1y;
227  if ((c4x - x) * (c4x - x) + (c4y - y) * (c4y - y) < rr) {
228  c0x = c1x;
229  c0y = c1y;
230  continue; // t4 is inside circle
231  }
232  double t3x = c3x - c1x;
233  double t3y = c3y - c1y;
234  if (ty * t3x - tx * t3y <= 0 ||
235  rr * (t3x * t3x + t3y * t3y) >
236  (ty * t3x - tx * t3y) * (ty * t3x - tx * t3y)) {
237  if (tx * t3x + ty * t3y > 0) {
238  if (std::abs(tx * t3x + ty * t3y) <= t3x * t3x + t3y * t3y ||
239  (x - c3x) * (c0x - c3x) + (y - c3y) * (c0y - c3y) >= 0) {
240  c2x = c1x;
241  c2y = c1y;
242  continue; // circle center is inside t0---t1---t3
243  }
244  } else if (-(tx * t3x + ty * t3y) <= t3x * t3x + t3y * t3y ||
245  (x - c4x) * (c2x - c4x) + (y - c4y) * (c2y - c4y) >= 0) {
246  c0x = c1x;
247  c0y = c1y;
248  continue; // circle center is inside t1---t2---t4
249  }
250  }
251  return false; // no collision possible
252  }
253  return false; // out of iterations so it is unsure if there was a
254  // collision. But have to return something.
255  }
256 
257  public:
258  // test for collision between an ellipse of horizontal radius w and vertical
259  // radius h at (x0, y0) and a circle of radius r at (x1, y1)
260  bool collide(double x0,
261  double y0,
262  double w,
263  double h,
264  double x1,
265  double y1,
266  double r) const
267  {
268  double x = std::abs(x1 - x0);
269  double y = std::abs(y1 - y0);
270  if (x * x + (h - y) * (h - y) <= r * r ||
271  (w - x) * (w - x) + y * y <= r * r ||
272  x * h + y * w <= w * h // collision with (0, h)
273  || ((x * h + y * w - w * h) * (x * h + y * w - w * h) <=
274  r * r * (w * w + h * h) &&
275  x * w - y * h >= -h * h &&
276  x * w - y * h <= w * w)) { // collision with (0, h)---(w, 0)
277  return true;
278  } else {
279  if ((x - w) * (x - w) + (y - h) * (y - h) <= r * r ||
280  (x <= w && y - r <= h) || (y <= h && x - r <= w)) {
281  return iterate(
282  x, y, w, 0, 0, h, r * r); // collision within triangle (0, h) (w, h)
283  // (0, 0) is possible
284  }
285  return false;
286  }
287  }
288  explicit EllipseCollisionTest(int maxIterations) : m_maxIterations(maxIterations)
289  {
290 
291  }
292  };
293 
295  // convert to cartesian coordinates
296  AmgMatrix(2, 2) covRotMatrix;
297  // cppcheck-suppress constStatement
298  covRotMatrix << scResult.cosC, -locpo(0, 0) * scResult.sinC, scResult.sinC,
299  locpo(0, 0) * scResult.cosC;
300  AmgMatrix(2, 2) lCovarianceCar =
301  covRotMatrix * bchk.lCovariance * covRotMatrix.transpose();
302  Amg::Vector2D locpoCar(covRotMatrix(1, 1), -covRotMatrix(0, 1));
303 
304  // ellipse is always at (0,0), surface is moved to ellipse position and then
305  // rotated
306  double w = bchk.nSigmas * sqrt(lCovarianceCar(0, 0));
307  double h = bchk.nSigmas * sqrt(lCovarianceCar(1, 1));
308  double x0 = 0;
309  double y0 = 0;
310  float theta =
311  (lCovarianceCar(1, 0) != 0 &&
312  (lCovarianceCar(1, 1) - lCovarianceCar(0, 0)) != 0)
313  ? .5 * bchk.FastArcTan(2 * lCovarianceCar(1, 0) /
314  (lCovarianceCar(1, 1) - lCovarianceCar(0, 0)))
315  : 0.;
316  scResult = bchk.FastSinCos(theta);
317  AmgMatrix(2, 2) rotMatrix;
318  rotMatrix << scResult.cosC, scResult.sinC, -scResult.sinC, scResult.cosC;
319  Amg::Vector2D tmp = rotMatrix * (-locpoCar);
320  double x1 = tmp(0, 0);
321  double y1 = tmp(1, 0);
322  double r = m_boundValues[DiscTrapezoidalBounds::bv_rMax];
323  // check if ellipse and circle overlap and return result
324  return test.collide(x0, y0, w, h, x1, y1, r);
325 }
326 
327 
328 double
330 {
331  const double pi2 = 2. * M_PI;
332  double alpha = std::abs(pos[locPhi]);
333  if (alpha > M_PI)
334  alpha = pi2 - alpha;
335 
336  double r = pos[locR];
337  if (r == 0.)
338  return m_boundValues[DiscTrapezoidalBounds::bv_rMin] * cos(m_boundValues[DiscTrapezoidalBounds::bv_halfPhiSector]) /
339  cos(alpha);
340 
341  // check if it is external in R
342  double sr0 = m_boundValues[DiscTrapezoidalBounds::bv_rMin] *
343  cos(m_boundValues[DiscTrapezoidalBounds::bv_halfPhiSector]) / cos(alpha) -
344  r;
345  if (sr0 > 0.)
346  return sr0;
347  double sr1 = r - m_boundValues[DiscTrapezoidalBounds::bv_rMax] *
348  cos(m_boundValues[DiscTrapezoidalBounds::bv_halfPhiSector]) / cos(alpha);
349  if (sr1 > 0.)
350  return sr1;
351 
352  // check if it is external in phi
353  if ((alpha - m_boundValues[DiscTrapezoidalBounds::bv_halfPhiSector]) > 0.)
354  return r * std::abs(sin(alpha - m_boundValues[DiscTrapezoidalBounds::bv_halfPhiSector]));
355 
356  // if here it is inside:
357  // Evaluate the distance from the 4 segments
358  double dist[4] = { sr0,
359  sr1,
360  -r * sin(m_boundValues[DiscTrapezoidalBounds::bv_halfPhiSector] - alpha),
361  -r * sin(m_boundValues[DiscTrapezoidalBounds::bv_halfPhiSector] + alpha) };
362 
363  return *std::max_element(dist, dist + 4);
364 }
365 
366 // ostream operator overload
367 
368 MsgStream&
370 {
371  sl << std::setiosflags(std::ios::fixed);
372  sl << std::setprecision(7);
373  sl << "Trk::DiscTrapezoidalBounds: (innerRadius, outerRadius, hMinX, hMaxX, hlengthY, hPhiSector, averagePhi, "
374  "rCenter, stereo) = ";
375  sl << "(" << this->rMin() << ", " << this->rMax() << ", " << this->minHalflengthX() << ", " << this->maxHalflengthX()
376  << ", " << this->halflengthY() << ", " << this->halfPhiSector() << ", " << this->averagePhi() << ", "
377  << this->rCenter() << ", " << this->stereo() << ")";
378  sl << std::setprecision(-1);
379  return sl;
380 }
381 
382 std::ostream&
383 Trk::DiscTrapezoidalBounds::dump(std::ostream& sl) const
384 {
385  sl << std::setiosflags(std::ios::fixed);
386  sl << std::setprecision(7);
387  sl << "Trk::DiscTrapezoidalBounds: (innerRadius, outerRadius, hMinX, hMaxX, hlengthY, hPhiSector, averagePhi, "
388  "rCenter, stereo) = ";
389  sl << "(" << this->rMin() << ", " << this->rMax() << ", " << this->minHalflengthX() << ", " << this->maxHalflengthX()
390  << ", " << this->halflengthY() << ", " << this->halfPhiSector() << ", " << this->averagePhi() << ", "
391  << this->rCenter() << ", " << this->stereo() << ")";
392  sl << std::setprecision(-1);
393  return sl;
394 }
Trk::DiscTrapezoidalBounds::bv_rMin
@ bv_rMin
Definition: DiscTrapezoidalBounds.h:48
Trk::y
@ y
Definition: ParamDefs.h:62
beamspotman.r
def r
Definition: beamspotman.py:676
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:216
python.CaloRecoConfig.f
f
Definition: CaloRecoConfig.py:127
Trk::AmgMatrix
AmgMatrix(3, 3) NeutralParticleParameterCalculator
Definition: NeutralParticleParameterCalculator.cxx:233
Trk::DiscTrapezoidalBounds::operator=
DiscTrapezoidalBounds & operator=(const DiscTrapezoidalBounds &disctrbo)
Assignment operator.
Definition: DiscTrapezoidalBounds.cxx:64
Trk::DiscTrapezoidalBounds::bv_stereo
@ bv_stereo
Definition: DiscTrapezoidalBounds.h:56
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
Trk::locX
@ locX
Definition: ParamDefs.h:43
Trk::locY
@ locY
local cartesian
Definition: ParamDefs.h:44
Trk::DiscTrapezoidalBounds::bv_rCenter
@ bv_rCenter
Definition: DiscTrapezoidalBounds.h:55
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::DiscTrapezoidalBounds::DiscTrapezoidalBounds
DiscTrapezoidalBounds()
Default Constructor.
Definition: DiscTrapezoidalBounds.cxx:17
Trk::DiscTrapezoidalBounds::bv_maxHalfX
@ bv_maxHalfX
Definition: DiscTrapezoidalBounds.h:51
Trk::locR
@ locR
Definition: ParamDefs.h:50
TrigInDetValidation_Base.test
test
Definition: TrigInDetValidation_Base.py:144
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
DiscTrapezoidalBounds.h
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
Trk::DiscTrapezoidalBounds
Definition: DiscTrapezoidalBounds.h:42
Trk::BoundaryCheck::toleranceLoc1
double toleranceLoc1
absolute tolerance in local 1 coordinate
Definition: BoundaryCheck.h:68
Trk::DiscTrapezoidalBounds::minDistance
virtual double minDistance(const Amg::Vector2D &pos) const override final
Minimal distance to boundary ( > 0 if outside and <=0 if inside)
Definition: DiscTrapezoidalBounds.cxx:329
Trk::theta
@ theta
Definition: ParamDefs.h:72
Trk::DiscTrapezoidalBounds::bv_rMax
@ bv_rMax
Definition: DiscTrapezoidalBounds.h:49
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::DiscTrapezoidalBounds::operator==
virtual bool operator==(const SurfaceBounds &sbo) const override
Equality operator.
Definition: DiscTrapezoidalBounds.cxx:72
Trk::DiscTrapezoidalBounds::m_boundValues
std::vector< TDD_real_t > m_boundValues
Internal members of the bounds (float/double)
Definition: DiscTrapezoidalBounds.h:144
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
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Trk::DiscTrapezoidalBounds::bv_minHalfX
@ bv_minHalfX
Definition: DiscTrapezoidalBounds.h:50
Trk::DiscTrapezoidalBounds::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: DiscTrapezoidalBounds.cxx:82
Trk::DiscTrapezoidalBounds::bv_halfY
@ bv_halfY
Definition: DiscTrapezoidalBounds.h:52
Trk::locPhi
@ locPhi
local polar
Definition: ParamDefs.h:51
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
Trk::BoundaryCheck::FastArcTan
double FastArcTan(double x) const
Trk::inside
@ inside
Definition: PropDirection.h:29
makeTRTBarrelCans.dy
tuple dy
Definition: makeTRTBarrelCans.py:21
Trk::DiscTrapezoidalBounds::dump
virtual MsgStream & dump(MsgStream &sl) const override
Output Method for MsgStream.
Definition: DiscTrapezoidalBounds.cxx:369
h
Trk::BoundaryCheck
Definition: BoundaryCheck.h:51
Trk::sincosCache::cosC
double cosC
Definition: BoundaryCheck.h:47
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
Trk::BoundaryCheck::FastSinCos
sincosCache FastSinCos(double x) const
hmax
double hmax(TH1 *&h)
Definition: listroot.cxx:115
extractSporadic.q
list q
Definition: extractSporadic.py:98
Trk::DiscTrapezoidalBounds::bv_averagePhi
@ bv_averagePhi
Definition: DiscTrapezoidalBounds.h:54
Trk::BoundaryCheck::toleranceLoc2
double toleranceLoc2
absolute tolerance in local 2 coordinate
Definition: BoundaryCheck.h:69
Trk::phi
@ phi
Definition: ParamDefs.h:81
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
Trk::DiscTrapezoidalBounds::stereo
double stereo() const
This method returns the stereo angle.
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
Trk::x
@ x
Definition: ParamDefs.h:61
TileDCSDataPlotter.tx
tx
Definition: TileDCSDataPlotter.py:878
Trk::DiscTrapezoidalBounds::bv_halfPhiSector
@ bv_halfPhiSector
Definition: DiscTrapezoidalBounds.h:53