ATLAS Offline Software
DiamondBounds.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // DiamondBounds.cxx, (c) ATLAS Detector Software
8 
9 // Trk
10 #include "TrkSurfaces/DiamondBounds.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::DiamondBounds::DiamondBounds()
20  : m_boundValues(DiamondBounds::bv_length, 0.)
21  , m_alpha1(0.)
22  , m_alpha2(0.)
23 {}
24 
25 // constructor from arguments I
26 Trk::DiamondBounds::DiamondBounds(double minhalex, double medhalex,
27  double maxhalex, double haley1, double haley2)
28  : m_boundValues(DiamondBounds::bv_length, 0.), m_alpha1(0.), m_alpha2(0.) {
29  m_boundValues[DiamondBounds::bv_minHalfX] = minhalex;
30  m_boundValues[DiamondBounds::bv_medHalfX] = medhalex;
31  m_boundValues[DiamondBounds::bv_maxHalfX] = maxhalex;
32  m_boundValues[DiamondBounds::bv_halfY1] = haley1;
33  m_boundValues[DiamondBounds::bv_halfY2] = haley2;
34  if (minhalex > maxhalex)
35  swap(m_boundValues[DiamondBounds::bv_minHalfX],
36  m_boundValues[DiamondBounds::bv_maxHalfX]);
37  DiamondBounds::initCache();
38 }
39 
40 void Trk::DiamondBounds::initCache() {
41  m_alpha1 = atan2(m_boundValues[DiamondBounds::bv_medHalfX] -
42  m_boundValues[DiamondBounds::bv_minHalfX],
43  2. * m_boundValues[DiamondBounds::bv_halfY1]);
44  m_alpha2 = atan2(m_boundValues[DiamondBounds::bv_medHalfX] -
45  m_boundValues[DiamondBounds::bv_maxHalfX],
46  2. * m_boundValues[DiamondBounds::bv_halfY2]);
47 }
48 
49 bool
50 Trk::DiamondBounds::operator==(const Trk::SurfaceBounds& sbo) const
51 {
52  // check the type first not to compare apples with oranges
53  const Trk::DiamondBounds* diabo = dynamic_cast<const Trk::DiamondBounds*>(&sbo);
54  if (!diabo)
55  return false;
56  return (m_boundValues == diabo->m_boundValues);
57 }
58 
59 // checking if inside bounds (Full symmetrical Diamond)
60 bool
61 Trk::DiamondBounds::inside(const Amg::Vector2D& locpo, double tol1, double tol2) const
62 {
63  return this->insideFull(locpo, tol1, tol2);
64 }
65 
66 bool
67 Trk::DiamondBounds::inside(const Amg::Vector2D& locpo,
68  const BoundaryCheck& bchk) const
69 {
70  if (bchk.bcType == 0)
71  return DiamondBounds::inside(locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
72 
73  // a fast FALSE
74  double max_ell = bchk.lCovariance(0, 0) > bchk.lCovariance(1, 1)
75  ? bchk.lCovariance(0, 0)
76  : bchk.lCovariance(1, 1);
77  double limit = bchk.nSigmas * sqrt(max_ell);
78  if (locpo[Trk::locY] < -2 * m_boundValues[DiamondBounds::bv_halfY1] - limit)
79  return false;
80  if (locpo[Trk::locY] > 2 * m_boundValues[DiamondBounds::bv_halfY2] + limit)
81  return false;
82  // a fast FALSE
83  double fabsX = std::abs(locpo[Trk::locX]);
84  if (fabsX > (m_boundValues[DiamondBounds::bv_medHalfX] + limit))
85  return false;
86  // a fast TRUE
87  double min_ell = bchk.lCovariance(0, 0) < bchk.lCovariance(1, 1)
88  ? bchk.lCovariance(0, 0)
89  : bchk.lCovariance(1, 1);
90  limit = bchk.nSigmas * sqrt(min_ell);
91  if (fabsX < (fmin(m_boundValues[DiamondBounds::bv_minHalfX],
92  m_boundValues[DiamondBounds::bv_maxHalfX]) -
93  limit))
94  return true;
95  // a fast TRUE
96  if (std::abs(locpo[Trk::locY]) < (fmin(m_boundValues[DiamondBounds::bv_halfY1],
97  m_boundValues[DiamondBounds::bv_halfY2]) -
98  limit))
99  return true;
100 
101  // compute KDOP and axes for surface polygon
102  std::vector<KDOP> elementKDOP(5);
103  std::vector<Amg::Vector2D> elementP(6);
104  float theta =
105  (bchk.lCovariance(1, 0) != 0 &&
106  (bchk.lCovariance(1, 1) - bchk.lCovariance(0, 0)) != 0)
107  ? .5 * bchk.FastArcTan(2 * bchk.lCovariance(1, 0) /
108  (bchk.lCovariance(1, 1) - bchk.lCovariance(0, 0)))
109  : 0.;
110  sincosCache scResult = bchk.FastSinCos(theta);
111  AmgMatrix(2, 2) rotMatrix;
112  rotMatrix << scResult.cosC, scResult.sinC, -scResult.sinC, scResult.cosC;
113  AmgMatrix(2, 2) normal;
114  // cppcheck-suppress constStatement
115  normal << 0, -1, 1, 0;
116  // ellipse is always at (0,0), surface is moved to ellipse position and then
117  // rotated
118  Amg::Vector2D p =
119  Amg::Vector2D(-m_boundValues[DiamondBounds::bv_minHalfX],
120  -2. * m_boundValues[DiamondBounds::bv_halfY1]);
121  elementP[0] = (rotMatrix * (p - locpo));
122  p = Amg::Vector2D (-m_boundValues[DiamondBounds::bv_medHalfX], 0.);
123  elementP[1] = (rotMatrix * (p - locpo));
124  p = Amg::Vector2D (-m_boundValues[DiamondBounds::bv_maxHalfX],
125  2. * m_boundValues[DiamondBounds::bv_halfY2]);
126  elementP[2] = (rotMatrix * (p - locpo));
127  p = Amg::Vector2D (m_boundValues[DiamondBounds::bv_maxHalfX],
128  2. * m_boundValues[DiamondBounds::bv_halfY2]);
129  elementP[3] = (rotMatrix * (p - locpo));
130  p = Amg::Vector2D (m_boundValues[DiamondBounds::bv_medHalfX], 0.);
131  elementP[4] = (rotMatrix * (p - locpo));
132  p = Amg::Vector2D (m_boundValues[DiamondBounds::bv_minHalfX],
133  -2. * m_boundValues[DiamondBounds::bv_halfY1]);
134  elementP[5] = (rotMatrix * (p - locpo));
135  std::vector<Amg::Vector2D> axis = { normal * (elementP[1] - elementP[0]),
136  normal * (elementP[2] - elementP[1]),
137  normal * (elementP[3] - elementP[2]),
138  normal * (elementP[4] - elementP[3]),
139  normal * (elementP[5] - elementP[4]) };
140  bchk.ComputeKDOP(elementP, axis, elementKDOP);
141  // compute KDOP for error ellipse
142  std::vector<KDOP> errelipseKDOP(5);
143  bchk.ComputeKDOP(bchk.EllipseToPoly(3), axis, errelipseKDOP);
144  // check if KDOPs overlap and return result
145  return bchk.TestKDOPKDOP(elementKDOP, errelipseKDOP);
146 }
147 
148 // checking if inside bounds (Full symmetrical Diamond)
149 bool
150 Trk::DiamondBounds::insideFull(const Amg::Vector2D& locpo, double tol1, double tol2) const
151 {
152  // validity check
153  if (!m_boundValues[DiamondBounds::bv_halfY1] && !m_boundValues[DiamondBounds::bv_minHalfX]) return false;
154 
155  // quick False (radial direction)
156  if (locpo[Trk::locY] < -2. * m_boundValues[DiamondBounds::bv_halfY1] - tol2) return false;
157  if (locpo[Trk::locY] > 2. * m_boundValues[DiamondBounds::bv_halfY2] + tol2) return false;
158 
159  double absX = std::abs(locpo[Trk::locX]);
160 
161  // quick False (transverse directon)
162  if (absX > m_boundValues[DiamondBounds::bv_medHalfX] + tol1) return false;
163 
164  // quick True
165  if (absX < std::min(m_boundValues[DiamondBounds::bv_minHalfX], m_boundValues[DiamondBounds::bv_maxHalfX]) + tol1) return true;
166 
168  const double& halfBaseUp = locpo[Trk::locY] < 0 ? m_boundValues[DiamondBounds::bv_medHalfX] : m_boundValues[DiamondBounds::bv_maxHalfX];
169  const double& halfBaseLo = locpo[Trk::locY] < 0 ? m_boundValues[DiamondBounds::bv_minHalfX] : m_boundValues[DiamondBounds::bv_medHalfX];
170  const double& halfH = locpo[Trk::locY] < 0 ? m_boundValues[DiamondBounds::bv_halfY1] : m_boundValues[DiamondBounds::bv_halfY2];
171  double k = halfH ? 0.5*(halfBaseUp - halfBaseLo)/halfH : 0.;
172  double sign = (k < 0) ? -1. : 1.;
173  return (absX - tol1 <= m_boundValues[DiamondBounds::bv_medHalfX] + k * (locpo[Trk::locY] + sign*tol2));
174 
175 }
176 
177 // opening angle in point A
178 double
179 Trk::DiamondBounds::alpha1() const
180 {
181  return m_alpha1;
182 }
183 
184 // opening angle in point A'
185 double
186 Trk::DiamondBounds::alpha2() const
187 {
188  return m_alpha2;
189 }
190 
191 double
192 Trk::DiamondBounds::minDistance(const Amg::Vector2D& pos) const
193 {
194  const int Np = 6;
195 
196  double y1 = 2. * m_boundValues[DiamondBounds::bv_halfY1];
197  double y2 = 2. * m_boundValues[DiamondBounds::bv_halfY2];
198 
199  double X[6] = { -m_boundValues[DiamondBounds::bv_minHalfX], -m_boundValues[DiamondBounds::bv_medHalfX],
200  -m_boundValues[DiamondBounds::bv_maxHalfX], m_boundValues[DiamondBounds::bv_maxHalfX],
201  m_boundValues[DiamondBounds::bv_medHalfX], m_boundValues[DiamondBounds::bv_minHalfX] };
202  double Y[6] = { -y1, 0., y2, y2, 0., -y1 };
203 
204  double dm = 1.e+20;
205  double Ao = 0.;
206  bool in = true;
207 
208  for (int i = 0; i != Np; ++i) {
209 
210  int j = (i == Np-1 ? 0 : i+1);
211 
212  double x = X[i] - pos[0];
213  double y = Y[i] - pos[1];
214  double dx = X[j] - X[i];
215  double dy = Y[j] - Y[i];
216  double A = x * dy - y * dx;
217  double S = -(x * dx + y * dy);
218 
219  if (S <= 0.) {
220  double d = x * x + y * y;
221  if (d < dm)
222  dm = d;
223  } else {
224  double a = dx * dx + dy * dy;
225  if (S <= a) {
226  double d = (A * A) / a;
227  if (d < dm)
228  dm = d;
229  }
230  }
231  if (i && in && Ao * A < 0.)
232  in = false;
233  Ao = A;
234  }
235  if (in)
236  return -sqrt(dm);
237  return sqrt(dm);
238 }
239 
240 // ostream operator overload
241 
242 MsgStream&
243 Trk::DiamondBounds::dump(MsgStream& sl) const
244 {
245  sl << std::setiosflags(std::ios::fixed);
246  sl << std::setprecision(7);
247  sl << "Trk::DiamondBounds: (minHlenghtX, medHlengthX, maxHlengthX, hlengthY1, hlengthY2 ) = ";
248  sl << "(" << m_boundValues[DiamondBounds::bv_minHalfX] << ", " << m_boundValues[DiamondBounds::bv_medHalfX] << ", "
249  << m_boundValues[DiamondBounds::bv_maxHalfX] << ", " << m_boundValues[DiamondBounds::bv_halfY1] << ", "
250  << m_boundValues[DiamondBounds::bv_halfY2] << ")";
251  sl << std::setprecision(-1);
252  return sl;
253 }
254 
255 std::ostream&
256 Trk::DiamondBounds::dump(std::ostream& sl) const
257 {
258  sl << std::setiosflags(std::ios::fixed);
259  sl << std::setprecision(7);
260  sl << "Trk::DiamondBounds: (minHlenghtX, medHlengthX, maxHlengthX, hlengthY1, hlengthY2 ) = ";
261  sl << "(" << m_boundValues[DiamondBounds::bv_minHalfX] << ", " << m_boundValues[DiamondBounds::bv_medHalfX] << ", "
262  << m_boundValues[DiamondBounds::bv_maxHalfX] << ", " << m_boundValues[DiamondBounds::bv_halfY1] << ", "
263  << m_boundValues[DiamondBounds::bv_halfY2] << ")";
264  sl << std::setprecision(-1);
265  return sl;
266 }
Trk::y
@ y
Definition: ParamDefs.h:56
Trk::AmgMatrix
AmgMatrix(3, 3) NeutralParticleParameterCalculator
Definition: NeutralParticleParameterCalculator.cxx:233
Trk::DiamondBounds::insideFull
bool insideFull(const Amg::Vector2D &locpo, double tol1=0., double tol2=0.) const
inside() method for a full symmetric diamond
Definition: DiamondBounds.cxx:150
Trk::locX
@ locX
Definition: ParamDefs.h:37
Trk::locY
@ locY
local cartesian
Definition: ParamDefs.h:38
Trk::DiamondBounds::minDistance
virtual double minDistance(const Amg::Vector2D &pos) const override final
Minimal distance to boundary ( > 0 if outside and <=0 if inside)
Definition: DiamondBounds.cxx:192
Trk::DiamondBounds::alpha2
double alpha2() const
This method returns the opening angle alpha in point A'
Definition: DiamondBounds.cxx:186
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
Trk::DiamondBounds::operator==
virtual bool operator==(const SurfaceBounds &diabo) const override
Equality operator.
Definition: DiamondBounds.cxx:50
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
Trk::DiamondBounds::initCache
virtual void initCache() override
initialize the alpha1/2 cache - needed also for object persistency
Definition: DiamondBounds.cxx:40
Trk::DiamondBounds::inside
virtual bool inside(const Amg::Vector2D &locpo, double tol1=0., double tol2=0.) const override final
The orientation of the Diamond is according to the figure.
Definition: DiamondBounds.cxx:61
JetTiledMap::S
@ S
Definition: TiledEtaPhiMap.h:44
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
Monitored::X
@ X
Definition: HistogramFillerUtils.h:24
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
DiamondBounds.h
Trk::DiamondBounds::m_boundValues
std::vector< TDD_real_t > m_boundValues
Internal parameters stored in the geometry.
Definition: DiamondBounds.h:140
Trk::DiamondBounds::bv_medHalfX
@ bv_medHalfX
Definition: DiamondBounds.h:44
A
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
Trk::DiamondBounds
Definition: DiamondBounds.h:37
lumiFormat.i
int i
Definition: lumiFormat.py:85
Trk::theta
@ theta
Definition: ParamDefs.h:66
makeTRTBarrelCans.y2
tuple y2
Definition: makeTRTBarrelCans.py:18
Trk::BoundaryCheck::bcType
BoundaryCheckType bcType
Definition: BoundaryCheck.h:70
Trk::DiamondBounds::bv_halfY1
@ bv_halfY1
Definition: DiamondBounds.h:46
sign
int sign(int a)
Definition: TRT_StrawNeighbourSvc.h:107
Trk::sincosCache
Definition: BoundaryCheck.h:45
Trk::BoundaryCheck::TestKDOPKDOP
bool TestKDOPKDOP(const std::vector< KDOP > &a, const std::vector< KDOP > &b) const
Trk::DiamondBounds::DiamondBounds
DiamondBounds()
Default Constructor, needed for persistency.
Definition: DiamondBounds.cxx:19
ReadCellNoiseFromCool.dm
dm
Definition: ReadCellNoiseFromCool.py:235
min
#define min(a, b)
Definition: cfImp.cxx:40
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::DiamondBounds::bv_maxHalfX
@ bv_maxHalfX
Definition: DiamondBounds.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
Trk::DiamondBounds::dump
virtual MsgStream & dump(MsgStream &sl) const override final
Output Method for MsgStream.
Definition: DiamondBounds.cxx:243
a
TList * a
Definition: liststreamerinfos.cxx:10
Trk::BoundaryCheck
Definition: BoundaryCheck.h:51
Trk::sincosCache::cosC
double cosC
Definition: BoundaryCheck.h:47
Trk::DiamondBounds::bv_minHalfX
@ bv_minHalfX
Definition: DiamondBounds.h:43
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
Trk::x
@ x
Definition: ParamDefs.h:55
updateCoolNtuple.limit
int limit
Definition: updateCoolNtuple.py:45
Trk::DiamondBounds::alpha1
double alpha1() const
This method returns the opening angle alpha in point A
Definition: DiamondBounds.cxx:179
Trk::DiamondBounds::bv_halfY2
@ bv_halfY2
Definition: DiamondBounds.h:47
fitman.k
k
Definition: fitman.py:528
Generated on Thu Nov 7 2024 21:13:41 for ATLAS Offline Software by   doxygen 1.8.18