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