Belle II Software  release-05-02-19
ThreeHitFilters.h
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Jakob Lettenbichler *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #pragma once
12 
13 #include <TVector3.h>
14 #include "tracking/vxdCaTracking/TwoHitFilters.h"
15 #include "tracking/vxdCaTracking/FilterExceptions.h"
16 #include <framework/logging/Logger.h>
17 
18 namespace Belle2 {
25  class ThreeHitFilters {
26  public:
27 
31  m_radiusCalculated(false),
32  m_radius(0.),
33  m_x2(0.),
34  m_y2(0.),
35  m_z2(0.)
36  {
37  m_hitA.SetXYZ(0., 0., 0.);
38  m_hitB.SetXYZ(0., 0., 0.);
39  m_hitC.SetXYZ(0., 0., 0.);
40  m_vecAB.SetXYZ(0., 0., 0.);
41  m_vecBC.SetXYZ(0., 0., 0.);
42  resetMagneticField(1.5);
43  }
44 
46  ThreeHitFilters(TVector3& outerHit, TVector3& centerHit, TVector3& innerHit, double magneticFieldStrength = 1.5):
48  m_radiusCalculated(false),
49  m_radius(0.),
50  m_hitA(outerHit),
51  m_hitB(centerHit),
52  m_hitC(innerHit),
53  m_vecAB(outerHit - centerHit),
54  m_vecBC(centerHit - innerHit)
55  {
56  m_x2 = m_vecAB[0] * m_vecBC[0];
57  m_y2 = m_vecAB[1] * m_vecBC[1];
58  m_z2 = m_vecAB[2] * m_vecBC[2];
59  resetMagneticField(magneticFieldStrength);
60  }
61 
62 
64  ~ThreeHitFilters() {}
65 
67  void resetValues(TVector3& outerHit, TVector3& centerHit, TVector3& innerHit)
68  {
69  m_radiusCalculated = false;
71  m_radius = 0.;
72  m_hitA = outerHit;
73  m_hitB = centerHit;
74  m_hitC = innerHit;
75  m_vecAB = outerHit - centerHit;
76  m_vecBC = centerHit - innerHit;
77 
78  m_x2 = m_vecAB[0] * m_vecBC[0];
79  m_y2 = m_vecAB[1] * m_vecBC[1];
80  m_z2 = m_vecAB[2] * m_vecBC[2];
81  }
82 
83 
84 
90  void resetMagneticField(double magneticFieldStrength = 1.5) { m_magneticFieldFactor = magneticFieldStrength * 0.00299710; }
91 
92 
93 
95  double getMagneticField() { return m_magneticFieldFactor / 0.00299710; }
96 
97 
98 
100  double filterNan(double value) { return m_twoHitFilter.filterNan(value); }
101 
102 
103 
105  double calcAngle3D()
106  {
107  double angle = ((m_x2 + m_y2 + m_z2) / (m_vecAB.Mag2() *
108  m_vecBC.Mag2())); // fullCalc would be acos(m_vecAB.Dot(m_vecBC) / m_vecAB.Mag()*m_vecBC.Mag())
109  return filterNan(angle);
110  } // return unit: none (calculation for degrees is incomplete, if you want readable numbers, use fullAngle3D instead)
111 
112 
113 
115  double fullAngle3D()
116  {
117  double angle = acos(m_vecAB.Dot(m_vecBC) / (m_vecAB.Mag() * m_vecBC.Mag())); // 0-pi
118  angle = (angle * (180. / M_PI));
119  return filterNan(angle);
120  } // return unit: ° (0 - 180°)
121 
122 
123 
125  double calcAngleXY()
126  {
127  double angle = ((m_x2 + m_y2) / (m_vecAB.Perp2() * m_vecBC.Perp2())); // fullAngle:
128  return filterNan(angle);
129  } // return unit: none (calculation for degrees is incomplete, if you want readable numbers, use fullAngleXY instead)
130 
131 
132 
134  double fullAngleXY()
135  {
136  double angle = fullAngle2D(m_vecAB, m_vecBC); // 0-pi
137  angle = (angle * (180. / M_PI));
138  return filterNan(angle);
139  } // return unit: ° (0 - 180°)
140 
141 
142 
144  double calcAngleRZ()
145  {
146  TVector3 rzVecAB(m_vecAB.Perp(), m_vecAB[2], 0.);
147  TVector3 rzVecBC(m_vecBC.Perp(), m_vecBC[2], 0.);
148  return calcAngle2D(rzVecAB, rzVecBC);
149  } // return unit: none (calculation for degrees is incomplete, if you want readable numbers, use fullAngleRZ instead)
150 
151 
152 
154  double fullAngleRZ()
155  {
156  TVector3 rzVecAB(m_vecAB.Perp(), m_vecAB[2], 0.);
157  TVector3 rzVecBC(m_vecBC.Perp(), m_vecBC[2], 0.);
158  double angle = fullAngle2D(rzVecAB, rzVecBC); // 0-pi
159  angle = (angle * (180. / M_PI));
160  return filterNan(angle);
161  } // return unit: ° (0 - 180°)
162 
163 
164 
166  double calcCircleDist2IP()
167  {
168  checkCalcRadius();
169 
170  return (fabs(m_centerABC.Perp() - m_radius)); // distance of closest approach of circle to the IP
171  } // return unit: cm
172 
173 
174 
176  double calcPt()
177  {
178  checkCalcRadius();
179 
180  return calcPt(m_radius);
181  } // return unit: GeV/c
182 
183 
184 
186  double calcPt(double radius)
187  {
188  sanityCheckRadius(radius);
189 
190  return (m_magneticFieldFactor * radius);
191  } // return unit: GeV/c
192 
193 
194 
196  double calcDeltaSlopeRZ()
197  {
199  double slopeAB = m_twoHitFilter.calcSlopeRZ();
201  double slopeBC = m_twoHitFilter.calcSlopeRZ();
202 
203  return filterNan(slopeBC - slopeAB);
204  } // return unit: none
205 
206 
207 
212  double calcDeltaSOverZ()
213  {
215 
216  TVector3 points2hitA = m_hitA - m_centerABC;
217  TVector3 points2hitB = m_hitB - m_centerABC;
218  TVector3 points2hitC = m_hitC - m_centerABC;
219  double alfaAB = fullAngle2D(points2hitA, points2hitB);
220  double alfaBC = fullAngle2D(points2hitB, points2hitC);
221  //return filterNan( (alfaAB * m_vecBC[2]) - (alfaBC *m_vecAB[2]) );
222  return (alfaAB * m_vecBC[2]) - (alfaBC *
223  m_vecAB[2]); // equals to alfaAB/dZAB and alfaBC/dZBC, but this solution here can not produce a division by zero
224  } // return unit: radians*cm
225 
226 
227 
229  double calcDeltaSlopeZOverS()
230  {
231  checkCalcRadius();
232 
233  TVector3 points2hitA = m_hitA - m_centerABC;
234  TVector3 points2hitB = m_hitB - m_centerABC;
235  TVector3 points2hitC = m_hitC - m_centerABC;
236  double alfaABr = fullAngle2D(points2hitA, points2hitB) * m_radius;
237  double alfaBCr = fullAngle2D(points2hitB, points2hitC) * m_radius;
238 
239  return filterNan((asin(m_vecAB[2] / alfaABr)) - asin(m_vecBC[2] / alfaBCr)); // Beware of z>r!
240  }
241 
242 
243 
245  double calcHelixParameterFit()
246  {
248  TVector3 points2hitA = m_hitA - m_centerABC;
249  TVector3 points2hitB = m_hitB - m_centerABC;
250  TVector3 points2hitC = m_hitC - m_centerABC;
251  double alfaAB = calcAngle2D(points2hitA, points2hitB);
252  double alfaBC = calcAngle2D(points2hitB, points2hitC);
253  // real calculation: ratio is (m_vecij[2] = deltaZ): alfaAB/deltaZab : alfaBC/deltaZbc, the following equation saves two times '/'
254  return filterNan((alfaAB * m_vecBC[2]) / (alfaBC * m_vecAB[2]));
255  } // return unit: none
256 
257 
258 
260  double calcHelixFit() { return calcHelixParameterFit(); }
261 
262 
263 
265  double calcAngle2D(TVector3& vecA, TVector3& vecB)
266  {
267  double angle = ((vecA[0] * vecB[0] + vecA[1] * vecB[1]) / sqrt(vecA.Perp2() * vecB.Perp2()));
268  return filterNan(angle);
269  }
270 
274  double fullAngle2D(TVector3& vecA, TVector3& vecB)
275  {
276  return acos(calcAngle2D(vecA, vecB));
277  //return filterNan(angle);
278  }
279 
280 
281 
283  double calcRadius(TVector3& a, TVector3& b, TVector3& c, TVector3& circleCenter)
284  {
285  return ((circleCenter - a).Perp() + (circleCenter - b).Perp() + (circleCenter - c).Perp()) /
286  3.; // = radius in [cm], sign here not needed. normally: signKappaAB/normAB1
287  } // used by calcPt() and calcCircleDist2IP()
288 
289 
290 
292  void calcCircleCenter(TVector3& a, TVector3& b, TVector3& c, TVector3& circleCenter)
293  {
294  // calculates the intersection point using Cramer's rule.
295  // x_1+s*n_1==x_2+t*n_2 --> n_1 *s - n_2 *t == x_2 - x_1 --> http://en.wikipedia.org/wiki/Cramer%27s_rule
296  double inX = b[0] - c[0]; // x value of the normal vector of the inner segment (b-c)
297  double inY = b[1] - c[1]; // y value of the normal vector of the inner segment (b-c)
298  double outX = a[0] - b[0]; // x value of the normal vector of the outer segment (a-b)
299  double outY = a[1] - b[1]; // y value of the normal vector of the outer segment (a-b)
300 
301  //searching solution for Ax = b, aij are the matrix elements of A, bi are elements of b
302  double a11 = inY;
303  double a12 = -inX;
304  double a21 = -outY;
305  double a22 = outX;
306  double b1 = b[0] + outX * 0.5 - (c[0] + inX * 0.5);
307  double b2 = b[1] + outY * 0.5 - (c[1] + inY * 0.5);
308 
309  if (a11 * a22 == a12 * a21) { throw FilterExceptions::Straight_Line(); }
310 
311  double s = (b1 * a22 - b2 * a21) / (a11 * a22 - a12 * a21); //the determinant is zero iff the three hits are on a line in (x,y).
312 
313  circleCenter.SetXYZ(c[0] + inX * 0.5 + s * inY, c[1] + inY * 0.5 - s * inX, 0.);
314  }
315 
316 
317 
319  int calcSign(TVector3& a, TVector3& b, TVector3& c);
320 
321 
322 
324  int calcSign(TVector3& a, TVector3& b, TVector3& c, TVector3& sigma_a, TVector3& sigma_b, TVector3& sigma_c)
325  {
326  TVector3 c2b = b - c; c2b.SetZ(0.);
327  TVector3 b2a = a - b; b2a.SetZ(0.);
328  double angle = atan2(b2a[0], b2a[1]) - atan2(c2b[0], c2b[1]);
329  double sigmaan = (sigma_a.Mag() + sigma_b.Mag() + sigma_c.Mag()) / (3.*(c2b.Mag() +
330  b2a.Mag())); //TODO 1/3...mean of the sigmas. Possible improvement: Use a parameter instead, and determine with simulated events.
331  if (angle < (-sigmaan)) { return -1; }
332  else if (angle > sigmaan) {return 1; }
333  else { return 0; }
334  }
335 
336 
337 
338  protected:
339 
340 
342  void checkCalcCircleCenter()
343  {
344  if (m_circleCenterCalculated == false) {
347  }
348  }
349 
350 
352  void checkCalcRadius()
353  {
355  if (m_radiusCalculated == false) {
358  m_radiusCalculated = true;
359  }
360  }
361 
362 
364  void sanityCheckRadius(double radius)
365  {
366  if (fabs(radius) <
367  0.0000001) { // WARNING hardcoded value, is there a quasi-global value for such cases (this case, minimal accepted radius)
368  m_radiusCalculated = false;
369  m_circleCenterCalculated = false;
370  m_radius = 0.;
371  B2ERROR("sanityCheckRadius: given radius is " << radius << ", which is below threshold of " << 0.0000001 << ", throwing exception");
372  throw FilterExceptions::Circle_too_small();
373  }
374  }
375 
376 
377 
378  TwoHitFilters m_twoHitFilter;
381  bool m_radiusCalculated;
382  double m_radius;
383  double m_x2;
384  double m_y2;
385  double m_z2;
387  TVector3 m_centerABC;
388  TVector3 m_hitA;
389  TVector3 m_hitB;
390  TVector3 m_hitC;
391  TVector3 m_vecAB;
392  TVector3 m_vecBC;
394  }; //end class ThreeHitFilters
396 } //end namespace Belle2
Belle2::ThreeHitFilters::resetValues
void resetValues(TVector3 &outerHit, TVector3 &centerHit, TVector3 &innerHit)
Overrides Constructor-Setup.
Definition: ThreeHitFilters.h:75
Belle2::ThreeHitFilters::calcAngleRZ
double calcAngleRZ()
calculates the angle between the hits/vectors (RZ), returning unit: none (calculation for degrees is ...
Definition: ThreeHitFilters.h:152
Belle2::ThreeHitFilters::m_magneticFieldFactor
double m_magneticFieldFactor
is factor containing speed of light (c), the magnetic field (b) and the scaling factor s for conversi...
Definition: ThreeHitFilters.h:394
Belle2::ThreeHitFilters::m_hitB
TVector3 m_hitB
center hit (position relevant for useful filter calculation) used for the filter calculation
Definition: ThreeHitFilters.h:397
Belle2::ThreeHitFilters::calcRadius
double calcRadius(TVector3 &a, TVector3 &b, TVector3 &c, TVector3 &circleCenter)
calculates an estimation of the radius of given hits and existing estimation of circleCenter,...
Definition: ThreeHitFilters.h:291
Belle2::ThreeHitFilters::m_z2
double m_z2
internal intermediate value storing z^2, no enduser-relevance
Definition: ThreeHitFilters.h:393
Belle2::ThreeHitFilters::m_twoHitFilter
TwoHitFilters m_twoHitFilter
instance of TwoHitFilters-class used for some internal calculations
Definition: ThreeHitFilters.h:386
Belle2::ThreeHitFilters::fullAngleRZ
double fullAngleRZ()
calculates the angle between the hits/vectors (RZ), returning unit: angle in degrees
Definition: ThreeHitFilters.h:162
Belle2::ThreeHitFilters::fullAngle2D
double fullAngle2D(TVector3 &vecA, TVector3 &vecB)
calculates the angle between the hits/vectors (2D), generalized, returning unit: angle in radians WAR...
Definition: ThreeHitFilters.h:282
Belle2::ThreeHitFilters::checkCalcCircleCenter
void checkCalcCircleCenter()
checks whether the calcCircleCenter()-Member has been executed already and executes it if not
Definition: ThreeHitFilters.h:350
Belle2::ThreeHitFilters::calcCircleDist2IP
double calcCircleDist2IP()
calculates the distance of the point of closest approach of circle to the IP, returning unit: cm
Definition: ThreeHitFilters.h:174
Belle2::ThreeHitFilters::filterNan
double filterNan(double value)
returns zero if value is nan or inf
Definition: ThreeHitFilters.h:108
Belle2::ThreeHitFilters::calcHelixFit
double calcHelixFit()
reverse compatibility, calls calcHelixParameterFit
Definition: ThreeHitFilters.h:268
Belle2::ThreeHitFilters::fullAngleXY
double fullAngleXY()
calculates the angle between the hits/vectors (XY), returning unit: angle in degrees
Definition: ThreeHitFilters.h:142
Belle2::ThreeHitFilters::calcDeltaSlopeRZ
double calcDeltaSlopeRZ()
calculates deviations in the slope of the inner segment and the outer segment, returning unit: none
Definition: ThreeHitFilters.h:204
Belle2::TwoHitFilters::calcSlopeRZ
double calcSlopeRZ() const
calculates the angle of the slope of the hits in RZ, returnValue = theta = atan(r/z)
Definition: TwoHitFilters.h:80
Belle2::ThreeHitFilters::m_x2
double m_x2
internal intermediate value storing x^2, no enduser-relevance
Definition: ThreeHitFilters.h:391
Belle2::ThreeHitFilters::calcSign
int calcSign(TVector3 &a, TVector3 &b, TVector3 &c)
calculates calculates the sign of the curvature of given 3-hit-tracklet.
Definition: ThreeHitFilters.cc:21
Belle2::ThreeHitFilters::calcAngleXY
double calcAngleXY()
calculates the angle between the hits/vectors (XY), returning unit: none (calculation for degrees is ...
Definition: ThreeHitFilters.h:133
Belle2::ThreeHitFilters::calcAngle2D
double calcAngle2D(TVector3 &vecA, TVector3 &vecB)
calculates the angle between the hits/vectors (2D), generalized, returning unit: none.
Definition: ThreeHitFilters.h:273
Belle2::ThreeHitFilters::calcAngle3D
double calcAngle3D()
calculates the angle between the hits/vectors (3D), returning unit: none (calculation for degrees is ...
Definition: ThreeHitFilters.h:113
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ThreeHitFilters::ThreeHitFilters
ThreeHitFilters()
Empty constructor.
Definition: ThreeHitFilters.h:37
Belle2::ThreeHitFilters::m_hitA
TVector3 m_hitA
outer hit (position relevant for useful filter calculation) used for the filter calculation
Definition: ThreeHitFilters.h:396
Belle2::ThreeHitFilters::m_centerABC
TVector3 m_centerABC
center position of a circle in r-phi-plane formed by the 3 hits
Definition: ThreeHitFilters.h:395
Belle2::ThreeHitFilters::~ThreeHitFilters
~ThreeHitFilters()
Destructor.
Definition: ThreeHitFilters.h:72
Belle2::ThreeHitFilters::fullAngle3D
double fullAngle3D()
calculates the angle between the hits/vectors (3D), returning unit: angle in degrees
Definition: ThreeHitFilters.h:123
Belle2::ThreeHitFilters::m_radius
double m_radius
radius[cm] of a circle in r-phi-plane formed by the 3 hits
Definition: ThreeHitFilters.h:390
Belle2::ThreeHitFilters::m_vecAB
TVector3 m_vecAB
vector pointing from center hit to outer hit (outer segment)
Definition: ThreeHitFilters.h:399
Belle2::ThreeHitFilters::calcHelixParameterFit
double calcHelixParameterFit()
calculates the helixparameter describing the deviation in z per unit angle, returning unit: none
Definition: ThreeHitFilters.h:253
Belle2::ThreeHitFilters::getMagneticField
double getMagneticField()
returns the set value of the magnetic field in Tesla
Definition: ThreeHitFilters.h:103
Belle2::ThreeHitFilters::calcDeltaSlopeZOverS
double calcDeltaSlopeZOverS()
compares the "slopes" z over arc length.
Definition: ThreeHitFilters.h:237
Belle2::ThreeHitFilters::resetMagneticField
void resetMagneticField(double magneticFieldStrength=1.5)
Overrides Constructor-Setup for magnetic field.
Definition: ThreeHitFilters.h:98
Belle2::ThreeHitFilters::m_y2
double m_y2
internal intermediate value storing y^2, no enduser-relevance
Definition: ThreeHitFilters.h:392
Belle2::TwoHitFilters::resetValues
void resetValues(TVector3 &outerHit, TVector3 &innerHit)
Overrides Constructor-Setup.
Definition: TwoHitFilters.h:56
Belle2::ThreeHitFilters::calcDeltaSOverZ
double calcDeltaSOverZ()
calculates the helixparameter describing the deviation in arc length per unit in z.
Definition: ThreeHitFilters.h:220
Belle2::ThreeHitFilters::calcCircleCenter
void calcCircleCenter(TVector3 &a, TVector3 &b, TVector3 &c, TVector3 &circleCenter)
calculates an estimation of circleCenter position, result is written into the 4th input-parameter
Definition: ThreeHitFilters.h:300
Belle2::ThreeHitFilters::m_vecBC
TVector3 m_vecBC
vector pointing from inner hit to center hit (inner segment)
Definition: ThreeHitFilters.h:400
Belle2::ThreeHitFilters::m_circleCenterCalculated
bool m_circleCenterCalculated
initially set to false, will be set true if calcCircleCenter() is used at least once
Definition: ThreeHitFilters.h:388
Belle2::ThreeHitFilters::sanityCheckRadius
void sanityCheckRadius(double radius)
check Radius for bad values and throw exception if the value is bad
Definition: ThreeHitFilters.h:372
Belle2::ThreeHitFilters::checkCalcRadius
void checkCalcRadius()
checks whether the calcRadius()-Member has been executed already and executes it if not
Definition: ThreeHitFilters.h:360
Belle2::ThreeHitFilters::m_hitC
TVector3 m_hitC
inner hit (position relevant for useful filter calculation) used for the filter calculation
Definition: ThreeHitFilters.h:398
Belle2::TwoHitFilters::filterNan
double filterNan(double value) const
nice little nanChecker returns 0 if value was nan or inf, else returns value itself
Definition: TwoHitFilters.cc:17
Belle2::ThreeHitFilters::calcPt
double calcPt()
calculates the estimation of the transverse momentum of the 3-hit-tracklet, returning unit: GeV/c
Definition: ThreeHitFilters.h:184
Belle2::ThreeHitFilters::m_radiusCalculated
bool m_radiusCalculated
initially set to false, will be set true if calcInvCurvature() is used at least once
Definition: ThreeHitFilters.h:389