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