Belle II Software  release-05-01-25
QualityEstimatorTripletFit.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 <math.h>
12 #include <TMath.h>
13 
14 #include <framework/gearbox/Const.h>
15 
16 #include <svd/geometry/SensorInfo.h>
17 
18 #include "tracking/trackFindingVXD/trackQualityEstimators/QualityEstimatorTripletFit.h"
19 #include <tracking/trackFindingVXD/utilities/CalcCurvatureSignum.h>
20 
21 
22 using namespace Belle2;
23 
24 double QualityEstimatorTripletFit::estimateQuality(std::vector<SpacePoint const*> const& measurements)
25 {
26  const int nTriplets = measurements.size() - 2;
27 
28  if (nTriplets < 1) return 0;
29 
30  double combinedChi2 = 0.;
31 
32  // values required for pt calculation
33  m_R3Ds.clear();
34  m_R3Ds.reserve(nTriplets);
35  m_sigmaR3DSquareds.clear();
36  m_sigmaR3DSquareds.reserve(nTriplets);
37  m_alphas.clear();
38  m_alphas.reserve(nTriplets);
39  m_thetas.clear();
40  m_thetas.reserve(nTriplets);
41 
42  // looping over all triplets
43  for (int i = 0; i < nTriplets; i++) {
44 
45  // Three hits relevant for curent triplet
46  const TVector3 hit0 = measurements.at(i)->getPosition();
47  const TVector3 hit1 = measurements.at(i + 1)->getPosition();
48  const TVector3 hit2 = measurements.at(i + 2)->getPosition();
49 
50  const double d01sq = pow(hit1.X() - hit0.X(), 2) + pow(hit1.Y() - hit0.Y(), 2);
51  const double d12sq = pow(hit2.X() - hit1.X(), 2) + pow(hit2.Y() - hit1.Y(), 2);
52  const double d02sq = pow(hit2.X() - hit0.X(), 2) + pow(hit2.Y() - hit0.Y(), 2);
53 
54  const double d01 = sqrt(d01sq);
55  const double d12 = sqrt(d12sq);
56  const double d02 = sqrt(d02sq);
57 
58  const double z01 = hit1.Z() - hit0.Z();
59  const double z12 = hit2.Z() - hit1.Z();
60 
61  const double R_C = (d01 * d12 * d02) / sqrt(-d01sq * d01sq - d12sq * d12sq - d02sq * d02sq + 2 * d01sq * d12sq + 2 * d12sq * d02sq +
62  2 *
63  d02sq * d01sq);
64 
65  const double Phi1C = 2. * asin(d01 / (2. * R_C));
66  const double Phi2C = 2. * asin(d12 / (2. * R_C));
67  // TODO Phi1C and Phi2C have 2 solutions (<Pi and >Pi), each, of which the correct one must be chosen!
68 
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  const double theta1C = acos(z01 / (Phi1C * R3D1C));
73  const double theta2C = acos(z12 / (Phi2C * R3D2C));
74  const double theta = (theta1C + theta2C) / 2.;
75 
76  double alpha1 = R_C * R_C * Phi1C * Phi1C + z01 * z01;
77  alpha1 *= 1. / (0.5 * R_C * R_C * Phi1C * Phi1C * Phi1C / tan(Phi1C / 2.) + z01 * z01);
78  double alpha2 = R_C * R_C * Phi2C * Phi2C + z12 * z12;
79  alpha2 *= 1. / (0.5 * R_C * R_C * Phi2C * Phi2C * Phi2C / tan(Phi2C / 2.) + z12 * z12);
80 
81  const double PhiTilde = - 0.5 * (Phi1C * alpha1 + Phi2C * alpha2);
82  const double eta = 0.5 * Phi1C * alpha1 / R3D1C + 0.5 * Phi2C * alpha2 / R3D2C;
83  const double ThetaTilde = theta2C - theta1C - (1 - alpha2) / tan(theta2C) + (1 - alpha1) / tan(theta1C);
84  const double beta = (1 - alpha2) / (R3D2C * tan(theta2C)) - (1 - alpha1) / (R3D1C * tan(theta1C));
85 
86  // Calculation of sigmaMS
87  // First get the orientation of the sensor plate for material budged calculation
88  double entranceAngle = theta;
89  VxdID::baseType sensorID = measurements.at(i + 1)->getVxdID();
90  int detID = measurements.at(i + 1)->getType();
91  if (sensorID != 0 and detID == VXD::SensorInfoBase::SVD) {
92  const SVD::SensorInfo& sensor = dynamic_cast<const SVD::SensorInfo&>(VXD::GeoCache::get(sensorID));
93  TVector3 sensorOrigin = sensor.pointToGlobal(TVector3(0, 0, 0), true);
94  TVector3 sensoru = sensor.pointToGlobal(TVector3(1, 0, 0), true);
95  TVector3 sensorv = sensor.pointToGlobal(TVector3(0, 1, 0), true);
96 
97  TVector3 globalu = sensoru - sensorOrigin;
98  TVector3 globalv = sensorv - sensorOrigin;
99  TVector3 normal = globalu.Cross(globalv);
100 
101  // Calculate the angle of incidence for the middle hit
102  if (sensorOrigin.Angle(normal) > M_PI * 0.5) { normal *= -1.; }
103  entranceAngle = (hit1 - hit0).Angle(normal);
104  }
105 
110  const double XoverX0 = m_materialBudgetFactor * 0.0057 / cos(entranceAngle);
111 
112  double R3D = - (eta * PhiTilde * sin(theta) * sin(theta) + beta * ThetaTilde);
113  R3D *= 1. / (eta * eta * sin(theta) * sin(theta) + beta * beta);
114  const double b = 4.5 / m_magneticFieldZ * sqrt(XoverX0);
115 
116  // Calculation of maximal 3D radius from a maximal p_t cut off value
117  // (which is a hyper parameter of the Triplet Fit QE) via conversion
118  // using the magnetic field at the origin and the speed of light in the
119  // default units cm/ns times 1e-4 due to the units of the magnetic field.
120  double R3DmaxCut = m_maxPt / (m_magneticFieldZ * 1e-4 * Const::speedOfLight);
121  double R3D_truncated = R3D > R3DmaxCut ? R3DmaxCut : R3D;
122  const double sigmaMS = b / R3D_truncated;
123 
124  double sigmaR3DSquared = pow(sigmaMS, 2) / (pow(eta * sin(theta), 2) + pow(beta, 2));
125 
126  double Chi2min = pow(beta * PhiTilde - eta * ThetaTilde, 2);
127  Chi2min *= 1. / (sigmaMS * sigmaMS * (eta * eta + beta * beta / pow(sin(theta), 2)));
128 
129  // bias correction as proposed in section 2.4 of the paper describing this fit. (see above)
130  double delta = (beta * PhiTilde - eta * ThetaTilde) / (eta * PhiTilde * sin(theta) * sin(theta) + beta * ThetaTilde);
131  if (8 * delta * delta * sin(theta) * sin(theta) <= 1) {
132  R3D *= (0.75 + sqrt(1 - 8 * delta * delta * sin(theta) * sin(theta)) / 4.);
133  }
134 
135  // required for optional pt calculation
136  m_thetas.push_back(theta);
137  m_alphas.push_back((alpha1 + alpha2) / 2.);
138 
139  // store values for combination
140  m_R3Ds.push_back(R3D);
141  m_sigmaR3DSquareds.push_back(sigmaR3DSquared);
142 
143  combinedChi2 += Chi2min;
144  }
145 
146  // Calculate average R3D
147  double numerator = 0;
148  double denominator = 0;
149  for (short i = 0; i < nTriplets; ++i) {
150  numerator += pow(m_R3Ds.at(i), 3) / m_sigmaR3DSquareds.at(i);
151  denominator += pow(m_R3Ds.at(i), 2) / m_sigmaR3DSquareds.at(i);
152  }
153  m_averageR3D = numerator / denominator;
154 
155  // Compare individual R3Ds with average R3D to improve chi2 as presented at:
156  // Connecting the Dots, Vienna, Feb 2016 by A. Schoening
157  double globalCompatibilityOfR3Ds = 0;
158  for (short i = 0; i < nTriplets; ++i) {
159  globalCompatibilityOfR3Ds += pow(m_R3Ds.at(i) - m_averageR3D, 2) / m_sigmaR3DSquareds.at(i);
160  }
161 
162  double const finalChi2 = combinedChi2 + globalCompatibilityOfR3Ds;
163 
164  m_results.chiSquared = finalChi2;
165 
166  return TMath::Prob(finalChi2, 2 * measurements.size() - 5);
167 }
168 
169 
171  measurements)
172 {
173  // calculate and store chiSquared(calcChiSquared) and curvature(calcCurvature) in m_reuslts.
175  if (measurements.size() < 3) return m_results;
176 
178 
179  // calculate pt and pt_sigma
180  double averageThetaPrime = 0;
181  for (unsigned short i = 0; i < m_thetas.size(); ++i) {
182  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));
183  }
184  averageThetaPrime /= m_thetas.size();
185 
186  m_results.pt = calcPt(m_averageR3D * sin(averageThetaPrime));
188  return m_results;
189 }
190 
Belle2::QualityEstimatorTripletFit::m_thetas
std::vector< double > m_thetas
angle theta
Definition: QualityEstimatorTripletFit.h:50
Belle2::QualityEstimationResults
Container for complete fit/estimation results.
Definition: QualityEstimatorBase.h:36
Belle2::QualityEstimationResults::pmag
boost::optional< double > pmag
momentum magnitute estimate from the QE
Definition: QualityEstimatorBase.h:41
Belle2::QualityEstimationResults::chiSquared
boost::optional< double > chiSquared
chi squared value obtained by the fit of the QE
Definition: QualityEstimatorBase.h:38
Belle2::VXD::GeoCache::get
static const SensorInfoBase & get(Belle2::VxdID id)
Return a reference to the SensorInfo of a given SensorID.
Definition: GeoCache.h:141
Belle2::calcCurvatureSignum
short calcCurvatureSignum(std::vector< SpacePoint const * > const &measurements)
Calculate curvature based on triplets of measurements.
Definition: CalcCurvatureSignum.h:32
Belle2::QualityEstimatorTripletFit::m_alphas
std::vector< double > m_alphas
angle alpha
Definition: QualityEstimatorTripletFit.h:49
Belle2::QualityEstimationResults::pt
boost::optional< double > pt
transverse momentum estimate from the QE
Definition: QualityEstimatorBase.h:40
Belle2::SVD::SensorInfo
Specific implementation of SensorInfo for SVD Sensors which provides additional sensor specific infor...
Definition: SensorInfo.h:35
Belle2::QualityEstimatorTripletFit::m_materialBudgetFactor
double m_materialBudgetFactor
Triplet Fit hyper parameters.
Definition: QualityEstimatorTripletFit.h:62
Belle2::QualityEstimatorTripletFit::estimateQualityAndProperties
virtual QualityEstimationResults estimateQualityAndProperties(std::vector< SpacePoint const * > const &measurements) final
perform quality estimation and extract additional information from the fit
Definition: QualityEstimatorTripletFit.cc:170
Belle2::Const::speedOfLight
static const double speedOfLight
[cm/ns]
Definition: Const.h:568
Belle2::VxdID::baseType
unsigned short baseType
The base integer type for VxdID.
Definition: VxdID.h:46
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::QualityEstimatorBase::m_magneticFieldZ
double m_magneticFieldZ
Member storing the z component of the magnetic field.
Definition: QualityEstimatorBase.h:95
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::VXD::SensorInfoBase::SVD
@ SVD
SVD Sensor.
Definition: SensorInfoBase.h:45
Belle2::QualityEstimatorTripletFit::estimateQuality
virtual double estimateQuality(std::vector< SpacePoint const * > const &measurements) final
Calculating the quality estimate using a triplet fit.
Definition: QualityEstimatorTripletFit.cc:24
Belle2::QualityEstimatorTripletFit::m_sigmaR3DSquareds
std::vector< double > m_sigmaR3DSquareds
squared error of 3D radius
Definition: QualityEstimatorTripletFit.h:52
Belle2::QualityEstimatorTripletFit::m_R3Ds
std::vector< double > m_R3Ds
3D radius
Definition: QualityEstimatorTripletFit.h:51
Belle2::QualityEstimatorTripletFit::m_averageR3D
double m_averageR3D
average 3D radius
Definition: QualityEstimatorTripletFit.h:54
Belle2::QualityEstimationResults::curvatureSign
boost::optional< short > curvatureSign
direction of curvature as obtained by the QE
Definition: QualityEstimatorBase.h:39
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:79
Belle2::QualityEstimatorTripletFit::m_maxPt
double m_maxPt
Cut off value for 3D Radius calculated in Triplet Fit.
Definition: QualityEstimatorTripletFit.h:67
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