Belle II Software  release-08-01-10
QualityEstimatorCircleFit Class Reference

Class containing the algorithm to perform the simple circle fit. More...

#include <QualityEstimatorCircleFit.h>

Inheritance diagram for QualityEstimatorCircleFit:
Collaboration diagram for QualityEstimatorCircleFit:

Public Member Functions

virtual double estimateQuality (std::vector< SpacePoint const * > const &measurements) final
 Minimal implementation of the quality estimation Calculates quality indicator in range [0,1]. More...
 
void setMagneticFieldStrength (double magneticFieldZ=1.5)
 Setter for z component of magnetic field. More...
 
virtual QualityEstimationResults estimateQualityAndProperties (std::vector< SpacePoint const * > const &measurements)
 Quality estimation providing additional quantities Calculates quality indicator in range [0,1] Optionally returns chi2 and additional informations. More...
 

Protected Member Functions

double calcPt (double const radius) const
 Returns a value for the transverse momentum in GeV calculated from a provided radius. More...
 

Protected Attributes

double m_magneticFieldZ = 1.5
 Member storing the z component of the magnetic field.
 
QualityEstimationResults m_results
 Result of the quality estimation This is stored as a member variable, because some values may be calculated by 'estimateQuality' anyways. More...
 

Detailed Description

Class containing the algorithm to perform the simple circle fit.

Definition at line 19 of file QualityEstimatorCircleFit.h.

Member Function Documentation

◆ calcPt()

double calcPt ( double const  radius) const
inlineprotectedinherited

Returns a value for the transverse momentum in GeV calculated from a provided radius.

Utilizing m_magneticFieldZ and hard coded speed of light

Definition at line 80 of file QualityEstimatorBase.h.

◆ estimateQuality()

double estimateQuality ( std::vector< SpacePoint const * > const &  measurements)
finalvirtual

Minimal implementation of the quality estimation Calculates quality indicator in range [0,1].

measurements - std::vector<SpacePoint const*> ordered from innermost to outermost measurement

WARNING hardcoded values!

Implements QualityEstimatorBase.

Definition at line 19 of file QualityEstimatorCircleFit.cc.

20 {
21  if (measurements.size() < 4) {
22  if (measurements.size() == 3) return 0.2; // Arbitrary value to prevent excluding measurements with 3 hits later on.
23  else return 0;
24  }
25  // Calculates Curvature: True means clockwise, False means counterclockwise.
26  // TODO this is not an optimized approach; just to get things to work.
27  // CalcCurvature could be integrated into the looping over the hits which CircleFit does anyhow.
29  bool clockwise = *(m_results.curvatureSign) >= 0;
30 
31  double stopper = 0.000000001;
32  double meanX = 0, meanY = 0, meanX2 = 0, meanY2 = 0, meanR2 = 0, meanR4 = 0, meanXR2 = 0, meanYR2 = 0, meanXY = 0; //mean values
33  double sumWeights = 0, divisor/*, weightNormalizer = 0*/; // sumWeights is sum of weights, divisor is 1/sumWeights;
34  double tuningParameter = 1.; //0.02; // this parameter is for internal tuning of the weights, 1 means no influence of parameter.
35 
36  // looping over all hits and do the division afterwards
37  for (const SpacePoint* hit : measurements) {
38  double weight = 1. / (sqrt(hit->getPositionError().X() * hit->getPositionError().X() + hit->getPositionError().Y() *
39  hit->getPositionError().Y()) * tuningParameter);
40  sumWeights += weight;
41  if (std::isnan(weight) or std::isinf(weight)) {
42  B2ERROR("QualityEstimators::circleFit, chosen sigma is 'nan': " << weight << ", setting arbitrary error: "
43  << stopper << ")"); weight = stopper;
44  }
45  double x = hit->getPosition().X();
46  double y = hit->getPosition().Y();
47  double x2 = x * x;
48  double y2 = y * y;
49  double r2 = x2 + y2;
50  meanX += x * weight;
51  meanY += y * weight;
52  meanXY += x * y * weight;
53  meanX2 += x2 * weight;
54  meanY2 += y2 * weight;
55  meanXR2 += x * r2 * weight;
56  meanYR2 += y * r2 * weight;
57  meanR2 += r2 * weight;
58  meanR4 += r2 * r2 * weight;
59  }
60  divisor = 1. / sumWeights;
61  meanX *= divisor;
62  meanY *= divisor;
63  meanXY *= divisor;
64  meanY2 *= divisor;
65  meanX2 *= divisor;
66  meanXR2 *= divisor;
67  meanYR2 *= divisor;
68  meanR2 *= divisor;
69  meanR4 *= divisor;
70 
71  // covariances:
72  double covXX = meanX2 - meanX * meanX;
73  double covXY = meanXY - meanX * meanY;
74  double covYY = meanY2 - meanY * meanY;
75  double covXR2 = meanXR2 - meanX * meanR2;
76  double covYR2 = meanYR2 - meanY * meanR2;
77  double covR2R2 = meanR4 - meanR2 * meanR2;
78 
79  if (covR2R2 == 0) {
80  return 0; // TODO could be problematic if it is pretty near to 0
81  }
82 
83  // q1, q2: helping variables, to make the code more readable
84  double q1 = covR2R2 * covXY - covXR2 * covYR2;
85  double q2 = covR2R2 * (covXX - covYY) - covXR2 * covXR2 + covYR2 * covYR2;
86 
87  // physical meaning: phi value of the point of closest approach of the fitted circle to the origin
88  double pocaPhi = 0.5 * atan2(2. * q1, q2);
89 
90  double sinPhi = sin(pocaPhi);
91  double cosPhi = cos(pocaPhi);
92  double kappa = (sinPhi * covXR2 - cosPhi * covYR2) / covR2R2;
93  double delta = -kappa * meanR2 + sinPhi * meanX - cosPhi * meanY;
94  double rootTerm = sqrt(1. - 4.*delta * kappa);
95  // rho = curvature in X-Y-plane = 1/radius of fitting circle, used for pT-calculation
96  double curvature = 2.*kappa / (rootTerm);
97  double pocaD = 2.*delta / (1. + rootTerm);
98 
99  // Checks if the random Curvature of CircleFit corresponds to CalcCurvature and adjust the results accordingly.
100  if ((curvature < 0 && clockwise) || (curvature > 0 && !clockwise)) {
101  // this is according to eq. 23 in the paper of Karimäki
102  curvature = -curvature;
103  pocaD = -pocaD;
104  //TODO ..and swap correlation Terms V_rho_phi and V_rho_d (which are not implemented anyway)
105  }
106  double radius = 1. / curvature;
107  double absRadius = fabs(radius);
108  double pt = calcPt(absRadius);
109  m_results.pt = pt;
110  m_results.pocaD = pocaD;
111 
112  // Estimating theta of the track with the theta of the innermost hit.
113  // Otherwise the track is solely a circle without theta information.
114  double innermostHitTheta = 0.;
115  double innermostHitRadiusSquared = 100000000;
116  for (const SpacePoint* hit : measurements) {
117  if (hit->getPosition().Perp2() < innermostHitRadiusSquared) {
118  innermostHitRadiusSquared = hit->getPosition().Perp2();
119  innermostHitTheta = hit->getPosition().Theta();
120  }
121  }
122  // Account for precision when checking equality of innermostHitTheta with 0
123  double pz = (innermostHitTheta <= 1e-6 ? 0 : pt / tan(innermostHitTheta));
124  // Both p and pmag are only rough estimates based and by far not accutate!
125  m_results.p = B2Vector3D(pt * cosPhi, pt * sinPhi, pz);
126  m_results.pmag = sqrt(pt * pt + pz * pz);
127 
128  double chi2 = sumWeights * (1. + curvature * pocaD) * (1. + curvature * pocaD)
129  * (sinPhi * sinPhi * covXX - 2.*sinPhi * cosPhi * covXY + cosPhi * cosPhi * covYY - kappa * kappa * covR2R2);
130 
131  m_results.chiSquared = chi2;
132 
133  return TMath::Prob(chi2, measurements.size() - 3);
134 }
QualityEstimationResults m_results
Result of the quality estimation This is stored as a member variable, because some values may be calc...
double calcPt(double const radius) const
Returns a value for the transverse momentum in GeV calculated from a provided radius.
SpacePoint typically is build from 1 PXDCluster or 1-2 SVDClusters.
Definition: SpacePoint.h:42
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
double tan(double a)
tan for double
Definition: beamHelpers.h:31
short calcCurvatureSignum(std::vector< SpacePoint const * > const &measurements)
Calculate curvature based on triplets of measurements.
std::optional< double > pt
transverse momentum estimate from the QE
std::optional< short > curvatureSign
direction of curvature as obtained by the QE
std::optional< double > pocaD
distance to the z-axis of the POCA
std::optional< double > chiSquared
chi squared value obtained by the fit of the QE
std::optional< B2Vector3D > p
momentum vector estimate from the QE
std::optional< double > pmag
momentum magnitute estimate from the QE

◆ estimateQualityAndProperties()

virtual QualityEstimationResults estimateQualityAndProperties ( std::vector< SpacePoint const * > const &  measurements)
inlinevirtualinherited

Quality estimation providing additional quantities Calculates quality indicator in range [0,1] Optionally returns chi2 and additional informations.

Eg. momentum estimation.

measurements - std::vector<SpacePoint const*> ordered from innermost to outermost measurement

Reimplemented in QualityEstimatorMC, and QualityEstimatorTripletFit.

Definition at line 68 of file QualityEstimatorBase.h.

◆ setMagneticFieldStrength()

void setMagneticFieldStrength ( double  magneticFieldZ = 1.5)
inlineinherited

Setter for z component of magnetic field.

Parameters
magneticFieldZ: value to set it to

Definition at line 53 of file QualityEstimatorBase.h.

Member Data Documentation

◆ m_results

QualityEstimationResults m_results
protectedinherited

Result of the quality estimation This is stored as a member variable, because some values may be calculated by 'estimateQuality' anyways.

Therefore they don't need to be calculated explicitly in 'estimateQualityAndProperties'.

Definition at line 90 of file QualityEstimatorBase.h.


The documentation for this class was generated from the following files: