Belle II Software development
KLMTrackFitter Class Reference

track fitting procedure More...

#include <KLMTrackFitter.h>

Public Member Functions

 KLMTrackFitter ()
 Default constructor.
 
 ~KLMTrackFitter ()
 Destructor.
 
double fit (std::list< KLMHit2d * > &listTrackPoint)
 do fit and returns chi square of the fit.
 
double globalDistanceToHit (KLMHit2d *hit, double &error, double &sigma)
 Distance from track to a hit in the global system.
 
double fit1dTrack (std::list< KLMHit2d * > hitList, CLHEP::HepVector &eta, CLHEP::HepSymMatrix &error, int depDir, int indDir)
 do fit in the global system
 
CLHEP::HepVector getTrackParam ()
 Get track parameters in the global system. y = p0 + p1 * x; z = p2 + p3 * x.
 
CLHEP::HepSymMatrix getTrackParamErr ()
 Get invariance matrix of track parameters in the global system.
 
bool isValid ()
 Is fit valid.
 
bool isGood ()
 Is fit good.
 
float getChi2 ()
 Chi square of the fit.
 
int getNumHit ()
 number of the hits on this track
 
void inValidate ()
 Invalidate track.
 

Private Attributes

bool m_Valid
 Is fit valid.
 
bool m_Good
 Is fit good.
 
float m_Chi2
 Chi square of fit.
 
int m_NumHit
 the number of hits on this track
 
CLHEP::HepVector m_GlobalPar
 track params in global system
 
CLHEP::HepSymMatrix m_GlobalErr
 track params errors in global system
 
Belle2::KLM::KLMGeometryParm_GeoPar
 pointer to GeometryPar singleton
 

Detailed Description

track fitting procedure

Definition at line 30 of file KLMTrackFitter.h.

Constructor & Destructor Documentation

◆ KLMTrackFitter()

Default constructor.

Constructor.

Definition at line 42 of file KLMTrackFitter.cc.

42 :
43// Description: "For the use of standalone KLM tracking, this module is the fitter component.";
44 m_Valid(false),
45 m_Good(false),
46 m_Chi2(0.0),
47 m_NumHit(0),
48 m_GlobalPar(4, 0),
49 m_GlobalErr(4, 0),
50 m_GeoPar(nullptr)
51{
52}
float m_Chi2
Chi square of fit.
int m_NumHit
the number of hits on this track
Belle2::KLM::KLMGeometryPar * m_GeoPar
pointer to GeometryPar singleton
CLHEP::HepSymMatrix m_GlobalErr
track params errors in global system
CLHEP::HepVector m_GlobalPar
track params in global system
bool m_Valid
Is fit valid.
bool m_Good
Is fit good.

◆ ~KLMTrackFitter()

Destructor.

Definition at line 55 of file KLMTrackFitter.cc.

56{
57}

Member Function Documentation

◆ fit()

double fit ( std::list< KLMHit2d * > &  listTrackPoint)

do fit and returns chi square of the fit.

Definition at line 59 of file KLMTrackFitter.cc.

60{
61
62 // We can only do fit if there are at least two hits
63 if (listHitSector.size() < 2) {
64 m_Valid = false;
65 m_Good = false;
66 return (false);
67 }
68
69 HepVector eta(2, 0); // ( a, b ) in y = a + bx, used as inbetween
70 HepSymMatrix error(2, 0);
71 HepVector gloEta(2, 0); // ( a, b ) in y = a + bx in global system
72 HepSymMatrix gloError(2, 0);
73 m_Chi2 = 0;
74
75 // Create temporary vector and matrix... so size can be set.
76 HepVector globalPar(4, 0);
77 HepSymMatrix globalErr(4, 0);
78
79 m_Chi2 = fit1dTrack(listHitSector, eta, error, VY, VX);
80 globalPar.sub(1, eta);
81 globalErr.sub(1, error);
82
83 //m_Chi2 += fit1dTrack(listHitSector, eta, error, VY, VZ);
84 m_Chi2 += fit1dTrack(listHitSector, eta, error, VZ, VX);
85 globalPar.sub(3, eta);
86 globalErr.sub(3, error);
87
88
89 m_Chi2 /= 2.0;
90
91 m_GlobalPar = globalPar;
92 m_GlobalErr = globalErr;
93
94 m_Valid = true;
95 m_Good = false;
96 m_NumHit = listHitSector.size();
97 if (m_Chi2 <= 20.0) {
98 m_Good = true;
99 }
100
101 return (m_Chi2);
102
103}
double fit1dTrack(std::list< KLMHit2d * > hitList, CLHEP::HepVector &eta, CLHEP::HepSymMatrix &error, int depDir, int indDir)
do fit in the global system

◆ fit1dTrack()

double fit1dTrack ( std::list< KLMHit2d * >  hitList,
CLHEP::HepVector &  eta,
CLHEP::HepSymMatrix &  error,
int  depDir,
int  indDir 
)

do fit in the global system

do fit in global system, handle tracks that go thrugh multi-sectors

Definition at line 253 of file KLMTrackFitter.cc.

257{
258// Fit d = a + bi, where d is dependent direction and i is independent
259// in global system we assume y = a + b*x and z = c + dx different from local fit
260// Check KLMTrackFitter::fit to see what modes are being used
261 HepMatrix globalHitErr(3, 3, 0);
262
263 double indPos = 0;
264 double depPos = 0;
265
266 // Matrix based solution
267 int noPoints = hitList.size();
268
269 // Derivative of y = a + bx, with respect to a and b evaluated at x.
270 HepMatrix A(noPoints, 2, 0);
271
272 // Measured data points.
273 HepVector y(noPoints, 0);
274
275 // Inverse of covariance (error) matrix, also known as the weight matrix.
276 // In plain English: V_y_inverse_nn = 1 / (error of nth measurement)^2
277 HepDiagMatrix V_y_inverse(noPoints, 0);
278
279 // Error or correlation matrix for coefficients (2x2 matrices)
280 HepSymMatrix V_A, V_A_inverse;
281
282
283 int n = 0;
284 for (std::list< KLMHit2d* >::iterator iHit = hitList.begin(); iHit != hitList.end(); ++iHit) {
285
286 KLMHit2d* hit = *iHit;
287
288 CLHEP::Hep3Vector globalPos(hit->getPositionX(), hit->getPositionY(), hit->getPositionZ());
289
290
291 if (hit->getSubdetector() == KLMElementNumbers::c_BKLM) {
292 const bklm::GeometryPar* bklmGeo = m_GeoPar->BarrelInstance();
293 const bklm::Module* corMod = bklmGeo->findModule(hit->getSection(), hit->getSector(), hit->getLayer());
294
295 double hit_localPhiErr = corMod->getPhiStripWidth() / sqrt(12);
296 double hit_localZErr = corMod->getZStripWidth() / sqrt(12);
297
298 if (hit->inRPC()) {
299 //+++ scale localErr based on measured-in-Belle resolution
300 int nStrips = hit->getPhiStripMax() - hit->getPhiStripMin() + 1;
301 double dn = nStrips - 1.5;
302 double factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.60;//measured-in-Belle resolution
303 hit_localPhiErr = hit_localPhiErr * sqrt(factor);
304
305 nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
306 dn = nStrips - 1.5;
307 factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.55;//measured-in-Belle resolution
308 hit_localZErr = hit_localZErr * sqrt(factor);
309 }
310
311 Hep3Vector globalOrigin = corMod->getGlobalOrigin();
312 double sinphi = globalOrigin[1] / globalOrigin.mag();
313 double cosphi = globalOrigin[0] / globalOrigin.mag();
314
315 globalHitErr[0][0] = pow(hit_localPhiErr * sinphi, 2); // x
316 globalHitErr[0][1] = (hit_localPhiErr * sinphi) * (hit_localPhiErr * cosphi);
317 globalHitErr[0][2] = 0.;
318 globalHitErr[1][1] = pow(hit_localPhiErr * cosphi, 2); // y
319 globalHitErr[1][0] = (hit_localPhiErr * sinphi) * (hit_localPhiErr * cosphi);
320 globalHitErr[1][2] = 0.;
321 globalHitErr[2][2] = pow(hit_localZErr, 2); // z
322 globalHitErr[2][0] = 0.;
323 globalHitErr[2][1] = 0.;
324
325 switch (indDir) {
326
327 case VX:
328 indPos = globalPos.x();
329 break;
330
331 case VY:
332 indPos = globalPos.y();
333 break;
334
335 case VZ:
336 indPos = globalPos.z();
337 break;
338
339 default:
340 B2DEBUG(20, "error in bklm in-TrackFit: illegal direction");
341
342 }
343
344 switch (depDir) {
345
346 case VX:
347 depPos = globalPos.x();
348 break;
349
350 case VY:
351 depPos = globalPos.y();
352 break;
353
354 case VZ:
355 depPos = globalPos.z();
356 break;
357
358 default:
359 B2DEBUG(20, "error in bklm dep-TrackFit: illegal direction");
360
361 }
362 } //end of BKLM section
363
364 else if (hit->getSubdetector() == KLMElementNumbers::c_EKLM) {
365 const EKLM::GeometryData* eklmGeo = m_GeoPar->EndcapInstance();
366
367 //here get the resolustion of a hit, repeated several times, ugly. should we store this in KLMHit2d object ?
368 double hit_xErr = (eklmGeo->getStripGeometry()->getWidth()) * (Unit::cm / CLHEP::cm) *
369 (hit->getXStripMax() - hit->getXStripMin() + 1) / sqrt(12);
370 double hit_yErr = (eklmGeo->getStripGeometry()->getWidth()) * (Unit::cm / CLHEP::cm) *
371 (hit->getYStripMax() - hit->getYStripMin() + 1) / sqrt(12);
372 double hit_zErr = 0; //KLMHit2d is always centred on the boundary between the x/y planes with ~0 uncertainty
373
374 B2DEBUG(28, "KLMTrackFitter" << " Width: " << eklmGeo->getStripGeometry()->getWidth() << " Vec_x: " << hit_xErr * sqrt(
375 12) << " Vec_y: " << hit_yErr * sqrt(12));
376
377 globalHitErr[0][0] = pow(hit_xErr, 2); //x
378 globalHitErr[0][1] = 0.;
379 globalHitErr[0][2] = 0.;
380 globalHitErr[1][1] = pow(hit_yErr, 2); //y
381 globalHitErr[1][0] = 0.;
382 globalHitErr[1][2] = 0.;
383 globalHitErr[2][2] = pow(hit_zErr, 2); //z
384 globalHitErr[2][0] = 0.;
385 globalHitErr[2][1] = 0.;
386
387
388
389 switch (indDir) {
390
391 case VX:
392 indPos = globalPos.x();
393 break;
394
395 case VY:
396 indPos = globalPos.y();
397 break;
398
399 case VZ:
400 indPos = globalPos.z();
401 break;
402
403 default:
404 B2DEBUG(20, "error in EKLM in-TrackFit: illegal direction");
405
406 }
407
408 switch (depDir) {
409
410 case VX:
411 depPos = globalPos.x();
412 break;
413
414 case VY:
415 depPos = globalPos.y();
416 break;
417
418 case VZ:
419 depPos = globalPos.z();
420 break;
421
422 default:
423 B2DEBUG(20, "error in EKLM dep-TrackFit: illegal direction");
424
425 }
426 } //end of EKLM section
427
428 A[ n ][ 0 ] = 1;
429 A[ n ][ 1 ] = indPos;
430
431 y[ n ] = depPos;
432
433 double error_raw = globalHitErr[indDir][indDir] + globalHitErr[depDir][depDir];
434 double weight = error_raw;
435 if (weight > 0) {
436 V_y_inverse[ n ][ n ] = 1.0 / weight;
437 } else {
438 V_y_inverse[ n ][ n ] = DBL_MAX;
439 }
440 ++n;//n is the index of measured points/hits
441 } // end of loop
442
443
444 V_A_inverse = V_y_inverse.similarityT(A);
445
446 int ierr = 0;
447 V_A = V_A_inverse.inverse(ierr);
448
449 eta = V_A * A.T() * V_y_inverse * y;
450 error = V_A;
451
452// Calculate residuals:
453 HepMatrix residual = y - A * eta;
454 HepMatrix chisqr = residual.T() * V_y_inverse * residual;
455
456 return (chisqr.trace());
457
458}
double getWidth() const
Get width.
const StripGeometry * getStripGeometry() const
Get strip geometry data.
EKLM geometry data.
Definition: GeometryData.h:38
KLM 2d hit.
Definition: KLMHit2d.h:33
static const bklm::GeometryPar * BarrelInstance()
Return a pointer to the bklm::GeometryPar instance.
static const EKLM::GeometryData * EndcapInstance()
Return a pointer to the EKLM::GeometryData instance.
static const double cm
Standard units with the value = 1.
Definition: Unit.h:47
Provides BKLM geometry parameters for simulation, reconstruction etc (from Gearbox or DataBase)
Definition: GeometryPar.h:37
const Module * findModule(int section, int sector, int layer) const
Get the pointer to the definition of a module.
Definition: GeometryPar.cc:721
Define the geometry of a BKLM module Each sector [octant] contains Modules.
Definition: Module.h:76
double getPhiStripWidth() const
Get phi-strip width.
Definition: Module.h:137
const CLHEP::Hep3Vector getGlobalOrigin() const
Return the position (in global coordinates) of this module's sensitive-volume origin.
Definition: Module.h:285
double getZStripWidth() const
Get z-strip width.
Definition: Module.h:155
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28

◆ getChi2()

float getChi2 ( )
inline

Chi square of the fit.

Definition at line 82 of file KLMTrackFitter.h.

83 {
84 return m_Chi2;
85 }

◆ getNumHit()

int getNumHit ( )
inline

number of the hits on this track

Definition at line 88 of file KLMTrackFitter.h.

89 {
90 return m_NumHit;
91 }

◆ getTrackParam()

CLHEP::HepVector getTrackParam ( )
inline

Get track parameters in the global system. y = p0 + p1 * x; z = p2 + p3 * x.

Definition at line 57 of file KLMTrackFitter.h.

58 {
59 return m_GlobalPar;
60 }

◆ getTrackParamErr()

CLHEP::HepSymMatrix getTrackParamErr ( )
inline

Get invariance matrix of track parameters in the global system.

Definition at line 63 of file KLMTrackFitter.h.

64 {
65 return m_GlobalErr;
66 }

◆ globalDistanceToHit()

double globalDistanceToHit ( KLMHit2d hit,
double &  error,
double &  sigma 
)

Distance from track to a hit in the global system.

Distance from track to a hit calculated in the global system.

Definition at line 108 of file KLMTrackFitter.cc.

111{
112
113 if (!m_Valid) {
114 error = DBL_MAX;
115 sigma = DBL_MAX;
116 return DBL_MAX;
117 }
118
119 //in global fit, we have two functions y = a + b*x and y = c + d*z
120 double x_mea = hit->getPositionX();
121 double y_mea = hit->getPositionY();
122 double z_mea = hit->getPositionZ();
123
124 //there is some subdetector dependence so let's define this first.
125 double x_pre, y_pre, z_pre, dx, dy, dz;
126
127 // Error is composed of four parts: error due to tracking;
128 // and error in hit, y(x) or z(x).
129 HepMatrix errors(2, 2, 0); // Matrix for errors
130 HepMatrix A(2, 4, 0); // Matrix for derivatives
131
132
133 //defining quanitites for hit errors
134 double hit_xErr; double hit_yErr; double hit_zErr;
135 HepMatrix globalHitErr(3, 3, 0);
136
137 // Derivatives of y (z) = a + bx with respect to a and b.
138 A[ MY ][ AY ] = 1. ;
139 A[ MY ][ BY ] = x_mea;
140 A[ MZ ][ AZ ] = 1;
141 A[ MZ ][ BZ ] = x_mea;
142
143
144 //error from trackPar is inclueded, error from y_mea is not included
145 errors = A * m_GlobalErr * A.T();
146
147 if (hit->getSubdetector() == KLMElementNumbers::c_BKLM) {
148
149 // defining terms relating to distance
150 x_pre = x_mea; //by definition
151 y_pre = m_GlobalPar[ AY ] + m_GlobalPar[ BY ] * x_mea;
152 z_pre = m_GlobalPar[ AZ ] + m_GlobalPar[ BZ ] * x_mea;
153
154 const bklm::GeometryPar* bklmGeo = m_GeoPar->BarrelInstance();
155 //here get the resolustion of a hit, repeated several times, ugly. should we store this in KLMHit2d object ?
156 const bklm::Module* corMod = bklmGeo->findModule(hit->getSection(), hit->getSector(), hit->getLayer());
157 double hit_localPhiErr = corMod->getPhiStripWidth() / sqrt(12);
158 double hit_localZErr = corMod->getZStripWidth() / sqrt(12);
159
160 if (hit->inRPC()) {
161 //+++ scale localErr based on measured-in-Belle resolution
162 int nStrips = hit->getPhiStripMax() - hit->getPhiStripMin() + 1;
163 double dn = nStrips - 1.5;
164 double factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.60;//measured-in-Belle resolution
165 hit_localPhiErr = hit_localPhiErr * sqrt(factor);
166
167 nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
168 dn = nStrips - 1.5;
169 factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.55;//measured-in-Belle resolution
170 hit_localZErr = hit_localZErr * sqrt(factor);
171 }
172
173 Hep3Vector globalOrigin = corMod->getGlobalOrigin();
174 double sinphi = globalOrigin[1] / globalOrigin.mag();
175 double cosphi = globalOrigin[0] / globalOrigin.mag();
176
177 hit_xErr = hit_localPhiErr * sinphi;
178 hit_yErr = hit_localPhiErr * cosphi;
179 hit_zErr = hit_localZErr;
180
181
182 globalHitErr[0][0] = pow(hit_xErr, 2); //x
183 globalHitErr[0][1] = (hit_xErr) * (hit_yErr);
184 globalHitErr[0][2] = 0;
185 globalHitErr[1][1] = pow(hit_yErr, 2);;
186 globalHitErr[1][0] = (hit_xErr) * (hit_yErr);
187 globalHitErr[1][2] = 0;
188 globalHitErr[2][2] = pow(hit_zErr, 2);
189 globalHitErr[2][0] = 0;
190 globalHitErr[2][1] = 0;
191
192
193 } // end of BKLM portion
194
195 else if (hit->getSubdetector() == KLMElementNumbers::c_EKLM) {
196 // distance related (choose z variable)
197 // z coordinate should correspond to where track passes through EKLM-plane
198 z_pre = z_mea; //by design, due to EKLM geometry
199 x_pre = (z_mea - m_GlobalPar[ AZ ]) / m_GlobalPar[ BZ ];
200 y_pre = m_GlobalPar[ AY ] + m_GlobalPar[ BY ] * x_pre;
201
202 const EKLM::GeometryData* eklmGeo = m_GeoPar->EndcapInstance();
203
204 //EKLM GeometryData is needed for strip widths
205 //here get the resolustion of a hit, repeated several times, ugly. should we store this in KLMHit2d object ?
206 hit_xErr = (eklmGeo->getStripGeometry()->getWidth()) * (Unit::cm / CLHEP::cm) *
207 (hit->getXStripMax() - hit->getXStripMin() + 1) / sqrt(12);
208 hit_yErr = (eklmGeo->getStripGeometry()->getWidth()) * (Unit::cm / CLHEP::cm) *
209 (hit->getYStripMax() - hit->getYStripMin() + 1) / sqrt(12);
210 hit_zErr = 0.; //KLMHit2d is always centred on the boundary between the x/y planes with ~0 uncertainty
211
212 globalHitErr[0][0] = pow(hit_xErr, 2); //x
213 globalHitErr[0][1] = 0.;
214 globalHitErr[0][2] = 0.;
215 globalHitErr[1][1] = pow(hit_yErr, 2); //y
216 globalHitErr[1][0] = 0.;
217 globalHitErr[1][2] = 0.;
218 globalHitErr[2][2] = pow(hit_zErr, 2);; //z
219 globalHitErr[2][0] = 0.;
220 globalHitErr[2][1] = 0.;
221
222
223 } // end of EKLM section
224
225 else
226 B2FATAL("In KLMTrackFitter globalDistanceToHit, hit is neither from E/B-KLM.");
227
228 dx = x_pre - x_mea;
229 dy = y_pre - y_mea;
230 dz = z_pre - z_mea;
231 double distance = sqrt(dx * dx + dy * dy + dz * dz);
232
233 // now that we have globalHitErr, compute error
234 error = sqrt(errors[ MY ][ MY ] +
235 errors[ MZ ][ MZ ] +
236 //errors_b[0] + errors_b[1] + //error of prediced point due to error of inPos (y_mea)
237 globalHitErr[0][0] +
238 globalHitErr[1][1] + //y_mea error is correlated to the error of predicted point, but we didn't consider that part in errors
239 globalHitErr[2][2]); //we ignore y_mea error?
240
241
242 if (error != 0.0) {
243 sigma = distance / error;
244 } else {
245 sigma = DBL_MAX;
246 }
247
248 return (distance);
249}

◆ inValidate()

void inValidate ( )
inline

Invalidate track.

Definition at line 94 of file KLMTrackFitter.h.

95 {
96 m_Valid = false;
97 }

◆ isGood()

bool isGood ( )
inline

Is fit good.

Definition at line 76 of file KLMTrackFitter.h.

77 {
78 return m_Good;
79 }

◆ isValid()

bool isValid ( )
inline

Is fit valid.

Definition at line 70 of file KLMTrackFitter.h.

71 {
72 return m_Valid;
73 }

Member Data Documentation

◆ m_Chi2

float m_Chi2
private

Chi square of fit.

Definition at line 109 of file KLMTrackFitter.h.

◆ m_GeoPar

Belle2::KLM::KLMGeometryPar* m_GeoPar
private

pointer to GeometryPar singleton

Definition at line 121 of file KLMTrackFitter.h.

◆ m_GlobalErr

CLHEP::HepSymMatrix m_GlobalErr
private

track params errors in global system

Definition at line 118 of file KLMTrackFitter.h.

◆ m_GlobalPar

CLHEP::HepVector m_GlobalPar
private

track params in global system

Definition at line 115 of file KLMTrackFitter.h.

◆ m_Good

bool m_Good
private

Is fit good.

Definition at line 106 of file KLMTrackFitter.h.

◆ m_NumHit

int m_NumHit
private

the number of hits on this track

Definition at line 112 of file KLMTrackFitter.h.

◆ m_Valid

bool m_Valid
private

Is fit valid.

Definition at line 103 of file KLMTrackFitter.h.


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