Belle II Software  release-05-01-25
QualityEstimatorCircleFit.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2017 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Felix Metzner, Jonas Wagner *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <tracking/trackFindingVXD/trackQualityEstimators/QualityEstimatorCircleFit.h>
12 #include <tracking/trackFindingVXD/utilities/CalcCurvatureSignum.h>
13 
14 #include <framework/logging/Logger.h>
15 
16 #include <TMath.h>
17 #include <cmath>
18 
19 using namespace Belle2;
20 
21 double QualityEstimatorCircleFit::estimateQuality(std::vector<SpacePoint const*> const& measurements)
22 {
23  if (measurements.size() < 4) {
24  if (measurements.size() == 3) return 0.2; // Arbitrary value to prevent excluding measurements with 3 hits later on.
25  else return 0;
26  }
27  // Calculates Curvature: True means clockwise, False means counterclockwise.
28  // TODO this is not an optimized approach; just to get things to work.
29  // CalcCurvature could be integrated into the looping over the hits which CircleFit does anyhow.
31  bool clockwise = *(m_results.curvatureSign) >= 0;
32 
33  double stopper = 0.000000001;
34  double meanX = 0, meanY = 0, meanX2 = 0, meanY2 = 0, meanR2 = 0, meanR4 = 0, meanXR2 = 0, meanYR2 = 0, meanXY = 0; //mean values
35  double sumWeights = 0, divisor/*, weightNormalizer = 0*/; // sumWeights is sum of weights, divisor is 1/sumWeights;
36  double tuningParameter = 1.; //0.02; // this parameter is for internal tuning of the weights, 1 means no influence of parameter.
37 
38  // looping over all hits and do the division afterwards
39  for (const SpacePoint* hit : measurements) {
40  double weight = 1. / (sqrt(hit->getPositionError().X() * hit->getPositionError().X() + hit->getPositionError().Y() *
41  hit->getPositionError().Y()) * tuningParameter);
42  sumWeights += weight;
43  if (std::isnan(weight) or std::isinf(weight)) {
44  B2ERROR("QualityEstimators::circleFit, chosen sigma is 'nan': " << weight << ", setting arbitrary error: "
45  << stopper << ")"); weight = stopper;
46  }
47  double x = hit->getPosition().X();
48  double y = hit->getPosition().Y();
49  double x2 = x * x;
50  double y2 = y * y;
51  double r2 = x2 + y2;
52  meanX += x * weight;
53  meanY += y * weight;
54  meanXY += x * y * weight;
55  meanX2 += x2 * weight;
56  meanY2 += y2 * weight;
57  meanXR2 += x * r2 * weight;
58  meanYR2 += y * r2 * weight;
59  meanR2 += r2 * weight;
60  meanR4 += r2 * r2 * weight;
61  }
62  divisor = 1. / sumWeights;
63  meanX *= divisor;
64  meanY *= divisor;
65  meanXY *= divisor;
66  meanY2 *= divisor;
67  meanX2 *= divisor;
68  meanXR2 *= divisor;
69  meanYR2 *= divisor;
70  meanR2 *= divisor;
71  meanR4 *= divisor;
72 
73  // covariances:
74  double covXX = meanX2 - meanX * meanX;
75  double covXY = meanXY - meanX * meanY;
76  double covYY = meanY2 - meanY * meanY;
77  double covXR2 = meanXR2 - meanX * meanR2;
78  double covYR2 = meanYR2 - meanY * meanR2;
79  double covR2R2 = meanR4 - meanR2 * meanR2;
80 
81  if (covR2R2 == 0) {
82  return 0; // TODO could be problematic if it is pretty near to 0
83  }
84 
85  // q1, q2: helping variables, to make the code more readable
86  double q1 = covR2R2 * covXY - covXR2 * covYR2;
87  double q2 = covR2R2 * (covXX - covYY) - covXR2 * covXR2 + covYR2 * covYR2;
88 
89  // physical meaning: phi value of the point of closest approach of the fitted circle to the origin
90  double pocaPhi = 0.5 * atan2(2. * q1, q2);
91 
92  double sinPhi = sin(pocaPhi);
93  double cosPhi = cos(pocaPhi);
94  double kappa = (sinPhi * covXR2 - cosPhi * covYR2) / covR2R2;
95  double delta = -kappa * meanR2 + sinPhi * meanX - cosPhi * meanY;
96  double rootTerm = sqrt(1. - 4.*delta * kappa);
97  // rho = curvature in X-Y-plane = 1/radius of fitting circle, used for pT-calculation
98  double curvature = 2.*kappa / (rootTerm);
99  double pocaD = 2.*delta / (1. + rootTerm);
100 
101  // Checks if the random Curvature of CircleFit corresponds to CalcCurvature and adjust the results accordingly.
102  if ((curvature < 0 && clockwise) || (curvature > 0 && !clockwise)) {
103  // this is according to eq. 23 in the paper of Karimäki
104  curvature = -curvature;
105  pocaD = -pocaD;
106  //TODO ..and swap correlation Terms V_rho_phi and V_rho_d (which are not implemented anyway)
107  }
108  double radius = 1. / curvature;
109  double absRadius = fabs(radius);
110 
111  m_results.pt = calcPt(absRadius);
112 
113  double chi2 = sumWeights * (1. + curvature * pocaD) * (1. + curvature * pocaD)
114  * (sinPhi * sinPhi * covXX - 2.*sinPhi * cosPhi * covXY + cosPhi * cosPhi * covYY - kappa * kappa * covR2R2);
115 
116  m_results.chiSquared = chi2;
117 
118  return TMath::Prob(chi2, measurements.size() - 3);
119 }
120 
Belle2::QualityEstimationResults::chiSquared
boost::optional< double > chiSquared
chi squared value obtained by the fit of the QE
Definition: QualityEstimatorBase.h:38
Belle2::calcCurvatureSignum
short calcCurvatureSignum(std::vector< SpacePoint const * > const &measurements)
Calculate curvature based on triplets of measurements.
Definition: CalcCurvatureSignum.h:32
Belle2::QualityEstimationResults::pt
boost::optional< double > pt
transverse momentum estimate from the QE
Definition: QualityEstimatorBase.h:40
Belle2::QualityEstimatorCircleFit::estimateQuality
virtual double estimateQuality(std::vector< SpacePoint const * > const &measurements) final
Minimal implementation of the quality estimation Calculates quality indicator in range [0,...
Definition: QualityEstimatorCircleFit.cc:21
Belle2::SpacePoint
SpacePoint typically is build from 1 PXDCluster or 1-2 SVDClusters.
Definition: SpacePoint.h:52
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::QualityEstimatorBase::calcPt
double calcPt(double const radius) const
Returns a value for the transverse momentum in GeV calculated from a provided radius.
Definition: QualityEstimatorBase.h:91
Belle2::QualityEstimationResults::curvatureSign
boost::optional< short > curvatureSign
direction of curvature as obtained by the QE
Definition: QualityEstimatorBase.h:39
Belle2::QualityEstimatorBase::m_results
QualityEstimationResults m_results
Result of the quality estimation This is stored as a member variable, because some values may be calc...
Definition: QualityEstimatorBase.h:101