ATLAS Offline Software
TrapezoidBounds.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 // TrapezoidBounds.cxx, (c) ATLAS Detector Software
8 
9 // Trk
10 #include "TrkSurfaces/TrapezoidBounds.h"
11 // Gaudi
12 #include "GaudiKernel/MsgStream.h"
13 // STD
14 #include <cmath>
15 #include <iomanip>
16 #include <iostream>
17 
18 // default constructor
19 Trk::TrapezoidBounds::TrapezoidBounds()
20  : m_boundValues(TrapezoidBounds::bv_length, 0.)
21  , m_alpha(0.)
22  , m_beta(0.)
23 {}
24 
25 // constructor from arguments I
26 Trk::TrapezoidBounds::TrapezoidBounds(double minhalex, double maxhalex, double haley)
27  : m_boundValues(TrapezoidBounds::bv_length, 0.)
28  , m_alpha(0.)
29  , m_beta(0.)
30 {
31  m_boundValues[TrapezoidBounds::bv_minHalfX] = std::abs(minhalex);
32  m_boundValues[TrapezoidBounds::bv_maxHalfX] = std::abs(maxhalex);
33  m_boundValues[TrapezoidBounds::bv_halfY] = std::abs(haley);
34  if (m_boundValues[TrapezoidBounds::bv_minHalfX] > m_boundValues[TrapezoidBounds::bv_maxHalfX])
35  swap(m_boundValues[TrapezoidBounds::bv_minHalfX], m_boundValues[TrapezoidBounds::bv_maxHalfX]);
36 }
37 
38 // constructor from arguments II
39 Trk::TrapezoidBounds::TrapezoidBounds(double minhalex, double haley, double alpha, double beta)
40  : m_boundValues(TrapezoidBounds::bv_length, 0.)
41  , m_alpha(alpha)
42  , m_beta(beta)
43 {
44  double gamma = (alpha > beta) ? (alpha - 0.5 * M_PI) : (beta - 0.5 * M_PI);
45  // now fill them
46  m_boundValues[TrapezoidBounds::bv_minHalfX] = std::abs(minhalex);
47  m_boundValues[TrapezoidBounds::bv_maxHalfX] = minhalex + (2. * m_boundValues[TrapezoidBounds::bv_halfY]) * tan(gamma);
48  m_boundValues[TrapezoidBounds::bv_halfY] = std::abs(haley);
49 }
50 
51 
52 bool
53 Trk::TrapezoidBounds::operator==(const Trk::SurfaceBounds& sbo) const
54 {
55  // check the type first not to compare apples with oranges
56  const Trk::TrapezoidBounds* trabo = dynamic_cast<const Trk::TrapezoidBounds*>(&sbo);
57  if (!trabo)
58  return false;
59  return (m_boundValues == trabo->m_boundValues);
60 }
61 
62 bool
63 Trk::TrapezoidBounds::inside(const Amg::Vector2D& locpo,
64  const BoundaryCheck& bchk) const
65 {
66  if (bchk.bcType == 0)
67  return TrapezoidBounds::inside(
68  locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
69 
70  // a fast FALSE
71  double fabsY = std::abs(locpo[Trk::locY]);
72  double max_ell = bchk.lCovariance(0, 0) > bchk.lCovariance(1, 1)
73  ? bchk.lCovariance(0, 0)
74  : bchk.lCovariance(1, 1);
75  double limit = bchk.nSigmas * sqrt(max_ell);
76  if (fabsY > (m_boundValues[TrapezoidBounds::bv_halfY] + limit))
77  return false;
78  // a fast FALSE
79  double fabsX = std::abs(locpo[Trk::locX]);
80  if (fabsX > (m_boundValues[TrapezoidBounds::bv_maxHalfX] + limit))
81  return false;
82  // a fast TRUE
83  double min_ell = bchk.lCovariance(0, 0) < bchk.lCovariance(1, 1)
84  ? bchk.lCovariance(0, 0)
85  : bchk.lCovariance(1, 1);
86  limit = bchk.nSigmas * sqrt(min_ell);
87  if (fabsX < (m_boundValues[TrapezoidBounds::bv_minHalfX] + limit) &&
88  fabsY < (m_boundValues[TrapezoidBounds::bv_halfY] + limit))
89  return true;
90 
91  // compute KDOP and axes for surface polygon
92  std::vector<KDOP> elementKDOP(3);
93  std::vector<Amg::Vector2D> elementP(4);
94  float theta =
95  (bchk.lCovariance(1, 0) != 0 &&
96  (bchk.lCovariance(1, 1) - bchk.lCovariance(0, 0)) != 0)
97  ? .5 * bchk.FastArcTan(2 * bchk.lCovariance(1, 0) /
98  (bchk.lCovariance(1, 1) - bchk.lCovariance(0, 0)))
99  : 0.;
100  sincosCache scResult = bchk.FastSinCos(theta);
101  AmgMatrix(2, 2) rotMatrix;
102  rotMatrix << scResult.cosC, scResult.sinC, -scResult.sinC, scResult.cosC;
103  AmgMatrix(2, 2) normal;
104  // cppcheck-suppress constStatement
105  normal << 0, -1, 1, 0;
106  // ellipse is always at (0,0), surface is moved to ellipse position and then
107  // rotated
108  Amg::Vector2D p = Amg::Vector2D(m_boundValues[TrapezoidBounds::bv_minHalfX],
109  -m_boundValues[TrapezoidBounds::bv_halfY]);
110  elementP[0] = (rotMatrix * (p - locpo));
111  p = Amg::Vector2D (-m_boundValues[TrapezoidBounds::bv_minHalfX],
112  -m_boundValues[TrapezoidBounds::bv_halfY]);
113  elementP[1] = (rotMatrix * (p - locpo));
114  scResult = bchk.FastSinCos(m_beta);
115  p = Amg::Vector2D (m_boundValues[TrapezoidBounds::bv_minHalfX] +
116  (2. * m_boundValues[TrapezoidBounds::bv_halfY]) *
117  (scResult.sinC / scResult.cosC),
118  m_boundValues[TrapezoidBounds::bv_halfY]);
119  elementP[2] = (rotMatrix * (p - locpo));
120  scResult = bchk.FastSinCos(m_alpha);
121  p = Amg::Vector2D (-(m_boundValues[TrapezoidBounds::bv_minHalfX] +
122  (2. * m_boundValues[TrapezoidBounds::bv_halfY]) *
123  (scResult.sinC / scResult.cosC)),
124  m_boundValues[TrapezoidBounds::bv_halfY]);
125  elementP[3] = (rotMatrix * (p - locpo));
126  std::vector<Amg::Vector2D> axis = { normal * (elementP[1] - elementP[0]),
127  normal * (elementP[3] - elementP[1]),
128  normal * (elementP[2] - elementP[0]) };
129  bchk.ComputeKDOP(elementP, axis, elementKDOP);
130  // compute KDOP for error ellipse
131  std::vector<KDOP> errelipseKDOP(3);
132  bchk.ComputeKDOP(bchk.EllipseToPoly(3), axis, errelipseKDOP);
133  // check if KDOPs overlap and return result
134  return bchk.TestKDOPKDOP(elementKDOP, errelipseKDOP);
135 }
136 
137 // checking if inside bounds
138 bool
139 Trk::TrapezoidBounds::inside(const Amg::Vector2D& locpo, double tol1, double tol2) const
140 {
141  if (m_alpha == 0.)
142  return insideFull(locpo, tol1, tol2);
143  return (insideFull(locpo, tol1, tol2) && !insideExclude(locpo, tol1, tol2));
144 }
145 
146 // checking if inside bounds (Full symmetrical Trapezoid)
147 bool
148 Trk::TrapezoidBounds::insideFull(const Amg::Vector2D& locpo, double tol1, double tol2) const
149 {
150  // the cases:
151  double fabsX = std::abs(locpo[Trk::locX]);
152  double fabsY = std::abs(locpo[Trk::locY]);
153  // (1) a fast FALSE
154  if (fabsY > (m_boundValues[TrapezoidBounds::bv_halfY] + tol2))
155  return false;
156  // (2) a fast FALSE
157  if (fabsX > (m_boundValues[TrapezoidBounds::bv_maxHalfX] + tol1))
158  return false;
159  // (3) a fast TRUE
160  if (fabsX < (m_boundValues[TrapezoidBounds::bv_minHalfX] + tol1))
161  return true;
162  // (4) particular case - a rectangle
163  if (m_boundValues[TrapezoidBounds::bv_maxHalfX] == m_boundValues[TrapezoidBounds::bv_minHalfX])
164  return true;
165  // (5) /** use basic calculation of a straight line */
166  double k = 2.0 * m_boundValues[TrapezoidBounds::bv_halfY] /
167  (m_boundValues[TrapezoidBounds::bv_maxHalfX] - m_boundValues[TrapezoidBounds::bv_minHalfX]) *
168  ((locpo[Trk::locX] > 0.) ? 1.0 : -1.0);
169  double d =
170  -std::abs(k) * 0.5 * (m_boundValues[TrapezoidBounds::bv_maxHalfX] + m_boundValues[TrapezoidBounds::bv_minHalfX]);
171  return (isAbove(locpo, tol1, tol2, k, d));
172 }
173 
174 // checking if local point is inside the exclude area
175 bool
176 Trk::TrapezoidBounds::insideExclude(const Amg::Vector2D& locpo, double tol1, double tol2) const
177 {
178 
179  // line a
180  bool alphaBiggerBeta(m_alpha > m_beta);
181  double ka = alphaBiggerBeta ? tan(M_PI - m_alpha) : tan(m_alpha);
182  double kb = alphaBiggerBeta ? tan(M_PI - m_beta) : tan(m_beta);
183  double sign = alphaBiggerBeta ? -1. : 1.;
184  double da = -m_boundValues[TrapezoidBounds::bv_halfY] + sign * ka * m_boundValues[TrapezoidBounds::bv_minHalfX];
185  double db = -m_boundValues[TrapezoidBounds::bv_halfY] + sign * kb * m_boundValues[TrapezoidBounds::bv_minHalfX];
186 
187  return (isAbove(locpo, tol1, tol2, ka, da) && isAbove(locpo, tol1, tol2, kb, db));
188 }
189 
190 // checking if local point lies above a line
191 bool
192 Trk::TrapezoidBounds::isAbove(const Amg::Vector2D& locpo, double tol1, double tol2, double k, double d)
193 {
194  // the most tolerant approach for tol1 and tol2
195  double sign = k > 0. ? -1. : +1.;
196  return (locpo[Trk::locY] + tol2 > (k * (locpo[Trk::locX] + sign * tol1) + d));
197 }
198 
199 double
200 Trk::TrapezoidBounds::minDistance(const Amg::Vector2D& pos) const
201 {
202  const int Np = 4;
203 
204  double xl = -m_boundValues[TrapezoidBounds::bv_maxHalfX];
205  double xr = m_boundValues[TrapezoidBounds::bv_maxHalfX];
206  if (m_alpha != 0.) {
207  xl = -m_boundValues[TrapezoidBounds::bv_minHalfX] - 2. * tan(m_alpha) * m_boundValues[TrapezoidBounds::bv_halfY];
208  } else if (m_beta != 0.) {
209  xr = m_boundValues[TrapezoidBounds::bv_minHalfX] + 2. * tan(m_beta) * m_boundValues[TrapezoidBounds::bv_halfY];
210  }
211  double X[4] = { -m_boundValues[TrapezoidBounds::bv_minHalfX], xl, xr, m_boundValues[TrapezoidBounds::bv_minHalfX] };
212  double Y[4] = { -m_boundValues[TrapezoidBounds::bv_halfY],
213  m_boundValues[TrapezoidBounds::bv_halfY],
214  m_boundValues[TrapezoidBounds::bv_halfY],
215  -m_boundValues[TrapezoidBounds::bv_halfY] };
216 
217  double dm = 1.e+20;
218  double Ao = 0.;
219  bool in = true;
220 
221  for (int i = 0; i != Np; ++i) {
222 
223  int j = (i == Np-1 ? 0 : i+1);
224 
225  double x = X[i] - pos[0];
226  double y = Y[i] - pos[1];
227  double dx = X[j] - X[i];
228  double dy = Y[j] - Y[i];
229  double A = x * dy - y * dx;
230  double S = -(x * dx + y * dy);
231 
232  if (S <= 0.) {
233  double d = x * x + y * y;
234  if (d < dm)
235  dm = d;
236  } else {
237  double a = dx * dx + dy * dy;
238  if (S <= a) {
239  double d = (A * A) / a;
240  if (d < dm)
241  dm = d;
242  }
243  }
244  if (i && in && Ao * A < 0.)
245  in = false;
246  Ao = A;
247  }
248  if (in){
249  return -sqrt(dm);
250  }
251  return sqrt(dm);
252 }
253 
254 // ostream operator overload
255 
256 MsgStream&
257 Trk::TrapezoidBounds::dump(MsgStream& sl) const
258 {
259  sl << std::setiosflags(std::ios::fixed);
260  sl << std::setprecision(7);
261  sl << "Trk::TrapezoidBounds: (minHlenghtX, maxHlengthX, hlengthY) = "
262  << "(" << m_boundValues[TrapezoidBounds::bv_minHalfX] << ", " << m_boundValues[TrapezoidBounds::bv_maxHalfX]
263  << ", " << m_boundValues[TrapezoidBounds::bv_halfY] << ")";
264  sl << std::setprecision(-1);
265  return sl;
266 }
267 
268 std::ostream&
269 Trk::TrapezoidBounds::dump(std::ostream& sl) const
270 {
271  sl << std::setiosflags(std::ios::fixed);
272  sl << std::setprecision(7);
273  sl << "Trk::TrapezoidBounds: (minHlenghtX, maxHlengthX, hlengthY) = "
274  << "(" << m_boundValues[TrapezoidBounds::bv_minHalfX] << ", " << m_boundValues[TrapezoidBounds::bv_maxHalfX]
275  << ", " << m_boundValues[TrapezoidBounds::bv_halfY] << ")";
276  sl << std::setprecision(-1);
277  return sl;
278 }
TrapezoidBounds.h
Trk::TrapezoidBounds::insideFull
bool insideFull(const Amg::Vector2D &locpo, double tol1=0., double tol2=0.) const
inside() method for a full symmetric trapezoid
Definition: TrapezoidBounds.cxx:148
Trk::y
@ y
Definition: ParamDefs.h:56
Trk::TrapezoidBounds::alpha
double alpha() const
This method returns the opening angle alpha in point A (negative local phi)
Trk::TrapezoidBounds::m_boundValues
std::vector< TDD_real_t > m_boundValues
Definition: TrapezoidBounds.h:177
Trk::AmgMatrix
AmgMatrix(3, 3) NeutralParticleParameterCalculator
Definition: NeutralParticleParameterCalculator.cxx:233
Trk::TrapezoidBounds::bv_minHalfX
@ bv_minHalfX
Definition: TrapezoidBounds.h:49
Trk::locX
@ locX
Definition: ParamDefs.h:37
Trk::locY
@ locY
local cartesian
Definition: ParamDefs.h:38
Trk::TrapezoidBounds::minDistance
virtual double minDistance(const Amg::Vector2D &pos) const override
Minimal distance to boundary ( > 0 if outside and <=0 if inside)
Definition: TrapezoidBounds.cxx:200
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition: GeoPrimitives.h:48
hist_file_dump.d
d
Definition: hist_file_dump.py:137
Trk::SurfaceBounds
Definition: SurfaceBounds.h:47
CaloCondBlobAlgs_fillNoiseFromASCII.db
db
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:43
Trk::SurfaceBounds::swap
void swap(double &b1, double &b2)
Swap method to be called from DiscBounds or TrapezoidalBounds.
Definition: SurfaceBounds.h:133
yodamerge_tmp.axis
list axis
Definition: yodamerge_tmp.py:241
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Trk::TrapezoidBounds::inside
virtual bool inside(const Amg::Vector2D &locpo, double tol1=0., double tol2=0.) const override
The orientation of the Trapezoid is according to the figure above, in words: the shorter of the two p...
Definition: TrapezoidBounds.cxx:139
Trk::TrapezoidBounds::bv_halfY
@ bv_halfY
Definition: TrapezoidBounds.h:51
JetTiledMap::S
@ S
Definition: TiledEtaPhiMap.h:44
Monitored::X
@ X
Definition: HistogramFillerUtils.h:24
dqt_zlumi_alleff_HIST.A
A
Definition: dqt_zlumi_alleff_HIST.py:110
Trk::BoundaryCheck::EllipseToPoly
std::vector< Amg::Vector2D > EllipseToPoly(int resolution=3) const
Trk::BoundaryCheck::toleranceLoc1
double toleranceLoc1
absolute tolerance in local 1 coordinate
Definition: BoundaryCheck.h:68
A
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
TrigVtx::gamma
@ gamma
Definition: TrigParticleTable.h:26
lumiFormat.i
int i
Definition: lumiFormat.py:85
Trk::theta
@ theta
Definition: ParamDefs.h:66
Trk::TrapezoidBounds::isAbove
static bool isAbove(const Amg::Vector2D &locpo, double tol1, double tol2, double k, double d)
isAbove() method for checking whether a point lies above or under a straight line
Definition: TrapezoidBounds.cxx:192
Trk::BoundaryCheck::bcType
BoundaryCheckType bcType
Definition: BoundaryCheck.h:70
sign
int sign(int a)
Definition: TRT_StrawNeighbourSvc.h:107
Trk::TrapezoidBounds::bv_maxHalfX
@ bv_maxHalfX
Definition: TrapezoidBounds.h:50
Trk::sincosCache
Definition: BoundaryCheck.h:45
Trk::BoundaryCheck::TestKDOPKDOP
bool TestKDOPKDOP(const std::vector< KDOP > &a, const std::vector< KDOP > &b) const
Trk::TrapezoidBounds::operator==
virtual bool operator==(const SurfaceBounds &trabo) const override
Equality operator.
Definition: TrapezoidBounds.cxx:53
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
checkxAOD.kb
kb
Definition: Tools/PyUtils/bin/checkxAOD.py:104
ReadCellNoiseFromCool.dm
dm
Definition: ReadCellNoiseFromCool.py:235
Monitored::Y
@ Y
Definition: HistogramFillerUtils.h:24
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::BoundaryCheck::ComputeKDOP
void ComputeKDOP(const std::vector< Amg::Vector2D > &v, const std::vector< Amg::Vector2D > &KDOPAxes, std::vector< KDOP > &kdop) const
Each Bounds has a method inside, which checks if a LocalPosition is inside the bounds.
Trk::TrapezoidBounds::TrapezoidBounds
TrapezoidBounds()
Default Constructor, needed for persistency.
Definition: TrapezoidBounds.cxx:19
Trk::TrapezoidBounds
Definition: TrapezoidBounds.h:43
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
Trk::BoundaryCheck::FastArcTan
double FastArcTan(double x) const
makeTRTBarrelCans.dy
tuple dy
Definition: makeTRTBarrelCans.py:21
Trk::TrapezoidBounds::dump
virtual MsgStream & dump(MsgStream &sl) const override
Output Method for MsgStream.
Definition: TrapezoidBounds.cxx:257
a
TList * a
Definition: liststreamerinfos.cxx:10
Trk::BoundaryCheck
Definition: BoundaryCheck.h:51
Trk::sincosCache::cosC
double cosC
Definition: BoundaryCheck.h:47
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
Trk::TrapezoidBounds::beta
double beta() const
This method returns the opening angle beta in point B (positive local phi)
Trk::BoundaryCheck::FastSinCos
sincosCache FastSinCos(double x) const
Trk::TrapezoidBounds::insideExclude
bool insideExclude(const Amg::Vector2D &locpo, double tol1=0., double tol2=0.) const
inside() method for the triangular exclude area for an arbitrary trapezoid
Definition: TrapezoidBounds.cxx:176
Trk::BoundaryCheck::toleranceLoc2
double toleranceLoc2
absolute tolerance in local 2 coordinate
Definition: BoundaryCheck.h:69
Trk::x
@ x
Definition: ParamDefs.h:55
MuonParameters::beta
@ beta
Definition: MuonParamDefs.h:144
updateCoolNtuple.limit
int limit
Definition: updateCoolNtuple.py:45
fitman.k
k
Definition: fitman.py:528
Generated on Sun Dec 22 2024 21:20:40 for ATLAS Offline Software by   doxygen 1.8.18