Belle II Software development
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
QualityEstimatorTripletFit Class Referencefinal
Packages » tracking

does a tripletFit of the given hits The filter is based on the paper 'A New Three-Dimensional Track Fit with Multiple Scattering' by Andre Schoening et al. More...

#include <QualityEstimatorTripletFit.h>

Inheritance diagram for QualityEstimatorTripletFit:
QualityEstimatorBase

Public Member Functions

virtual double estimateQuality (std::vector< SpacePoint const * > const &measurements) final
 Calculating the quality estimate using a triplet fit.
 
virtual QualityEstimationResults estimateQualityAndProperties (std::vector< SpacePoint const * > const &measurements) final
 perform quality estimation and extract additional information from the fit
 
void setMagneticFieldStrength (double magneticFieldZ=1.5)
 Setter for z component of magnetic field.
 

Protected Member Functions

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

Protected Attributes

std::vector< double > m_alphas
 angle alpha
 
std::vector< double > m_thetas
 angle theta
 
std::vector< double > m_R3Ds
 3D radius
 
std::vector< double > m_sigmaR3DSquareds
 squared error of 3D radius
 
double m_averageR3D = 0
 average 3D radius
 
double m_materialBudgetFactor = 1.45
 Triplet Fit hyper parameters.
 
double m_maxPt = 0.01
 Cut off value for 3D Radius calculated in Triplet Fit.
 
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.
 

Detailed Description

does a tripletFit of the given hits The filter is based on the paper 'A New Three-Dimensional Track Fit with Multiple Scattering' by Andre Schoening et al.

https://arxiv.org/abs/1606.04990

Definition at line 21 of file QualityEstimatorTripletFit.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.

80{ return m_magneticFieldZ * radius * 0.00299792458; }
Belle2::QualityEstimatorBase::m_magneticFieldZ
double m_magneticFieldZ
Member storing the z component of the magnetic field.
Definition: QualityEstimatorBase.h:84

◆ estimateQuality()

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

Calculating the quality estimate using a triplet fit.

Parameters
measurements: vector of SPs of the track candidate
Returns
quality indicator value

Using average material budged of SVD sensors for approximation of radiation length Belle II TDR page 156 states a value of 0.57% X_0. This approximation is a first approach to the problem and must be checked.

Implements QualityEstimatorBase.

Definition at line 22 of file QualityEstimatorTripletFit.cc.

23{
24 const int nTriplets = measurements.size() - 2;
25
26 if (nTriplets < 1) return 0;
27
28 double combinedChi2 = 0.;
29
30 // values required for pt calculation
31 m_R3Ds.clear();
32 m_R3Ds.reserve(nTriplets);
33 m_sigmaR3DSquareds.clear();
34 m_sigmaR3DSquareds.reserve(nTriplets);
35 m_alphas.clear();
36 m_alphas.reserve(nTriplets);
37 m_thetas.clear();
38 m_thetas.reserve(nTriplets);
39
40 // looping over all triplets
41 for (int i = 0; i < nTriplets; i++) {
42
43 // Three hits relevant for curent triplet
44 const B2Vector3D& hit0 = measurements.at(i)->getPosition();
45 const B2Vector3D& hit1 = measurements.at(i + 1)->getPosition();
46 const B2Vector3D& hit2 = measurements.at(i + 2)->getPosition();
47
48 const double d01sq = pow(hit1.X() - hit0.X(), 2) + pow(hit1.Y() - hit0.Y(), 2);
49 const double d12sq = pow(hit2.X() - hit1.X(), 2) + pow(hit2.Y() - hit1.Y(), 2);
50 const double d02sq = pow(hit2.X() - hit0.X(), 2) + pow(hit2.Y() - hit0.Y(), 2);
51
52 const double d01 = sqrt(d01sq);
53 const double d12 = sqrt(d12sq);
54 const double d02 = sqrt(d02sq);
55
56 const double z01 = hit1.Z() - hit0.Z();
57 const double z12 = hit2.Z() - hit1.Z();
58
59 // equation (8)
60 const double R_C = (d01 * d12 * d02) / sqrt(-d01sq * d01sq - d12sq * d12sq - d02sq * d02sq +
61 2 * d01sq * d12sq + 2 * d12sq * d02sq + 2 * d02sq * d01sq);
62
63 // equations (9)
64 const double Phi1C = 2. * asin(d01 / (2. * R_C));
65 const double Phi2C = 2. * asin(d12 / (2. * R_C));
66 // TODO Phi1C and Phi2C have 2 solutions (<Pi and >Pi), each, of which the correct one must be chosen!
67
68 // equations (10)
69 const double R3D1C = sqrt(R_C * R_C + (z01 * z01) / (Phi1C * Phi1C));
70 const double R3D2C = sqrt(R_C * R_C + (z12 * z12) / (Phi2C * Phi2C));
71
72 // equations (11)
73 const double theta1C = acos(z01 / (Phi1C * R3D1C));
74 const double theta2C = acos(z12 / (Phi2C * R3D2C));
75 const double theta = (theta1C + theta2C) / 2.;
76
77 // equations (15)
78 double alpha1 = R_C * R_C * Phi1C * Phi1C + z01 * z01;
79 alpha1 *= 1. / (0.5 * R_C * R_C * Phi1C * Phi1C * Phi1C / tan(Phi1C / 2.) + z01 * z01);
80 double alpha2 = R_C * R_C * Phi2C * Phi2C + z12 * z12;
81 alpha2 *= 1. / (0.5 * R_C * R_C * Phi2C * Phi2C * Phi2C / tan(Phi2C / 2.) + z12 * z12);
82
83 // equation (18)
84 const double PhiTilde = - 0.5 * (Phi1C * alpha1 + Phi2C * alpha2);
85 // equation (19)
86 const double eta = 0.5 * (Phi1C * alpha1 / R3D1C + Phi2C * alpha2 / R3D2C);
87
88 const double tanTheta1C = tan(theta1C);
89 const double tanTheta2C = tan(theta2C);
90 // equation (21)
91 const double ThetaTilde = theta2C - theta1C - (1 - alpha2) / tanTheta2C + (1 - alpha1) / tanTheta1C;
92 // equation (22)
93 const double beta = (1 - alpha2) / (R3D2C * tanTheta2C) - (1 - alpha1) / (R3D1C * tanTheta1C);
94
95 // Calculation of sigmaMS
96 // First get the orientation of the sensor plate for material budged calculation
97 double entranceAngle = theta;
98 VxdID::baseType sensorID = measurements.at(i + 1)->getVxdID();
99 int detID = measurements.at(i + 1)->getType();
100 if (sensorID != 0 and detID == VXD::SensorInfoBase::SVD) {
101 const SVD::SensorInfo& sensor = dynamic_cast<const SVD::SensorInfo&>(VXD::GeoCache::getInstance().getSensorInfo(sensorID));
102 const B2Vector3D& sensorOrigin = sensor.pointToGlobal(B2Vector3D(0, 0, 0), true);
103 const B2Vector3D& sensoru = sensor.pointToGlobal(B2Vector3D(1, 0, 0), true);
104 const B2Vector3D& sensorv = sensor.pointToGlobal(B2Vector3D(0, 1, 0), true);
105
106 B2Vector3D globalu = sensoru - sensorOrigin;
107 B2Vector3D globalv = sensorv - sensorOrigin;
108 B2Vector3D normal = globalu.Cross(globalv);
109
110 // Calculate the angle of incidence for the middle hit
111 if (sensorOrigin.Angle(normal) > M_PI * 0.5) { normal *= -1.; }
112 entranceAngle = (hit1 - hit0).Angle(normal);
113 }
114
119 const double XoverX0 = m_materialBudgetFactor * 0.0057 / cos(entranceAngle);
120
121 const double sinTheta = sin(theta);
122 // equation (23)
123 double R3D = - (eta * PhiTilde * sinTheta * sinTheta + beta * ThetaTilde);
124 R3D *= 1. / (eta * eta * sinTheta * sinTheta + beta * beta);
125 const double b = 4.5 / m_magneticFieldZ * sqrt(XoverX0);
126
127 // Calculation of maximal 3D radius from a maximal p_t cut off value
128 // (which is a hyper parameter of the Triplet Fit QE) via conversion
129 // using the magnetic field at the origin and the speed of light in the
130 // default units cm/ns times 1e-4 due to the units of the magnetic field.
131 double R3DmaxCut = m_maxPt / (m_magneticFieldZ * 1e-4 * Const::speedOfLight);
132 double R3D_truncated = R3D > R3DmaxCut ? R3DmaxCut : R3D;
133 const double sigmaMS = b / R3D_truncated;
134
135 // equation (24)
136 double Chi2min = pow(beta * PhiTilde - eta * ThetaTilde, 2);
137 Chi2min *= 1. / (sigmaMS * sigmaMS * (eta * eta + beta * beta / (sinTheta * sinTheta)));
138
139 // equation (25)
140 double sigmaR3DSquared = sigmaMS * sigmaMS / (eta * eta * sinTheta * sinTheta + beta * beta);
141
142 // bias correction as proposed in section 2.4 of the paper describing this fit. (see above)
143 // equation (30)
144 double delta = (beta * PhiTilde - eta * ThetaTilde) / (eta * PhiTilde * sinTheta * sinTheta + beta * ThetaTilde);
145 if (8 * delta * delta * sinTheta * sinTheta <= 1) {
146 // equation (29), the first part is the same as in equation (23)
147 R3D *= (0.75 + sqrt(1 - 8 * delta * delta * sinTheta * sinTheta) / 4.);
148 }
149
150 // required for optional pt calculation
151 m_thetas.push_back(theta);
152 m_alphas.push_back((alpha1 + alpha2) / 2.);
153
154 // store values for combination
155 m_R3Ds.push_back(R3D);
156 m_sigmaR3DSquareds.push_back(sigmaR3DSquared);
157
158 // equation (39)
159 combinedChi2 += Chi2min;
160 }
161
162 // Calculate average R3D
163 double numerator = 0;
164 double denominator = 0;
165 // c.f. equation (40)
166 for (short i = 0; i < nTriplets; ++i) {
167 // What is meant:
168 // numerator += pow(m_R3Ds.at(i), 3) / m_sigmaR3DSquareds.at(i);
169 // denominator += pow(m_R3Ds.at(i), 2) / m_sigmaR3DSquareds.at(i);
170 // To save clock cycles, only calculate R^2 / sigma^2 once, and only multiply with the missing R
171 double addValue = m_R3Ds.at(i) * m_R3Ds.at(i) / m_sigmaR3DSquareds.at(i);
172 numerator += addValue * m_R3Ds.at(i);
173 denominator += addValue;
174 }
175 m_averageR3D = numerator / denominator;
176
177 // Compare individual R3Ds with average R3D to improve chi2 as presented at:
178 // Connecting the Dots, Vienna, Feb 2016 by A. Schoening
179 double globalCompatibilityOfR3Ds = 0;
180 for (short i = 0; i < nTriplets; ++i) {
181 globalCompatibilityOfR3Ds += pow(m_R3Ds.at(i) - m_averageR3D, 2) / m_sigmaR3DSquareds.at(i);
182 }
183
184 double const finalChi2 = combinedChi2 + globalCompatibilityOfR3Ds;
185
186 m_results.chiSquared = finalChi2;
187
188 return TMath::Prob(finalChi2, 2 * measurements.size() - 5);
189}
Belle2::B2Vector3::Z
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:435
Belle2::B2Vector3::Cross
B2Vector3< DataType > Cross(const B2Vector3< DataType > &p) const
Cross product.
Definition: B2Vector3.h:296
Belle2::B2Vector3::X
DataType X() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:431
Belle2::B2Vector3::Y
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:433
Belle2::Const::speedOfLight
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
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:90
Belle2::QualityEstimatorTripletFit::m_R3Ds
std::vector< double > m_R3Ds
3D radius
Definition: QualityEstimatorTripletFit.h:41
Belle2::QualityEstimatorTripletFit::m_materialBudgetFactor
double m_materialBudgetFactor
Triplet Fit hyper parameters.
Definition: QualityEstimatorTripletFit.h:52
Belle2::QualityEstimatorTripletFit::m_thetas
std::vector< double > m_thetas
angle theta
Definition: QualityEstimatorTripletFit.h:40
Belle2::QualityEstimatorTripletFit::m_maxPt
double m_maxPt
Cut off value for 3D Radius calculated in Triplet Fit.
Definition: QualityEstimatorTripletFit.h:57
Belle2::QualityEstimatorTripletFit::m_sigmaR3DSquareds
std::vector< double > m_sigmaR3DSquareds
squared error of 3D radius
Definition: QualityEstimatorTripletFit.h:42
Belle2::QualityEstimatorTripletFit::m_alphas
std::vector< double > m_alphas
angle alpha
Definition: QualityEstimatorTripletFit.h:39
Belle2::SVD::SensorInfo
Specific implementation of SensorInfo for SVD Sensors which provides additional sensor specific infor...
Definition: SensorInfo.h:25
Belle2::VXD::GeoCache::getSensorInfo
const SensorInfoBase & getSensorInfo(Belle2::VxdID id) const
Return a referecne to the SensorInfo of a given SensorID.
Definition: GeoCache.cc:67
Belle2::VXD::GeoCache::getInstance
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition: GeoCache.cc:214
Belle2::VxdID::baseType
unsigned short baseType
The base integer type for VxdID.
Definition: VxdID.h:36
Belle2::B2Vector3D
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
Belle2::B2Vector3::at
DataType at(unsigned i) const
safe member access (with boundary check!)
Definition: B2Vector3.h:751
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2::QualityEstimationResults::chiSquared
std::optional< double > chiSquared
chi squared value obtained by the fit of the QE
Definition: QualityEstimatorBase.h:26

◆ estimateQualityAndProperties()

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

perform quality estimation and extract additional information from the fit

Parameters
measurements: vector of SPs of the track candidate
Returns
: QualityEstimationResults struct

Reimplemented from QualityEstimatorBase.

Definition at line 192 of file QualityEstimatorTripletFit.cc.

194{
195 // calculate and store chiSquared(calcChiSquared) and curvature(calcCurvature) in m_reuslts.
196 QualityEstimatorBase::estimateQualityAndProperties(measurements);
197 if (measurements.size() < 3) return m_results;
198
199 m_results.curvatureSign = calcCurvatureSignum(measurements);
200
201 // calculate pt and pt_sigma
202 double averageThetaPrime = 0;
203 // equation (43)
204 for (unsigned short i = 0; i < m_thetas.size(); ++i) {
205 averageThetaPrime += m_thetas.at(i) - (m_averageR3D - m_R3Ds.at(i)) * (1 - m_alphas.at(i)) / (tan(m_thetas.at(i)) * m_R3Ds.at(i));
206 }
207 averageThetaPrime /= m_thetas.size();
208
209 m_results.pt = calcPt(m_averageR3D * sin(averageThetaPrime));
210 m_results.pmag = calcPt(m_averageR3D);
211 return m_results;
212}
Belle2::QualityEstimatorBase::estimateQualityAndProperties
virtual QualityEstimationResults estimateQualityAndProperties(std::vector< SpacePoint const * > const &measurements)
Quality estimation providing additional quantities Calculates quality indicator in range [0,...
Definition: QualityEstimatorBase.h:68
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:80
Belle2::calcCurvatureSignum
short calcCurvatureSignum(std::vector< SpacePoint const * > const &measurements)
Calculate curvature based on triplets of measurements.
Definition: CalcCurvatureSignum.h:22
Belle2::QualityEstimationResults::pt
std::optional< double > pt
transverse momentum estimate from the QE
Definition: QualityEstimatorBase.h:29
Belle2::QualityEstimationResults::curvatureSign
std::optional< short > curvatureSign
direction of curvature as obtained by the QE
Definition: QualityEstimatorBase.h:27
Belle2::QualityEstimationResults::pmag
std::optional< double > pmag
momentum magnitute estimate from the QE
Definition: QualityEstimatorBase.h:30

◆ 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.

53{m_magneticFieldZ = magneticFieldZ;}

Member Data Documentation

◆ m_alphas

std::vector<double> m_alphas
protected

angle alpha

Definition at line 39 of file QualityEstimatorTripletFit.h.

◆ m_averageR3D

double m_averageR3D = 0
protected

average 3D radius

Definition at line 44 of file QualityEstimatorTripletFit.h.

◆ m_magneticFieldZ

double m_magneticFieldZ = 1.5
protectedinherited

Member storing the z component of the magnetic field.

Definition at line 84 of file QualityEstimatorBase.h.

◆ m_materialBudgetFactor

double m_materialBudgetFactor = 1.45
protected

Triplet Fit hyper parameters.

Scaling factor for material budget which is applied to the radiation length value X/X_0 = 0.57% which is taken from the Belle II TDR page 156. This scaling factor is optimized to achieve the best performance on MC.

Definition at line 52 of file QualityEstimatorTripletFit.h.

◆ m_maxPt

double m_maxPt = 0.01
protected

Cut off value for 3D Radius calculated in Triplet Fit.

This value is a hyper parameter of the Triplet Fit which is tuned on MC.

Definition at line 57 of file QualityEstimatorTripletFit.h.

◆ m_R3Ds

std::vector<double> m_R3Ds
protected

3D radius

Definition at line 41 of file QualityEstimatorTripletFit.h.

◆ 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.

◆ m_sigmaR3DSquareds

std::vector<double> m_sigmaR3DSquareds
protected

squared error of 3D radius

Definition at line 42 of file QualityEstimatorTripletFit.h.

◆ m_thetas

std::vector<double> m_thetas
protected

angle theta

Definition at line 40 of file QualityEstimatorTripletFit.h.

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