Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GblPoint.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GblPoint.cc
1 /*
2  * GblPoint.cpp
3  *
4  * Created on: Aug 18, 2011
5  * Author: kleinwrt
6  */
7 
30 #include "GblPoint.h"
31 
33 namespace gbl {
34 
36 
40 GblPoint::GblPoint(const TMatrixD &aJacobian) :
41  theLabel(0), theOffset(0), measDim(0), transFlag(false), measTransformation(), scatFlag(
42  false), localDerivatives(), globalLabels(), globalDerivatives() {
43 
44  for (unsigned int i = 0; i < 5; ++i) {
45  for (unsigned int j = 0; j < 5; ++j) {
46  p2pJacobian(i, j) = aJacobian(i, j);
47  }
48  }
49 }
50 
51 GblPoint::GblPoint(const SMatrix55 &aJacobian) :
52  theLabel(0), theOffset(0), p2pJacobian(aJacobian), measDim(0), transFlag(
53  false), measTransformation(), scatFlag(false), localDerivatives(), globalLabels(), globalDerivatives() {
54 
55 }
56 
58 }
59 
61 
69 void GblPoint::addMeasurement(const TMatrixD &aProjection,
70  const TVectorD &aResiduals, const TVectorD &aPrecision,
71  double minPrecision) {
72  measDim = aResiduals.GetNrows();
73  unsigned int iOff = 5 - measDim;
74  for (unsigned int i = 0; i < measDim; ++i) {
75  measResiduals(iOff + i) = aResiduals[i];
76  measPrecision(iOff + i) = (
77  aPrecision[i] >= minPrecision ? aPrecision[i] : 0.);
78  for (unsigned int j = 0; j < measDim; ++j) {
79  measProjection(iOff + i, iOff + j) = aProjection(i, j);
80  }
81  }
82 }
83 
85 
94 void GblPoint::addMeasurement(const TMatrixD &aProjection,
95  const TVectorD &aResiduals, const TMatrixDSym &aPrecision,
96  double minPrecision) {
97  measDim = aResiduals.GetNrows();
98  TMatrixDSymEigen measEigen(aPrecision);
100  measTransformation = measEigen.GetEigenVectors();
101  measTransformation.T();
102  transFlag = true;
103  TVectorD transResiduals = measTransformation * aResiduals;
104  TVectorD transPrecision = measEigen.GetEigenValues();
105  TMatrixD transProjection = measTransformation * aProjection;
106  unsigned int iOff = 5 - measDim;
107  for (unsigned int i = 0; i < measDim; ++i) {
108  measResiduals(iOff + i) = transResiduals[i];
109  measPrecision(iOff + i) = (
110  transPrecision[i] >= minPrecision ? transPrecision[i] : 0.);
111  for (unsigned int j = 0; j < measDim; ++j) {
112  measProjection(iOff + i, iOff + j) = transProjection(i, j);
113  }
114  }
115 }
116 
118 
125 void GblPoint::addMeasurement(const TVectorD &aResiduals,
126  const TVectorD &aPrecision, double minPrecision) {
127  measDim = aResiduals.GetNrows();
128  unsigned int iOff = 5 - measDim;
129  for (unsigned int i = 0; i < measDim; ++i) {
130  measResiduals(iOff + i) = aResiduals[i];
131  measPrecision(iOff + i) = (
132  aPrecision[i] >= minPrecision ? aPrecision[i] : 0.);
133  }
134  measProjection = ROOT::Math::SMatrixIdentity();
135 }
136 
138 
146 void GblPoint::addMeasurement(const TVectorD &aResiduals,
147  const TMatrixDSym &aPrecision, double minPrecision) {
148  measDim = aResiduals.GetNrows();
149  TMatrixDSymEigen measEigen(aPrecision);
151  measTransformation = measEigen.GetEigenVectors();
152  measTransformation.T();
153  transFlag = true;
154  TVectorD transResiduals = measTransformation * aResiduals;
155  TVectorD transPrecision = measEigen.GetEigenValues();
156  unsigned int iOff = 5 - measDim;
157  for (unsigned int i = 0; i < measDim; ++i) {
158  measResiduals(iOff + i) = transResiduals[i];
159  measPrecision(iOff + i) = (
160  transPrecision[i] >= minPrecision ? transPrecision[i] : 0.);
161  for (unsigned int j = 0; j < measDim; ++j) {
162  measProjection(iOff + i, iOff + j) = measTransformation(i, j);
163  }
164  }
165 }
166 
168 
172 unsigned int GblPoint::hasMeasurement() const {
173  return measDim;
174 }
175 
177 
182 void GblPoint::getMeasurement(SMatrix55 &aProjection, SVector5 &aResiduals,
183  SVector5 &aPrecision) const {
184  aProjection = measProjection;
185  aResiduals = measResiduals;
186  aPrecision = measPrecision;
187 }
188 
190 
193 void GblPoint::getMeasTransformation(TMatrixD &aTransformation) const {
194  aTransformation.ResizeTo(measDim, measDim);
195  if (transFlag) {
196  aTransformation = measTransformation;
197  } else {
198  aTransformation.UnitMatrix();
199  }
200 }
201 
203 
210 void GblPoint::addScatterer(const TVectorD &aResiduals,
211  const TVectorD &aPrecision) {
212  scatFlag = true;
213  scatResiduals(0) = aResiduals[0];
214  scatResiduals(1) = aResiduals[1];
215  scatPrecision(0) = aPrecision[0];
216  scatPrecision(1) = aPrecision[1];
217  scatTransformation = ROOT::Math::SMatrixIdentity();
218 }
219 
221 
236 void GblPoint::addScatterer(const TVectorD &aResiduals,
237  const TMatrixDSym &aPrecision) {
238  scatFlag = true;
239  TMatrixDSymEigen scatEigen(aPrecision);
240  TMatrixD aTransformation = scatEigen.GetEigenVectors();
241  aTransformation.T();
242  TVectorD transResiduals = aTransformation * aResiduals;
243  TVectorD transPrecision = scatEigen.GetEigenValues();
244  for (unsigned int i = 0; i < 2; ++i) {
245  scatResiduals(i) = transResiduals[i];
246  scatPrecision(i) = transPrecision[i];
247  for (unsigned int j = 0; j < 2; ++j) {
248  scatTransformation(i, j) = aTransformation(i, j);
249  }
250  }
251 }
252 
255  return scatFlag;
256 }
257 
259 
264 void GblPoint::getScatterer(SMatrix22 &aTransformation, SVector2 &aResiduals,
265  SVector2 &aPrecision) const {
266  aTransformation = scatTransformation;
267  aResiduals = scatResiduals;
268  aPrecision = scatPrecision;
269 }
270 
272 
275 void GblPoint::getScatTransformation(TMatrixD &aTransformation) const {
276  aTransformation.ResizeTo(2, 2);
277  if (scatFlag) {
278  for (unsigned int i = 0; i < 2; ++i) {
279  for (unsigned int j = 0; j < 2; ++j) {
280  aTransformation(i, j) = scatTransformation(i, j);
281  }
282  }
283  } else {
284  aTransformation.UnitMatrix();
285  }
286 }
287 
289 
293 void GblPoint::addLocals(const TMatrixD &aDerivatives) {
294  if (measDim) {
295  localDerivatives.ResizeTo(aDerivatives);
296  if (transFlag) {
297  localDerivatives = measTransformation * aDerivatives;
298  } else {
299  localDerivatives = aDerivatives;
300  }
301  }
302 }
303 
305 unsigned int GblPoint::getNumLocals() const {
306  return localDerivatives.GetNcols();
307 }
308 
310 const TMatrixD& GblPoint::getLocalDerivatives() const {
311  return localDerivatives;
312 }
313 
315 
320 void GblPoint::addGlobals(const std::vector<int> &aLabels,
321  const TMatrixD &aDerivatives) {
322  if (measDim) {
323  globalLabels = aLabels;
324  globalDerivatives.ResizeTo(aDerivatives);
325  if (transFlag) {
326  globalDerivatives = measTransformation * aDerivatives;
327  } else {
328  globalDerivatives = aDerivatives;
329  }
330 
331  }
332 }
333 
335 unsigned int GblPoint::getNumGlobals() const {
336  return globalDerivatives.GetNcols();
337 }
338 
340 std::vector<int> GblPoint::getGlobalLabels() const {
341  return globalLabels;
342 }
343 
345 const TMatrixD& GblPoint::getGlobalDerivatives() const {
346  return globalDerivatives;
347 }
348 
350 
353 void GblPoint::setLabel(unsigned int aLabel) {
354  theLabel = aLabel;
355 }
356 
358 unsigned int GblPoint::getLabel() const {
359  return theLabel;
360 }
361 
363 
366 void GblPoint::setOffset(int anOffset) {
367  theOffset = anOffset;
368 }
369 
371 int GblPoint::getOffset() const {
372  return theOffset;
373 }
374 
377  return p2pJacobian;
378 }
379 
381 
385  int ifail = 0;
386 // to optimize: need only two last rows of inverse
387 // prevJacobian = aJac.InverseFast(ifail);
388 // block matrix algebra
389  SMatrix23 CA = aJac.Sub<SMatrix23>(3, 0)
390  * aJac.Sub<SMatrix33>(0, 0).InverseFast(ifail); // C*A^-1
391  SMatrix22 DCAB = aJac.Sub<SMatrix22>(3, 3) - CA * aJac.Sub<SMatrix32>(0, 3); // D - C*A^-1 *B
392  DCAB.InvertFast();
393  prevJacobian.Place_at(DCAB, 3, 3);
394  prevJacobian.Place_at(-DCAB * CA, 3, 0);
395 }
396 
398 
402  nextJacobian = aJac;
403 }
404 
406 
415 void GblPoint::getDerivatives(int aDirection, SMatrix22 &matW, SMatrix22 &matWJ,
416  SVector2 &vecWd) const {
417 
418  SMatrix22 matJ;
419  SVector2 vecd;
420  if (aDirection < 1) {
421  matJ = prevJacobian.Sub<SMatrix22>(3, 3);
422  matW = -prevJacobian.Sub<SMatrix22>(3, 1);
423  vecd = prevJacobian.SubCol<SVector2>(0, 3);
424  } else {
425  matJ = nextJacobian.Sub<SMatrix22>(3, 3);
426  matW = nextJacobian.Sub<SMatrix22>(3, 1);
427  vecd = nextJacobian.SubCol<SVector2>(0, 3);
428  }
429 
430  if (!matW.InvertFast()) {
431  std::cout << " GblPoint::getDerivatives failed to invert matrix: "
432  << matW << std::endl;
433  std::cout
434  << " Possible reason for singular matrix: multiple GblPoints at same arc-length"
435  << std::endl;
436  throw std::overflow_error("Singular matrix inversion exception");
437  }
438  matWJ = matW * matJ;
439  vecWd = matW * vecd;
440 
441 }
442 
444 
447 void GblPoint::printPoint(unsigned int level) const {
448  std::cout << " GblPoint";
449  if (theLabel) {
450  std::cout << ", label " << theLabel;
451  if (theOffset >= 0) {
452  std::cout << ", offset " << theOffset;
453  }
454  }
455  if (measDim) {
456  std::cout << ", " << measDim << " measurements";
457  }
458  if (scatFlag) {
459  std::cout << ", scatterer";
460  }
461  if (transFlag) {
462  std::cout << ", diagonalized";
463  }
464  if (localDerivatives.GetNcols()) {
465  std::cout << ", " << localDerivatives.GetNcols()
466  << " local derivatives";
467  }
468  if (globalDerivatives.GetNcols()) {
469  std::cout << ", " << globalDerivatives.GetNcols()
470  << " global derivatives";
471  }
472  std::cout << std::endl;
473  if (level > 0) {
474  if (measDim) {
475  std::cout << " Measurement" << std::endl;
476  std::cout << " Projection: " << std::endl << measProjection
477  << std::endl;
478  std::cout << " Residuals: " << measResiduals << std::endl;
479  std::cout << " Precision: " << measPrecision << std::endl;
480  }
481  if (scatFlag) {
482  std::cout << " Scatterer" << std::endl;
483  std::cout << " Residuals: " << scatResiduals << std::endl;
484  std::cout << " Precision: " << scatPrecision << std::endl;
485  }
486  if (localDerivatives.GetNcols()) {
487  std::cout << " Local Derivatives:" << std::endl;
488  localDerivatives.Print();
489  }
490  if (globalDerivatives.GetNcols()) {
491  std::cout << " Global Labels:";
492  for (unsigned int i = 0; i < globalLabels.size(); ++i) {
493  std::cout << " " << globalLabels[i];
494  }
495  std::cout << std::endl;
496  std::cout << " Global Derivatives:" << std::endl;
497  globalDerivatives.Print();
498  }
499  std::cout << " Jacobian " << std::endl;
500  std::cout << " Point-to-point " << std::endl << p2pJacobian
501  << std::endl;
502  if (theLabel) {
503  std::cout << " To previous offset " << std::endl << prevJacobian
504  << std::endl;
505  std::cout << " To next offset " << std::endl << nextJacobian
506  << std::endl;
507  }
508  }
509 }
510 
511 }