Belle II Software  release-05-01-25
BKLMTrackFitter.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Yinghui GUAN *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 /* Own header. */
12 #include <klm/bklm/modules/bklmTracking/BKLMTrackFitter.h>
13 
14 /* KLM headers. */
15 #include <klm/bklm/geometry/GeometryPar.h>
16 #include <klm/bklm/geometry/Module.h>
17 
18 /* Belle 2 headers. */
19 #include <framework/logging/Logger.h>
20 
21 /* CLHEP headers. */
22 #include <CLHEP/Matrix/DiagMatrix.h>
23 #include <CLHEP/Matrix/Matrix.h>
24 
25 /* C++ headers. */
26 #include <cfloat>
27 
28 using namespace std;
29 using namespace CLHEP;
30 using namespace Belle2;
31 using namespace Belle2::bklm;
32 
34 enum { VX = 0, VY = 1, VZ = 2 };
35 
37 enum { AY = 0, BY = 1, AZ = 2, BZ = 3 };
38 
40 enum { MY = 0, MZ = 1 };
41 
43 BKLMTrackFitter::BKLMTrackFitter():
44  m_Valid(false),
45  m_Good(false),
46  m_Chi2(0.0),
47  m_NumHit(0),
48  m_globalFit(false),
49  m_SectorPar(4, 0),
50  m_SectorErr(4, 0),
51  m_GlobalPar(4, 0),
52  m_GlobalErr(4, 0),
53  m_GeoPar(nullptr)
54 {
55 }
56 
58 BKLMTrackFitter::~BKLMTrackFitter()
59 {
60 }
61 
63 double BKLMTrackFitter::fit(std::list<BKLMHit2d* >& listHitSector)
64 {
65 
66  // We can only do fit if there are at least two hits
67  if (listHitSector.size() < 2) {
68  m_Valid = false;
69  m_Good = false;
70  return (false);
71  }
72 
73  HepVector eta(2, 0); // ( a, b ) in y = a + bx
74  HepSymMatrix error(2, 0);
75  HepVector gloEta(2, 0); // ( a, b ) in y = a + bx in global system
76  HepSymMatrix gloError(2, 0);
77  m_Chi2 = 0;
78 
79  // Create temporary vector and matrix... so size can be set.
80  HepVector sectorPar(4, 0);
81  HepSymMatrix sectorErr(4, 0);
82  HepVector globalPar(4, 0);
83  HepSymMatrix globalErr(4, 0);
84 
85  if (m_globalFit) {
86  m_Chi2 = fit1dTrack(listHitSector, eta, error, VY, VX);
87  globalPar.sub(1, eta);
88  globalErr.sub(1, error);
89  } else {
90  m_Chi2 = fit1dSectorTrack(listHitSector, eta, error, VY, VX);
91  sectorPar.sub(1, eta);
92  sectorErr.sub(1, error);
93  }
94 
95  if (m_globalFit) {
96  m_Chi2 += fit1dTrack(listHitSector, eta, error, VY, VZ);
97  globalPar.sub(3, eta);
98  globalErr.sub(3, error);
99  } else {
100  m_Chi2 += fit1dSectorTrack(listHitSector, eta, error, VZ, VX);
101  sectorPar.sub(3, eta);
102  sectorErr.sub(3, error);
103  }
104 
105  if (!m_globalFit) {
106  //transfer to the global system, choose two abitrary points on track within in the sector jpionts on this track
107  const Belle2::bklm::Module* refMod = m_GeoPar->findModule((*listHitSector.begin())->getSection(),
108  (*listHitSector.begin())->getSector(), 1);
109 
110  Hep3Vector p1(0, 0, 0); Hep3Vector p2(0, 0, 0);
111  double x1 = 5; // sector localX
112  double x2 = 10;
113  double y1 = sectorPar[0] + sectorPar[1] * x1;
114  double y2 = sectorPar[0] + sectorPar[1] * x2;
115  double z1 = sectorPar[2] + sectorPar[3] * x1;
116  double z2 = sectorPar[2] + sectorPar[3] * x2;
117  p1.setX(x1); p1.setY(y1); p1.setZ(z1);
118  p2.setX(x2); p2.setY(y2); p2.setZ(z2);
119  Hep3Vector gl1 = refMod->localToGlobal(p1);
120  Hep3Vector gl2 = refMod->localToGlobal(p2);
121 
122  //transfer the trackParameters to global system y = p0 + p1 * x and z= p2 + p3*x
123  if (gl2[0] != gl1[0]) {
124  globalPar[1] = (gl2[1] - gl1[1]) / (gl2[0] - gl1[0]);
125  globalPar[0] = gl1[1] - globalPar[1] * gl1[0];
126  globalPar[3] = (gl2[2] - gl1[2]) / (gl2[0] - gl1[0]);
127  globalPar[2] = gl1[2] - globalPar[3] * gl1[0];
128  globalErr = sectorErr;
129  } else {
130  globalPar[1] = DBL_MAX;
131  globalPar[0] = DBL_MAX;
132  globalPar[3] = DBL_MAX;
133  globalPar[2] = DBL_MAX;
134  globalErr = sectorErr; //?
135  }
136  } else { //transfer to the local system, can not do. One global track might go through two different sectors.
137  }
138 
139  m_Chi2 /= 2.0;
140 
141  m_SectorPar = sectorPar;
142  m_SectorErr = sectorErr;
143  m_GlobalPar = globalPar;
144  m_GlobalErr = globalErr;
145 
146  m_Valid = true;
147  m_Good = false;
148  m_NumHit = listHitSector.size();
149  if (m_Chi2 <= 20.0) {
150  m_Good = true;
151  }
152 
153  return (m_Chi2);
154 
155 }
156 
158 double BKLMTrackFitter::distanceToHit(BKLMHit2d* hit,
159  double& error,
160  double& sigma)
161 {
162 
163  double x, y, z, dy, dz;
164 
165  if (!m_Valid) {
166  error = DBL_MAX;
167  sigma = DBL_MAX;
168  return DBL_MAX;
169  }
170 
171  m_GeoPar = GeometryPar::instance();
172  const Belle2::bklm::Module* refMod = m_GeoPar->findModule(hit->getSection(), hit->getSector(), 1);
173  const Belle2::bklm::Module* corMod = m_GeoPar->findModule(hit->getSection(), hit->getSector(), hit->getLayer());
174 
175  CLHEP::Hep3Vector globalPos(hit->getGlobalPosition()[0], hit->getGlobalPosition()[1], hit->getGlobalPosition()[2]);
176 
177  //+++ local coordinates of the hit
178  CLHEP::Hep3Vector local = refMod->globalToLocal(globalPos);
179 
180  x = local[0] ;
181 
182  y = m_SectorPar[ AY ] + x * m_SectorPar[ BY ];
183  z = m_SectorPar[ AZ ] + x * m_SectorPar[ BZ ];
184 
185  dy = y - local[1];
186  dz = z - local[2];
187 
188  double distance = sqrt(dy * dy + dz * dz);
189 
190  // Error is composed of four parts: error due to tracking, y and z;
191  // and error in hit, y and z. We know the latter two, got to find
192  // the first two. We could calculate this from simple equations or
193  // using matrices. I choose the later because it is extendable.
194  HepMatrix errors(2, 2, 0); // Matrix for errors
195  HepMatrix A(2, 4, 0); // Matrix for derivatives
196 
197  // Derivatives of y (z) = a + bx with respect to a and b.
198  A[ MY ][ AY ] = 1.0;
199  A[ MY ][ BY ] = x;
200  A[ MZ ][ AZ ] = 1.0;
201  A[ MZ ][ BZ ] = x;
202 
203  errors = A * m_SectorErr * A.T();
204 
205  double hit_localPhiErr = corMod->getPhiStripWidth() / sqrt(12);
206  double hit_localZErr = corMod->getZStripWidth() / sqrt(12);
207 
208  if (hit->inRPC()) {
209  //+++ scale localErr based on measured-in-Belle resolution
210  int nStrips = hit->getPhiStripMax() - hit->getPhiStripMin() + 1;
211  double dn = nStrips - 1.5;
212  double factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.60;//measured-in-Belle resolution
213  hit_localPhiErr = hit_localPhiErr * sqrt(factor);
214 
215  nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
216  dn = nStrips - 1.5;
217  factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.55;//measured-in-Belle resolution
218  hit_localZErr = hit_localZErr * sqrt(factor);
219  }
220 
221  error = sqrt(errors[ MY ][ MY ] +
222  errors[ MZ ][ MZ ] +
223  pow(hit_localPhiErr, 2) +
224  pow(hit_localZErr, 2));
225 
226  if (error != 0.0) {
227  sigma = distance / error;
228  } else {
229  sigma = DBL_MAX;
230  }
231 
232  return (distance);
233 }
234 
236 double BKLMTrackFitter::globalDistanceToHit(BKLMHit2d* hit,
237  double& error,
238  double& sigma)
239 {
240 
241  if (!m_Valid) {
242  error = DBL_MAX;
243  sigma = DBL_MAX;
244  return DBL_MAX;
245  }
246 
247  //in global fit, we have two functions y = a + b*x and y = c + d*z
248  double x_mea = hit->getGlobalPosition()[0];
249  double y_mea = hit->getGlobalPosition()[1];
250  double z_mea = hit->getGlobalPosition()[2];
251 
252  double x_pre = (y_mea - m_GlobalPar[ AY ]) / m_GlobalPar[ BY ]; //y_mea has uncertainties actually
253  double z_pre = (y_mea - m_GlobalPar[ AZ ]) / m_GlobalPar[ BZ ];
254 
255  double dx = x_pre - x_mea;
256  double dz = z_pre - z_mea;
257 
258  double distance = sqrt(dx * dx + dz * dz);
259 
260  // Error is composed of four parts: error due to tracking;
261  // and error in hit, y, x or y, z.
262  HepMatrix errors(2, 2, 0); // Matrix for errors
263  HepMatrix A(2, 4, 0); // Matrix for derivatives
264 
265  // Derivatives of x (z) = y/b - a/b with respect to a and b.
266  A[ MY ][ AY ] = -1. / m_GlobalPar[BY];
267  A[ MY ][ BY ] = -1 * (y_mea - m_GlobalPar[ AY ]) / (m_GlobalPar[ BY ] * m_GlobalPar[ BY ]);
268  A[ MZ ][ AZ ] = -1. / m_GlobalPar[ BZ ];
269  A[ MZ ][ BZ ] = -1 * (y_mea - m_GlobalPar[ AZ ]) / (m_GlobalPar[ BZ ] * m_GlobalPar[ BZ ]);
270 
271  //error from trackPar is inclueded, error from y_mea is not included
272  errors = A * m_GlobalErr * A.T();
273 
274  //here get the resolustion of a hit, repeated several times, ugly. should we store this in BKLMHit2d object ?
275  const Belle2::bklm::Module* corMod = m_GeoPar->findModule(hit->getSection(), hit->getSector(), hit->getLayer());
276  double hit_localPhiErr = corMod->getPhiStripWidth() / sqrt(12);
277  double hit_localZErr = corMod->getZStripWidth() / sqrt(12);
278 
279  if (hit->inRPC()) {
280  //+++ scale localErr based on measured-in-Belle resolution
281  int nStrips = hit->getPhiStripMax() - hit->getPhiStripMin() + 1;
282  double dn = nStrips - 1.5;
283  double factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.60;//measured-in-Belle resolution
284  hit_localPhiErr = hit_localPhiErr * sqrt(factor);
285 
286  nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
287  dn = nStrips - 1.5;
288  factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.55;//measured-in-Belle resolution
289  hit_localZErr = hit_localZErr * sqrt(factor);
290  }
291 
292  Hep3Vector globalOrigin = corMod->getGlobalOrigin();
293  double sinphi = globalOrigin[1] / globalOrigin.mag();
294  double cosphi = globalOrigin[0] / globalOrigin.mag();
295 
296  HepMatrix globalHitErr(3, 3, 0);
297  globalHitErr[0][0] = pow(hit_localPhiErr * sinphi, 2); //x
298  globalHitErr[0][1] = (hit_localPhiErr * sinphi) * (hit_localPhiErr * cosphi);
299  globalHitErr[0][2] = 0;
300  globalHitErr[1][1] = pow(hit_localPhiErr * cosphi, 2);;
301  globalHitErr[1][0] = (hit_localPhiErr * sinphi) * (hit_localPhiErr * cosphi);
302  globalHitErr[1][2] = 0;
303  globalHitErr[2][2] = pow(hit_localZErr, 2);
304  globalHitErr[2][0] = 0;
305  globalHitErr[2][1] = 0;
306 
307  //HepMatrix B(2, 2, 0); // Matrix for derivatives
308 
309  // Derivatives of x (z) = y/b - a/b with respect to y.
310  //B[ MY ][ AY ] = 1./m_GlobalPar[BY];
311  //B[ MZ ][ BY ] = 1./m_GlobalPar[BZ] ;
312  //double errors_b[0] = B[ MY ][ AY ]*globalHitErr[1][1];
313  //double errors_b[1] = B[ MZ ][ BY ]*globalHitErr[1][1];
314 
315  HepMatrix error_mea(2, 2, 0);
316  error = sqrt(errors[ MY ][ MY ] +
317  errors[ MZ ][ MZ ] +
318  //errors_b[0] + errors_b[1] + //error of prediced point due to error of inPos (y_mea)
319  globalHitErr[0][0] +
320  globalHitErr[1][1] + //y_mea error is correlated to the error of predicted point, but we didn't consider that part in errors
321  globalHitErr[2][2]); //we ignore y_mea error?
322 
323  if (error != 0.0) {
324  sigma = distance / error;
325  } else {
326  sigma = DBL_MAX;
327  }
328 
329  return (distance);
330 }
331 
332 
334 double BKLMTrackFitter::fit1dSectorTrack(std::list< BKLMHit2d* > hitList,
335  HepVector& eta,
336  HepSymMatrix& error,
337  int depDir, int indDir)
338 {
339 
340 // Fit d = a + bi, where d is dependent direction and i is independent
341 
342  std::list< BKLMHit2d* >::iterator iHit;
343 
344  Hep3Vector localHitPos;
345  HepMatrix localHitErr(3, 3, 0);
346 
347  double indPos = 0;
348  double depPos = 0;
349 
350  // Matrix based solution
351 
352  int noPoints = hitList.size();
353 
354  // Derivative of y = a + bx, with respect to a and b evaluated at x.
355  HepMatrix A(noPoints, 2, 0);
356 
357  // Measured data points.
358  HepVector y(noPoints, 0);
359 
360  // Inverse of covariance (error) matrix, also known as the weight matrix.
361  // In plain English: V_y_inverse_nn = 1 / (error of nth measurement)^2
362  HepDiagMatrix V_y_inverse(noPoints, 0);
363  //HepSymMatrix V_y_inverse(noPoints, 0);
364 
365  // Error or correlation matrix for coefficients (2x2 matrices)
366  HepSymMatrix V_A, V_A_inverse;
367 
368  iHit = hitList.begin();
369  BKLMHit2d* hit = *iHit;
370  int section = hit->getSection();
371  int sector = hit->getSector();
372 
373  m_GeoPar = GeometryPar::instance();
374  const Belle2::bklm::Module* refMod = m_GeoPar->findModule((*hitList.begin())->getSection(), (*hitList.begin())->getSector(), 1);
375 
376  int n = 0;
377  for (iHit = hitList.begin(); iHit != hitList.end(); ++iHit) {
378 
379  hit = *iHit;
380  if (hit->getSection() != section || hit->getSector() != sector) {
381  continue;
382  }
383 
384  // m_GeoPar = GeometryPar::instance();
385  //const Belle2::bklm::Module* refMod = m_GeoPar->findModule(hit->getSection(), hit->getSector(), 1);
386  const Belle2::bklm::Module* corMod = m_GeoPar->findModule(hit->getSection(), hit->getSector(), hit->getLayer());
387 
388  CLHEP::Hep3Vector globalPos;
389  globalPos[0] = hit->getGlobalPosition()[0];
390  globalPos[1] = hit->getGlobalPosition()[1];
391  globalPos[2] = hit->getGlobalPosition()[2];
392 
393  localHitPos = refMod->globalToLocal(globalPos);
394  double hit_localPhiErr = corMod->getPhiStripWidth() / sqrt(12);
395  double hit_localZErr = corMod->getZStripWidth() / sqrt(12);
396 
397  if (hit->inRPC()) {
398  //+++ scale localErr based on measured-in-Belle resolution
399  int nStrips = hit->getPhiStripMax() - hit->getPhiStripMin() + 1;
400  double dn = nStrips - 1.5;
401  double factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.60;//measured-in-Belle resolution
402  hit_localPhiErr = hit_localPhiErr * sqrt(factor);
403 
404  nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
405  dn = nStrips - 1.5;
406  factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.55;//measured-in-Belle resolution
407  hit_localZErr = hit_localZErr * sqrt(factor);
408  }
409 
410  localHitErr[0][0] = 0.0; //x
411  localHitErr[0][1] = 0;
412  localHitErr[0][2] = 0;
413  localHitErr[1][1] = hit_localPhiErr;
414  localHitErr[1][0] = 0;
415  localHitErr[1][2] = 0;
416  localHitErr[2][2] = hit_localZErr;
417  localHitErr[2][0] = 0;
418  localHitErr[2][1] = 0;
419 
420  switch (indDir) {
421 
422  case VX:
423  indPos = localHitPos.x();
424  break;
425 
426  case VY:
427  indPos = localHitPos.y();
428  break;
429 
430  case VZ:
431  indPos = localHitPos.z();
432  break;
433 
434  default:
435  B2DEBUG(20, "error in klm_trackSectorFit: illegal direction");
436 
437  }
438 
439  switch (depDir) {
440 
441  case VX:
442  depPos = localHitPos.x();
443  break;
444 
445  case VY:
446  depPos = localHitPos.y();
447  break;
448 
449  case VZ:
450  depPos = localHitPos.z();
451  break;
452 
453  default:
454  B2DEBUG(20, "error in klm_trackSectorFit: illegal direction");
455 
456  }
457 
458 
459  A[ n ][ 0 ] = 1;
460  A[ n ][ 1 ] = indPos;
461 
462  y[ n ] = depPos;
463 
464  if (localHitErr[ depDir ][ depDir ] > 0.0) {
465  V_y_inverse[ n ][ n ] = 1.0 / localHitErr[ depDir ][ depDir ];
466  } else {
467  V_y_inverse[ n ][ n ] = DBL_MAX;
468  }
469  ++n;//n is the index of measured points/hits
470  }
471 
472  V_A_inverse = V_y_inverse.similarityT(A);
473 
474  int ierr = 0;
475  V_A = V_A_inverse.inverse(ierr);
476 
477  eta = V_A * A.T() * V_y_inverse * y;
478  error = V_A;
479 
480 // Calculate residuals:
481  HepMatrix residual = y - A * eta;
482  HepMatrix chisqr = residual.T() * V_y_inverse * residual;
483 
484  return (chisqr.trace());
485 
486 }
487 
489 double BKLMTrackFitter::fit1dTrack(std::list< BKLMHit2d* > hitList,
490  HepVector& eta,
491  HepSymMatrix& error,
492  int depDir, int indDir)
493 {
494 // Fit d = a + bi, where d is dependent direction and i is independent
495 // in global system we assume y = a + b*x and y = c + d*z different from local fit
496 
497  HepMatrix globalHitErr(3, 3, 0);
498 
499  double indPos = 0;
500  double depPos = 0;
501 
502  // Matrix based solution
503  int noPoints = hitList.size();
504 
505  // Derivative of y = a + bx, with respect to a and b evaluated at x.
506  HepMatrix A(noPoints, 2, 0);
507 
508  // Measured data points.
509  HepVector y(noPoints, 0);
510 
511  // Inverse of covariance (error) matrix, also known as the weight matrix.
512  // In plain English: V_y_inverse_nn = 1 / (error of nth measurement)^2
513  HepDiagMatrix V_y_inverse(noPoints, 0);
514 
515  // Error or correlation matrix for coefficients (2x2 matrices)
516  HepSymMatrix V_A, V_A_inverse;
517 
518  m_GeoPar = GeometryPar::instance();
519  const Belle2::bklm::Module* corMod;
520 
521  int n = 0;
522  for (std::list< BKLMHit2d* >::iterator iHit = hitList.begin(); iHit != hitList.end(); ++iHit) {
523 
524  BKLMHit2d* hit = *iHit;
525 
526  // m_GeoPar = GeometryPar::instance();
527  //const Belle2::bklm::Module* refMod = m_GeoPar->findModule(hit->getSection(), hit->getSector(), 1);
528  corMod = m_GeoPar->findModule(hit->getSection(), hit->getSector(), hit->getLayer());
529 
530  CLHEP::Hep3Vector globalPos;
531  globalPos[0] = hit->getGlobalPosition()[0];
532  globalPos[1] = hit->getGlobalPosition()[1];
533  globalPos[2] = hit->getGlobalPosition()[2];
534 
535  //localHitPos = refMod->globalToLocal(globalPos);
536  double hit_localPhiErr = corMod->getPhiStripWidth() / sqrt(12);
537  double hit_localZErr = corMod->getZStripWidth() / sqrt(12);
538 
539  if (hit->inRPC()) {
540  //+++ scale localErr based on measured-in-Belle resolution
541  int nStrips = hit->getPhiStripMax() - hit->getPhiStripMin() + 1;
542  double dn = nStrips - 1.5;
543  double factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.60;//measured-in-Belle resolution
544  hit_localPhiErr = hit_localPhiErr * sqrt(factor);
545 
546  nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
547  dn = nStrips - 1.5;
548  factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.55;//measured-in-Belle resolution
549  hit_localZErr = hit_localZErr * sqrt(factor);
550  }
551 
552  Hep3Vector globalOrigin = corMod->getGlobalOrigin();
553  double sinphi = globalOrigin[1] / globalOrigin.mag();
554  double cosphi = globalOrigin[0] / globalOrigin.mag();
555 
556  globalHitErr[0][0] = pow(hit_localPhiErr * sinphi, 2); //x
557  globalHitErr[0][1] = (hit_localPhiErr * sinphi) * (hit_localPhiErr * cosphi);
558  globalHitErr[0][2] = 0;
559  globalHitErr[1][1] = pow(hit_localPhiErr * cosphi, 2);;
560  globalHitErr[1][0] = (hit_localPhiErr * sinphi) * (hit_localPhiErr * cosphi);
561  globalHitErr[1][2] = 0;
562  globalHitErr[2][2] = pow(hit_localZErr, 2);;
563  globalHitErr[2][0] = 0;
564  globalHitErr[2][1] = 0;
565 
566  switch (indDir) {
567 
568  case VX:
569  indPos = globalPos.x();
570  break;
571 
572  case VY:
573  indPos = globalPos.y();
574  break;
575 
576  case VZ:
577  indPos = globalPos.z();
578  break;
579 
580  default:
581  B2DEBUG(20, "error in bklm trackFit: illegal direction");
582 
583  }
584 
585  switch (depDir) {
586 
587  case VX:
588  depPos = globalPos.x();
589  break;
590 
591  case VY:
592  depPos = globalPos.y();
593  break;
594 
595  case VZ:
596  depPos = globalPos.z();
597  break;
598 
599  default:
600  B2DEBUG(20, "error in bklm trackFit: illegal direction");
601 
602  }
603 
604 
605  A[ n ][ 0 ] = 1;
606  A[ n ][ 1 ] = indPos;
607 
608  y[ n ] = depPos;
609 
610  double error_raw = globalHitErr[indDir][indDir] + globalHitErr[depDir][depDir];
611  //double correlate_ = 2.0*globalHitErr[indDir][depDir]; //?
612  //double weight= error_raw - correlate_;
613  double weight = error_raw;
614  if (weight > 0) {
615  V_y_inverse[ n ][ n ] = 1.0 / weight;
616  } else {
617  V_y_inverse[ n ][ n ] = DBL_MAX;
618  }
619  ++n;//n is the index of measured points/hits
620  }
621 
622  //HepMatrix AT = A.T();
623  //HepMatrix tmp = AT*V_y_inverse;
624  //V_A_inverse = tmp*A;
625  V_A_inverse = V_y_inverse.similarityT(A);
626 
627  int ierr = 0;
628  V_A = V_A_inverse.inverse(ierr);
629 
630  eta = V_A * A.T() * V_y_inverse * y;
631  error = V_A;
632 
633 // Calculate residuals:
634  HepMatrix residual = y - A * eta;
635  HepMatrix chisqr = residual.T() * V_y_inverse * residual;
636 
637  return (chisqr.trace());
638 
639 }
Belle2::BKLMTrackFitter::distanceToHit
double distanceToHit(BKLMHit2d *hit, double &error, double &sigma)
Distance from track to a hit in the plane of the module.
Definition: BKLMTrackFitter.cc:158
Belle2::BKLMTrackFitter::m_Chi2
float m_Chi2
Chi square of fit.
Definition: BKLMTrackFitter.h:143
Belle2::bklm::Module::getPhiStripWidth
double getPhiStripWidth() const
Get phi-strip width.
Definition: Module.h:147
Belle2::BKLMTrackFitter::m_SectorErr
CLHEP::HepSymMatrix m_SectorErr
track params errors in the sector local system
Definition: BKLMTrackFitter.h:155
Belle2::BKLMTrackFitter::m_SectorPar
CLHEP::HepVector m_SectorPar
track params in the sector local system
Definition: BKLMTrackFitter.h:152
Belle2::bklm::Module
Define the geometry of a BKLM module Each sector [octant] contains Modules.
Definition: Module.h:86
Belle2::bklm::Module::getGlobalOrigin
const CLHEP::Hep3Vector getGlobalOrigin() const
Return the position (in global coordinates) of this module's sensitive-volume origin.
Definition: Module.h:287
Belle2::BKLMTrackFitter::m_GlobalPar
CLHEP::HepVector m_GlobalPar
track params in global system
Definition: BKLMTrackFitter.h:158
Belle2::BKLMTrackFitter::m_Good
bool m_Good
Is fit good.
Definition: BKLMTrackFitter.h:140
Belle2::BKLMTrackFitter::fit1dTrack
double fit1dTrack(std::list< BKLMHit2d * > hitList, CLHEP::HepVector &eta, CLHEP::HepSymMatrix &error, int depDir, int indDir)
do fit in the global system
Definition: BKLMTrackFitter.cc:489
Belle2::bklm::Module::globalToLocal
const CLHEP::Hep3Vector globalToLocal(const CLHEP::Hep3Vector &v, bool reco=false) const
Transform space-point within this module from global to local coordinates.
Definition: Module.cc:320
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::BKLMTrackFitter::m_NumHit
int m_NumHit
the number of hits on this track
Definition: BKLMTrackFitter.h:146
Belle2::BKLMTrackFitter::m_globalFit
bool m_globalFit
do fit in the local system or global system false: local sys; true: global sys.
Definition: BKLMTrackFitter.h:149
Belle2::BKLMTrackFitter::m_GeoPar
bklm::GeometryPar * m_GeoPar
pointer to GeometryPar singleton
Definition: BKLMTrackFitter.h:164
Belle2::bklm::Module::localToGlobal
const CLHEP::Hep3Vector localToGlobal(const CLHEP::Hep3Vector &v, bool reco=false) const
Transform space-point within this module from local to global coordinates.
Definition: Module.cc:306
Belle2::BKLMTrackFitter::~BKLMTrackFitter
~BKLMTrackFitter()
Destructor.
Definition: BKLMTrackFitter.cc:58
Belle2::BKLMTrackFitter::m_GlobalErr
CLHEP::HepSymMatrix m_GlobalErr
track params errors in global system
Definition: BKLMTrackFitter.h:161
Belle2::BKLMTrackFitter::globalDistanceToHit
double globalDistanceToHit(BKLMHit2d *hit, double &error, double &sigma)
Distance from track to a hit in the global system.
Definition: BKLMTrackFitter.cc:236
Belle2::BKLMTrackFitter::fit1dSectorTrack
double fit1dSectorTrack(std::list< BKLMHit2d * > hitList, CLHEP::HepVector &eta, CLHEP::HepSymMatrix &error, int depDir, int indDir)
do fit in the y-x plane or z-x plane
Definition: BKLMTrackFitter.cc:334
Belle2::bklm::Module::getZStripWidth
double getZStripWidth() const
Get z-strip width.
Definition: Module.h:165
Belle2::BKLMTrackFitter::fit
double fit(std::list< BKLMHit2d * > &listTrackPoint)
do fit and returns chi square of the fit.
Definition: BKLMTrackFitter.cc:63
Belle2::bklm::GeometryPar::findModule
const Module * findModule(int section, int sector, int layer) const
Get the pointer to the definition of a module.
Definition: GeometryPar.cc:716
Belle2::BKLMTrackFitter::m_Valid
bool m_Valid
Is fit valid.
Definition: BKLMTrackFitter.h:137
Belle2::BKLMHit2d
Store one BKLM strip hit as a ROOT object.
Definition: BKLMHit2d.h:42