Belle II Software development
QualityEstimatorLineFit3D.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include "tracking/trackFindingVXD/trackQualityEstimators/QualityEstimatorLineFit3D.h"
10#include <math.h>
11#include <TMath.h>
12
13using namespace Belle2;
14
15double QualityEstimatorLineFit3D::estimateQuality(std::vector<SpacePoint const*> const& measurements)
16{
17 double sumWyi = 0, // sum of weights for Yi
18 sumWzi = 0, // sum of weights for Zi
19 sumWyiXi = 0, // sum of (y-weights times x-values)
20 sumWziXi = 0, // sum of (z-weights times x-values)
21 sumWyiYi = 0, // sum of (y-weights times y-values)
22 sumWziZi = 0, // sum of (z-weights times z-values)
23 sumWyiXiYi = 0, // sum of (y-weights times x-values times y-values)
24 sumWziXiZi = 0, // sum of (z-weights times x-values times z-values)
25 sumWyiXi2 = 0, // sum of (y-weights times x-values^2)
26 sumWziXi2 = 0, // sum of (z-weights times x-values^2)
27 detValY = 0, // determinant for norming values - y
28 detValZ = 0, // determinant for norming values - z
29 slopeY = 0, // = a of model
30 slopeZ = 0, // = c of model
31 chi2 = 0, // final chi2-value of fit
32 interceptY = 0, // b of model, needed only for chi2-calculation
33 interceptZ = 0; // d of model, needed only for chi2-calculation
34
35 // NOTE: this approach is not optimal. Maybe can be optimized for less redundancy
36 for (const SpacePoint* aHit : measurements) {
37 double Wyi = (1. / (aHit->getPositionError().Y() * aHit->getPositionError().Y()));
38 double Wzi = (1. / (aHit->getPositionError().Z() * aHit->getPositionError().Z()));
39
40 sumWyi += Wyi;
41 sumWzi += Wzi;
42
43 sumWyiXi += Wyi * aHit->getPosition().X();
44 sumWziXi += Wzi * aHit->getPosition().X();
45
46 sumWyiYi += Wyi * aHit->getPosition().Y();
47 sumWziZi += Wzi * aHit->getPosition().Z();
48
49 sumWyiXiYi += Wyi * aHit->getPosition().X() * aHit->getPosition().Y();
50 sumWziXiZi += Wzi * aHit->getPosition().X() * aHit->getPosition().Z();
51
52 sumWyiXi2 += Wyi * aHit->getPosition().X() * aHit->getPosition().X();
53 sumWziXi2 += Wzi * aHit->getPosition().X() * aHit->getPosition().X();
54 }
55
56 detValY = sumWyiXi2 * sumWyi - sumWyiXi * sumWyiXi;
57 if (detValY == 0) {
58 return 0;
59 }
60 detValY = 1. / detValY; // invert
61
62 detValZ = sumWziXi2 * sumWzi - sumWziXi * sumWziXi;
63 if (detValZ == 0) {
64 return 0;
65 }
66 detValZ = 1. / detValZ; // invert
67
68 slopeY = detValY * (sumWyi * sumWyiXiYi - sumWyiXi * sumWyiYi);
69 slopeZ = detValZ * (sumWzi * sumWziXiZi - sumWziXi * sumWziZi);
70
71 interceptY = detValY * (- sumWyiXi * sumWyiXiYi + sumWyiXi2 * sumWyiYi);
72 interceptZ = detValZ * (- sumWziXi * sumWziXiZi + sumWziXi2 * sumWziZi);
73
74 for (const SpacePoint* aHit : measurements) { // chi2 of xy-fit and of xz-fit can be combined by adding their values
75 chi2 += pow(((aHit->getPosition().Y() - slopeY * aHit->getPosition().X() - interceptY) / aHit->getPositionError().Y()), 2)
76 + pow(((aHit->getPosition().Z() - slopeZ * aHit->getPosition().X() - interceptZ) / aHit->getPositionError().Z()), 2);
77 }
78 m_results.chiSquared = chi2;
79
80 return TMath::Prob(chi2, measurements.size() - 1);
81}
82
QualityEstimationResults m_results
Result of the quality estimation This is stored as a member variable, because some values may be calc...
virtual double estimateQuality(std::vector< SpacePoint const * > const &measurements) final
Minimal implementation of the quality estimation Calculates quality indicator in range [0,...
SpacePoint typically is build from 1 PXDCluster or 1-2 SVDClusters.
Definition: SpacePoint.h:42
Abstract base class for different kinds of events.
std::optional< double > chiSquared
chi squared value obtained by the fit of the QE