Belle II Software  release-05-01-25
GblTrajectory.cc
Go to the documentation of this file.
1 /*
2  * GblTrajectory.cpp
3  *
4  * Created on: Aug 18, 2011
5  * Author: kleinwrt
6  */
7 
133 #include "GblTrajectory.h"
134 
136 namespace gbl {
137 
139 
147 GblTrajectory::GblTrajectory(const std::vector<GblPoint> &aPointList,
148  bool flagCurv, bool flagU1dir, bool flagU2dir) :
149  numAllPoints(aPointList.size()), numPoints(), numOffsets(0), numInnerTrans(
150  0), numCurvature(flagCurv ? 1 : 0), numParameters(0), numLocals(
151  0), numMeasurements(0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
152 
153  if (flagU1dir)
154  theDimension.push_back(0);
155  if (flagU2dir)
156  theDimension.push_back(1);
157  // simple (single) trajectory
158  thePoints.push_back(aPointList);
159  numPoints.push_back(numAllPoints);
160  construct(); // construct trajectory
161 }
162 
164 
175 GblTrajectory::GblTrajectory(const std::vector<GblPoint> &aPointList,
176  unsigned int aLabel, const TMatrixDSym &aSeed, bool flagCurv,
177  bool flagU1dir, bool flagU2dir) :
178  numAllPoints(aPointList.size()), numPoints(), numOffsets(0), numInnerTrans(
179  0), numCurvature(flagCurv ? 1 : 0), numParameters(0), numLocals(
180  0), numMeasurements(0), externalPoint(aLabel), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(
181  aSeed), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
182 
183  if (flagU1dir)
184  theDimension.push_back(0);
185  if (flagU2dir)
186  theDimension.push_back(1);
187  // simple (single) trajectory
188  thePoints.push_back(aPointList);
189  numPoints.push_back(numAllPoints);
190  construct(); // construct trajectory
191 }
192 
194 
198 GblTrajectory::GblTrajectory(
199  const std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > &aPointsAndTransList) :
200  numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
201  aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
202  0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
203 
204  for (unsigned int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
205  thePoints.push_back(aPointsAndTransList[iTraj].first);
206  numPoints.push_back(thePoints.back().size());
207  numAllPoints += numPoints.back();
208  innerTransformations.push_back(aPointsAndTransList[iTraj].second);
209  }
210  theDimension.push_back(0);
211  theDimension.push_back(1);
212  numCurvature = innerTransformations[0].GetNcols();
213  construct(); // construct (composed) trajectory
214 }
215 
217 
224 GblTrajectory::GblTrajectory(
225  const std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > &aPointsAndTransList,
226  const TMatrixD &extDerivatives, const TVectorD &extMeasurements,
227  const TVectorD &extPrecisions) :
228  numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
229  aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
230  0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(
231  extDerivatives), externalMeasurements(extMeasurements), externalPrecisions(
232  extPrecisions) {
233 
234  for (unsigned int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
235  thePoints.push_back(aPointsAndTransList[iTraj].first);
236  numPoints.push_back(thePoints.back().size());
237  numAllPoints += numPoints.back();
238  innerTransformations.push_back(aPointsAndTransList[iTraj].second);
239  }
240  theDimension.push_back(0);
241  theDimension.push_back(1);
242  numCurvature = innerTransformations[0].GetNcols();
243  construct(); // construct (composed) trajectory
244 }
245 
247 
254 GblTrajectory::GblTrajectory(
255  const std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > &aPointsAndTransList,
256  const TMatrixD &extDerivatives, const TVectorD &extMeasurements,
257  const TMatrixDSym &extPrecisions) :
258  numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
259  aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
260  0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations() {
261 
262  // diagonalize external measurement
263  TMatrixDSymEigen extEigen(extPrecisions);
264  TMatrixD extTransformation = extEigen.GetEigenVectors();
265  extTransformation.T();
266  externalDerivatives.ResizeTo(extDerivatives);
267  externalDerivatives = extTransformation * extDerivatives;
268  externalMeasurements.ResizeTo(extMeasurements);
269  externalMeasurements = extTransformation * extMeasurements;
270  externalPrecisions.ResizeTo(extMeasurements);
271  externalPrecisions = extEigen.GetEigenValues();
272 
273  for (unsigned int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
274  thePoints.push_back(aPointsAndTransList[iTraj].first);
275  numPoints.push_back(thePoints.back().size());
276  numAllPoints += numPoints.back();
277  innerTransformations.push_back(aPointsAndTransList[iTraj].second);
278  }
279  theDimension.push_back(0);
280  theDimension.push_back(1);
281  numCurvature = innerTransformations[0].GetNcols();
282  construct(); // construct (composed) trajectory
283 }
284 
285 GblTrajectory::~GblTrajectory() {
286 }
287 
289 bool GblTrajectory::isValid() const {
290  return constructOK;
291 }
292 
294 unsigned int GblTrajectory::getNumPoints() const {
295  return numAllPoints;
296 }
297 
299 
302 void GblTrajectory::construct() {
303 
304  constructOK = false;
305  fitOK = false;
306  unsigned int aLabel = 0;
307  if (numAllPoints < 2) {
308  std::cout << " GblTrajectory construction failed: too few GblPoints "
309  << std::endl;
310  return;
311  }
312  // loop over trajectories
313  numTrajectories = thePoints.size();
314  //std::cout << " numTrajectories: " << numTrajectories << ", " << innerTransformations.size() << std::endl;
315  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
316  std::vector<GblPoint>::iterator itPoint;
317  for (itPoint = thePoints[iTraj].begin();
318  itPoint < thePoints[iTraj].end(); ++itPoint) {
319  numLocals = std::max(numLocals, itPoint->getNumLocals());
320  numMeasurements += itPoint->hasMeasurement();
321  itPoint->setLabel(++aLabel);
322  }
323  }
324  defineOffsets();
325  calcJacobians();
326  try {
327  prepare();
328  } catch (std::overflow_error &e) {
329  std::cout << " GblTrajectory construction failed: " << e.what()
330  << std::endl;
331  return;
332  }
333  constructOK = true;
334  // number of fit parameters
335  numParameters = (numOffsets - 2 * numInnerTrans) * theDimension.size()
336  + numCurvature + numLocals;
337 }
338 
340 
344 void GblTrajectory::defineOffsets() {
345 
346  // loop over trajectories
347  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
348  // first point is offset
349  thePoints[iTraj].front().setOffset(numOffsets++);
350  // intermediate scatterers are offsets
351  std::vector<GblPoint>::iterator itPoint;
352  for (itPoint = thePoints[iTraj].begin() + 1;
353  itPoint < thePoints[iTraj].end() - 1; ++itPoint) {
354  if (itPoint->hasScatterer()) {
355  itPoint->setOffset(numOffsets++);
356  } else {
357  itPoint->setOffset(-numOffsets);
358  }
359  }
360  // last point is offset
361  thePoints[iTraj].back().setOffset(numOffsets++);
362  }
363 }
364 
366 void GblTrajectory::calcJacobians() {
367 
368  SMatrix55 scatJacobian;
369  // loop over trajectories
370  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
371  // forward propagation (all)
372  GblPoint* previousPoint = &thePoints[iTraj].front();
373  unsigned int numStep = 0;
374  std::vector<GblPoint>::iterator itPoint;
375  for (itPoint = thePoints[iTraj].begin() + 1;
376  itPoint < thePoints[iTraj].end(); ++itPoint) {
377  if (numStep == 0) {
378  scatJacobian = itPoint->getP2pJacobian();
379  } else {
380  scatJacobian = itPoint->getP2pJacobian() * scatJacobian;
381  }
382  numStep++;
383  itPoint->addPrevJacobian(scatJacobian); // iPoint -> previous scatterer
384  if (itPoint->getOffset() >= 0) {
385  previousPoint->addNextJacobian(scatJacobian); // lastPoint -> next scatterer
386  numStep = 0;
387  previousPoint = &(*itPoint);
388  }
389  }
390  // backward propagation (without scatterers)
391  for (itPoint = thePoints[iTraj].end() - 1;
392  itPoint > thePoints[iTraj].begin(); --itPoint) {
393  if (itPoint->getOffset() >= 0) {
394  scatJacobian = itPoint->getP2pJacobian();
395  continue; // skip offsets
396  }
397  itPoint->addNextJacobian(scatJacobian); // iPoint -> next scatterer
398  scatJacobian = scatJacobian * itPoint->getP2pJacobian();
399  }
400  }
401 }
402 
404 
412 std::pair<std::vector<unsigned int>, TMatrixD> GblTrajectory::getJacobian(
413  int aSignedLabel) const {
414 
415  unsigned int nDim = theDimension.size();
416  unsigned int nCurv = numCurvature;
417  unsigned int nLocals = numLocals;
418  unsigned int nBorder = nCurv + nLocals;
419  unsigned int nParBRL = nBorder + 2 * nDim;
420  unsigned int nParLoc = nLocals + 5;
421  std::vector<unsigned int> anIndex;
422  anIndex.reserve(nParBRL);
423  TMatrixD aJacobian(nParLoc, nParBRL);
424  aJacobian.Zero();
425 
426  unsigned int aLabel = abs(aSignedLabel);
427  unsigned int firstLabel = 1;
428  unsigned int lastLabel = 0;
429  unsigned int aTrajectory = 0;
430  // loop over trajectories
431  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
432  aTrajectory = iTraj;
433  lastLabel += numPoints[iTraj];
434  if (aLabel <= lastLabel)
435  break;
436  if (iTraj < numTrajectories - 1)
437  firstLabel += numPoints[iTraj];
438  }
439  int nJacobian; // 0: prev, 1: next
440  // check consistency of (index, direction)
441  if (aSignedLabel > 0) {
442  nJacobian = 1;
443  if (aLabel >= lastLabel) {
444  aLabel = lastLabel;
445  nJacobian = 0;
446  }
447  } else {
448  nJacobian = 0;
449  if (aLabel <= firstLabel) {
450  aLabel = firstLabel;
451  nJacobian = 1;
452  }
453  }
454  const GblPoint aPoint = thePoints[aTrajectory][aLabel - firstLabel];
455  std::vector<unsigned int> labDer(5);
456  SMatrix55 matDer;
457  getFitToLocalJacobian(labDer, matDer, aPoint, 5, nJacobian);
458 
459  // from local parameters
460  for (unsigned int i = 0; i < nLocals; ++i) {
461  aJacobian(i + 5, i) = 1.0;
462  anIndex.push_back(i + 1);
463  }
464  // from trajectory parameters
465  unsigned int iCol = nLocals;
466  for (unsigned int i = 0; i < 5; ++i) {
467  if (labDer[i] > 0) {
468  anIndex.push_back(labDer[i]);
469  for (unsigned int j = 0; j < 5; ++j) {
470  aJacobian(j, iCol) = matDer(j, i);
471  }
472  ++iCol;
473  }
474  }
475  return std::make_pair(anIndex, aJacobian);
476 }
477 
479 
488 void GblTrajectory::getFitToLocalJacobian(std::vector<unsigned int> &anIndex,
489  SMatrix55 &aJacobian, const GblPoint &aPoint, unsigned int measDim,
490  unsigned int nJacobian) const {
491 
492  unsigned int nDim = theDimension.size();
493  unsigned int nCurv = numCurvature;
494  unsigned int nLocals = numLocals;
495 
496  int nOffset = aPoint.getOffset();
497 
498  if (nOffset < 0) // need interpolation
499  {
500  SMatrix22 prevW, prevWJ, nextW, nextWJ, matN;
501  SVector2 prevWd, nextWd;
502  int ierr;
503  aPoint.getDerivatives(0, prevW, prevWJ, prevWd); // W-, W- * J-, W- * d-
504  aPoint.getDerivatives(1, nextW, nextWJ, nextWd); // W-, W- * J-, W- * d-
505  const SMatrix22 sumWJ(prevWJ + nextWJ);
506  matN = sumWJ.Inverse(ierr); // N = (W- * J- + W+ * J+)^-1
507  // derivatives for u_int
508  const SMatrix22 prevNW(matN * prevW); // N * W-
509  const SMatrix22 nextNW(matN * nextW); // N * W+
510  const SVector2 prevNd(matN * prevWd); // N * W- * d-
511  const SVector2 nextNd(matN * nextWd); // N * W+ * d+
512 
513  unsigned int iOff = nDim * (-nOffset - 1) + nLocals + nCurv + 1; // first offset ('i' in u_i)
514 
515  // local offset
516  if (nCurv > 0) {
517  aJacobian.Place_in_col(-prevNd - nextNd, 3, 0); // from curvature
518  anIndex[0] = nLocals + 1;
519  }
520  aJacobian.Place_at(prevNW, 3, 1); // from 1st Offset
521  aJacobian.Place_at(nextNW, 3, 3); // from 2nd Offset
522  for (unsigned int i = 0; i < nDim; ++i) {
523  anIndex[1 + theDimension[i]] = iOff + i;
524  anIndex[3 + theDimension[i]] = iOff + nDim + i;
525  }
526 
527  // local slope and curvature
528  if (measDim > 2) {
529  // derivatives for u'_int
530  const SMatrix22 prevWPN(nextWJ * prevNW); // W+ * J+ * N * W-
531  const SMatrix22 nextWPN(prevWJ * nextNW); // W- * J- * N * W+
532  const SVector2 prevWNd(nextWJ * prevNd); // W+ * J+ * N * W- * d-
533  const SVector2 nextWNd(prevWJ * nextNd); // W- * J- * N * W+ * d+
534  if (nCurv > 0) {
535  aJacobian(0, 0) = 1.0;
536  aJacobian.Place_in_col(prevWNd - nextWNd, 1, 0); // from curvature
537  }
538  aJacobian.Place_at(-prevWPN, 1, 1); // from 1st Offset
539  aJacobian.Place_at(nextWPN, 1, 3); // from 2nd Offset
540  }
541  } else { // at point
542  // anIndex must be sorted
543  // forward : iOff2 = iOff1 + nDim, index1 = 1, index2 = 3
544  // backward: iOff2 = iOff1 - nDim, index1 = 3, index2 = 1
545  unsigned int iOff1 = nDim * nOffset + nCurv + nLocals + 1; // first offset ('i' in u_i)
546  unsigned int index1 = 3 - 2 * nJacobian; // index of first offset
547  unsigned int iOff2 = iOff1 + nDim * (nJacobian * 2 - 1); // second offset ('i' in u_i)
548  unsigned int index2 = 1 + 2 * nJacobian; // index of second offset
549  // local offset
550  aJacobian(3, index1) = 1.0; // from 1st Offset
551  aJacobian(4, index1 + 1) = 1.0;
552  for (unsigned int i = 0; i < nDim; ++i) {
553  anIndex[index1 + theDimension[i]] = iOff1 + i;
554  }
555 
556  // local slope and curvature
557  if (measDim > 2) {
558  SMatrix22 matW, matWJ;
559  SVector2 vecWd;
560  aPoint.getDerivatives(nJacobian, matW, matWJ, vecWd); // W, W * J, W * d
561  double sign = (nJacobian > 0) ? 1. : -1.;
562  if (nCurv > 0) {
563  aJacobian(0, 0) = 1.0;
564  aJacobian.Place_in_col(-sign * vecWd, 1, 0); // from curvature
565  anIndex[0] = nLocals + 1;
566  }
567  aJacobian.Place_at(-sign * matWJ, 1, index1); // from 1st Offset
568  aJacobian.Place_at(sign * matW, 1, index2); // from 2nd Offset
569  for (unsigned int i = 0; i < nDim; ++i) {
570  anIndex[index2 + theDimension[i]] = iOff2 + i;
571  }
572  }
573  }
574 }
575 
577 
583 void GblTrajectory::getFitToKinkJacobian(std::vector<unsigned int> &anIndex,
584  SMatrix27 &aJacobian, const GblPoint &aPoint) const {
585 
586  unsigned int nDim = theDimension.size();
587  unsigned int nCurv = numCurvature;
588  unsigned int nLocals = numLocals;
589 
590  int nOffset = aPoint.getOffset();
591 
592  SMatrix22 prevW, prevWJ, nextW, nextWJ;
593  SVector2 prevWd, nextWd;
594  aPoint.getDerivatives(0, prevW, prevWJ, prevWd); // W-, W- * J-, W- * d-
595  aPoint.getDerivatives(1, nextW, nextWJ, nextWd); // W-, W- * J-, W- * d-
596  const SMatrix22 sumWJ(prevWJ + nextWJ); // W- * J- + W+ * J+
597  const SVector2 sumWd(prevWd + nextWd); // W+ * d+ + W- * d-
598 
599  unsigned int iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1; // first offset ('i' in u_i)
600 
601  // local offset
602  if (nCurv > 0) {
603  aJacobian.Place_in_col(-sumWd, 0, 0); // from curvature
604  anIndex[0] = nLocals + 1;
605  }
606  aJacobian.Place_at(prevW, 0, 1); // from 1st Offset
607  aJacobian.Place_at(-sumWJ, 0, 3); // from 2nd Offset
608  aJacobian.Place_at(nextW, 0, 5); // from 1st Offset
609  for (unsigned int i = 0; i < nDim; ++i) {
610  anIndex[1 + theDimension[i]] = iOff + i;
611  anIndex[3 + theDimension[i]] = iOff + nDim + i;
612  anIndex[5 + theDimension[i]] = iOff + nDim * 2 + i;
613  }
614 }
615 
617 
631 unsigned int GblTrajectory::getResults(int aSignedLabel, TVectorD &localPar,
632  TMatrixDSym &localCov) const {
633  if (not fitOK)
634  return 1;
635  std::pair<std::vector<unsigned int>, TMatrixD> indexAndJacobian =
636  getJacobian(aSignedLabel);
637  unsigned int nParBrl = indexAndJacobian.first.size();
638  TVectorD aVec(nParBrl); // compressed vector
639  for (unsigned int i = 0; i < nParBrl; ++i) {
640  aVec[i] = theVector(indexAndJacobian.first[i] - 1);
641  }
642  TMatrixDSym aMat = theMatrix.getBlockMatrix(indexAndJacobian.first); // compressed matrix
643  localPar = indexAndJacobian.second * aVec;
644  localCov = aMat.Similarity(indexAndJacobian.second);
645  return 0;
646 }
647 
649 
661 unsigned int GblTrajectory::getMeasResults(unsigned int aLabel,
662  unsigned int &numData, TVectorD &aResiduals, TVectorD &aMeasErrors,
663  TVectorD &aResErrors, TVectorD &aDownWeights) {
664  numData = 0;
665  if (not fitOK)
666  return 1;
667 
668  unsigned int firstData = measDataIndex[aLabel - 1]; // first data block with measurement
669  numData = measDataIndex[aLabel] - firstData; // number of data blocks
670  for (unsigned int i = 0; i < numData; ++i) {
671  getResAndErr(firstData + i, aResiduals[i], aMeasErrors[i],
672  aResErrors[i], aDownWeights[i]);
673  }
674  return 0;
675 }
676 
678 
690 unsigned int GblTrajectory::getScatResults(unsigned int aLabel,
691  unsigned int &numData, TVectorD &aResiduals, TVectorD &aMeasErrors,
692  TVectorD &aResErrors, TVectorD &aDownWeights) {
693  numData = 0;
694  if (not fitOK)
695  return 1;
696 
697  unsigned int firstData = scatDataIndex[aLabel - 1]; // first data block with scatterer
698  numData = scatDataIndex[aLabel] - firstData; // number of data blocks
699  for (unsigned int i = 0; i < numData; ++i) {
700  getResAndErr(firstData + i, aResiduals[i], aMeasErrors[i],
701  aResErrors[i], aDownWeights[i]);
702  }
703  return 0;
704 }
705 
707 
711 unsigned int GblTrajectory::getLabels(std::vector<unsigned int> &aLabelList) {
712  if (not constructOK)
713  return 1;
714 
715  unsigned int aLabel = 0;
716  unsigned int nPoint = thePoints[0].size();
717  aLabelList.resize(nPoint);
718  for (unsigned i = 0; i < nPoint; ++i) {
719  aLabelList[i] = ++aLabel;
720  }
721  return 0;
722 }
723 
725 
729 unsigned int GblTrajectory::getLabels(
730  std::vector<std::vector<unsigned int> > &aLabelList) {
731  if (not constructOK)
732  return 1;
733 
734  unsigned int aLabel = 0;
735  aLabelList.resize(numTrajectories);
736  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
737  unsigned int nPoint = thePoints[iTraj].size();
738  aLabelList[iTraj].resize(nPoint);
739  for (unsigned i = 0; i < nPoint; ++i) {
740  aLabelList[iTraj][i] = ++aLabel;
741  }
742  }
743  return 0;
744 }
745 
747 
756 void GblTrajectory::getResAndErr(unsigned int aData, double &aResidual,
757  double &aMeasError, double &aResError, double &aDownWeight) {
758 
759  double aMeasVar;
760  std::vector<unsigned int>* indLocal;
761  std::vector<double>* derLocal;
762  theData[aData].getResidual(aResidual, aMeasVar, aDownWeight, indLocal,
763  derLocal);
764  unsigned int nParBrl = (*indLocal).size();
765  TVectorD aVec(nParBrl); // compressed vector of derivatives
766  for (unsigned int j = 0; j < nParBrl; ++j) {
767  aVec[j] = (*derLocal)[j];
768  }
769  TMatrixDSym aMat = theMatrix.getBlockMatrix(*indLocal); // compressed (covariance) matrix
770  double aFitVar = aMat.Similarity(aVec); // variance from track fit
771  aMeasError = sqrt(aMeasVar); // error of measurement
772  aResError = (aFitVar < aMeasVar ? sqrt(aMeasVar - aFitVar) : 0.); // error of residual
773 }
774 
776 void GblTrajectory::buildLinearEquationSystem() {
777  unsigned int nBorder = numCurvature + numLocals;
778  theVector.resize(numParameters);
779  theMatrix.resize(numParameters, nBorder);
780  double aValue, aWeight;
781  std::vector<unsigned int>* indLocal;
782  std::vector<double>* derLocal;
783  std::vector<GblData>::iterator itData;
784  for (itData = theData.begin(); itData < theData.end(); ++itData) {
785  itData->getLocalData(aValue, aWeight, indLocal, derLocal);
786  for (unsigned int j = 0; j < indLocal->size(); ++j) {
787  theVector((*indLocal)[j] - 1) += (*derLocal)[j] * aWeight * aValue;
788  }
789  theMatrix.addBlockMatrix(aWeight, indLocal, derLocal);
790  }
791 }
792 
794 
797 void GblTrajectory::prepare() {
798  unsigned int nDim = theDimension.size();
799  // upper limit
800  unsigned int maxData = numMeasurements + nDim * (numOffsets - 2)
801  + externalSeed.GetNrows();
802  theData.reserve(maxData);
803  measDataIndex.resize(numAllPoints + 3); // include external seed and measurements
804  scatDataIndex.resize(numAllPoints + 1);
805  unsigned int nData = 0;
806  std::vector<TMatrixD> innerTransDer;
807  std::vector<std::vector<unsigned int> > innerTransLab;
808  // composed trajectory ?
809  if (numInnerTrans > 0) {
810  //std::cout << "composed trajectory" << std::endl;
811  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
812  // innermost point
813  GblPoint* innerPoint = &thePoints[iTraj].front();
814  // transformation fit to local track parameters
815  std::vector<unsigned int> firstLabels(5);
816  SMatrix55 matFitToLocal, matLocalToFit;
817  getFitToLocalJacobian(firstLabels, matFitToLocal, *innerPoint, 5);
818  // transformation local track to fit parameters
819  int ierr;
820  matLocalToFit = matFitToLocal.Inverse(ierr);
821  TMatrixD localToFit(5, 5);
822  for (unsigned int i = 0; i < 5; ++i) {
823  for (unsigned int j = 0; j < 5; ++j) {
824  localToFit(i, j) = matLocalToFit(i, j);
825  }
826  }
827  // transformation external to fit parameters at inner (first) point
828  innerTransDer.push_back(localToFit * innerTransformations[iTraj]);
829  innerTransLab.push_back(firstLabels);
830  }
831  }
832  // measurements
833  SMatrix55 matP;
834  // loop over trajectories
835  std::vector<GblPoint>::iterator itPoint;
836  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
837  for (itPoint = thePoints[iTraj].begin();
838  itPoint < thePoints[iTraj].end(); ++itPoint) {
839  SVector5 aMeas, aPrec;
840  unsigned int nLabel = itPoint->getLabel();
841  unsigned int measDim = itPoint->hasMeasurement();
842  if (measDim) {
843  const TMatrixD localDer = itPoint->getLocalDerivatives();
844  const std::vector<int> globalLab = itPoint->getGlobalLabels();
845  const TMatrixD globalDer = itPoint->getGlobalDerivatives();
846  TMatrixD transDer;
847  itPoint->getMeasurement(matP, aMeas, aPrec);
848  unsigned int iOff = 5 - measDim; // first active component
849  std::vector<unsigned int> labDer(5);
850  SMatrix55 matDer, matPDer;
851  unsigned int nJacobian =
852  (itPoint < thePoints[iTraj].end() - 1) ? 1 : 0; // last point needs backward propagation
853  getFitToLocalJacobian(labDer, matDer, *itPoint, measDim,
854  nJacobian);
855  if (measDim > 2) {
856  matPDer = matP * matDer;
857  } else { // 'shortcut' for position measurements
858  matPDer.Place_at(
859  matP.Sub<SMatrix22>(3, 3)
860  * matDer.Sub<SMatrix25>(3, 0), 3, 0);
861  }
862 
863  if (numInnerTrans > 0) {
864  // transform for external parameters
865  TMatrixD proDer(measDim, 5);
866  // match parameters
867  unsigned int ifirst = 0;
868  unsigned int ilabel = 0;
869  while (ilabel < 5) {
870  if (labDer[ilabel] > 0) {
871  while (innerTransLab[iTraj][ifirst]
872  != labDer[ilabel] and ifirst < 5) {
873  ++ifirst;
874  }
875  if (ifirst >= 5) {
876  labDer[ilabel] -= 2 * nDim * (iTraj + 1); // adjust label
877  } else {
878  // match
879  labDer[ilabel] = 0; // mark as related to external parameters
880  for (unsigned int k = iOff; k < 5; ++k) {
881  proDer(k - iOff, ifirst) = matPDer(k,
882  ilabel);
883  }
884  }
885  }
886  ++ilabel;
887  }
888  transDer.ResizeTo(measDim, numCurvature);
889  transDer = proDer * innerTransDer[iTraj];
890  }
891  for (unsigned int i = iOff; i < 5; ++i) {
892  if (aPrec(i) > 0.) {
893  GblData aData(nLabel, aMeas(i), aPrec(i));
894  aData.addDerivatives(i, labDer, matPDer, iOff, localDer,
895  globalLab, globalDer, numLocals, transDer);
896  theData.push_back(aData);
897  nData++;
898  }
899  }
900 
901  }
902  measDataIndex[nLabel] = nData;
903  }
904  }
905 
906  // pseudo measurements from kinks
907  SMatrix22 matT;
908  scatDataIndex[0] = nData;
909  scatDataIndex[1] = nData;
910  // loop over trajectories
911  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
912  for (itPoint = thePoints[iTraj].begin() + 1;
913  itPoint < thePoints[iTraj].end() - 1; ++itPoint) {
914  SVector2 aMeas, aPrec;
915  unsigned int nLabel = itPoint->getLabel();
916  if (itPoint->hasScatterer()) {
917  itPoint->getScatterer(matT, aMeas, aPrec);
918  TMatrixD transDer;
919  std::vector<unsigned int> labDer(7);
920  SMatrix27 matDer, matTDer;
921  getFitToKinkJacobian(labDer, matDer, *itPoint);
922  matTDer = matT * matDer;
923  if (numInnerTrans > 0) {
924  // transform for external parameters
925  TMatrixD proDer(nDim, 5);
926  // match parameters
927  unsigned int ifirst = 0;
928  unsigned int ilabel = 0;
929  while (ilabel < 7) {
930  if (labDer[ilabel] > 0) {
931  while (innerTransLab[iTraj][ifirst]
932  != labDer[ilabel] and ifirst < 5) {
933  ++ifirst;
934  }
935  if (ifirst >= 5) {
936  labDer[ilabel] -= 2 * nDim * (iTraj + 1); // adjust label
937  } else {
938  // match
939  labDer[ilabel] = 0; // mark as related to external parameters
940  for (unsigned int k = 0; k < nDim; ++k) {
941  proDer(k, ifirst) = matTDer(k, ilabel);
942  }
943  }
944  }
945  ++ilabel;
946  }
947  transDer.ResizeTo(nDim, numCurvature);
948  transDer = proDer * innerTransDer[iTraj];
949  }
950  for (unsigned int i = 0; i < nDim; ++i) {
951  unsigned int iDim = theDimension[i];
952  if (aPrec(iDim) > 0.) {
953  GblData aData(nLabel, aMeas(iDim), aPrec(iDim));
954  aData.addDerivatives(iDim, labDer, matTDer, numLocals,
955  transDer);
956  theData.push_back(aData);
957  nData++;
958  }
959  }
960  }
961  scatDataIndex[nLabel] = nData;
962  }
963  scatDataIndex[thePoints[iTraj].back().getLabel()] = nData;
964  }
965 
966  // external seed
967  if (externalPoint > 0) {
968  std::pair<std::vector<unsigned int>, TMatrixD> indexAndJacobian =
969  getJacobian(externalPoint);
970  std::vector<unsigned int> externalSeedIndex = indexAndJacobian.first;
971  std::vector<double> externalSeedDerivatives(externalSeedIndex.size());
972  const TMatrixDSymEigen externalSeedEigen(externalSeed);
973  const TVectorD valEigen = externalSeedEigen.GetEigenValues();
974  TMatrixD vecEigen = externalSeedEigen.GetEigenVectors();
975  vecEigen = vecEigen.T() * indexAndJacobian.second;
976  for (int i = 0; i < externalSeed.GetNrows(); ++i) {
977  if (valEigen(i) > 0.) {
978  for (int j = 0; j < externalSeed.GetNcols(); ++j) {
979  externalSeedDerivatives[j] = vecEigen(i, j);
980  }
981  GblData aData(externalPoint, 0., valEigen(i));
982  aData.addDerivatives(externalSeedIndex,
983  externalSeedDerivatives);
984  theData.push_back(aData);
985  nData++;
986  }
987  }
988  }
989  measDataIndex[numAllPoints + 1] = nData;
990  // external measurements
991  unsigned int nExt = externalMeasurements.GetNrows();
992  if (nExt > 0) {
993  std::vector<unsigned int> index(numCurvature);
994  std::vector<double> derivatives(numCurvature);
995  for (unsigned int iExt = 0; iExt < nExt; ++iExt) {
996  for (unsigned int iCol = 0; iCol < numCurvature; ++iCol) {
997  index[iCol] = numLocals + iCol + 1;
998  derivatives[iCol] = externalDerivatives(iExt, iCol);
999  }
1000  GblData aData(1U, externalMeasurements(iExt),
1001  externalPrecisions(iExt));
1002  aData.addDerivatives(index, derivatives);
1003  theData.push_back(aData);
1004  nData++;
1005  }
1006  }
1007  measDataIndex[numAllPoints + 2] = nData;
1008 }
1009 
1011 void GblTrajectory::predict() {
1012  std::vector<GblData>::iterator itData;
1013  for (itData = theData.begin(); itData < theData.end(); ++itData) {
1014  itData->setPrediction(theVector);
1015  }
1016 }
1017 
1019 
1022 double GblTrajectory::downWeight(unsigned int aMethod) {
1023  double aLoss = 0.;
1024  std::vector<GblData>::iterator itData;
1025  for (itData = theData.begin(); itData < theData.end(); ++itData) {
1026  aLoss += (1. - itData->setDownWeighting(aMethod));
1027  }
1028  return aLoss;
1029 }
1030 
1032 
1043 unsigned int GblTrajectory::fit(double &Chi2, int &Ndf, double &lostWeight,
1044  std::string optionList) {
1045  const double normChi2[4] = { 1.0, 0.8737, 0.9326, 0.8228 };
1046  const std::string methodList = "TtHhCc";
1047 
1048  Chi2 = 0.;
1049  Ndf = -1;
1050  lostWeight = 0.;
1051  if (not constructOK)
1052  return 10;
1053 
1054  unsigned int aMethod = 0;
1055 
1056  buildLinearEquationSystem();
1057  lostWeight = 0.;
1058  unsigned int ierr = 0;
1059  try {
1060 
1061  theMatrix.solveAndInvertBorderedBand(theVector, theVector);
1062  predict();
1063 
1064  for (unsigned int i = 0; i < optionList.size(); ++i) // down weighting iterations
1065  {
1066  size_t aPosition = methodList.find(optionList[i]);
1067  if (aPosition != std::string::npos) {
1068  aMethod = aPosition / 2 + 1;
1069  lostWeight = downWeight(aMethod);
1070  buildLinearEquationSystem();
1071  theMatrix.solveAndInvertBorderedBand(theVector, theVector);
1072  predict();
1073  }
1074  }
1075  Ndf = theData.size() - numParameters;
1076  Chi2 = 0.;
1077  for (unsigned int i = 0; i < theData.size(); ++i) {
1078  Chi2 += theData[i].getChi2();
1079  }
1080  Chi2 /= normChi2[aMethod];
1081  fitOK = true;
1082 
1083  } catch (int e) {
1084  std::cout << " GblTrajectory::fit exception " << e << std::endl;
1085  ierr = e;
1086  }
1087  return ierr;
1088 }
1089 
1091 void GblTrajectory::milleOut(MilleBinary &aMille) {
1092  double aValue;
1093  double aErr;
1094  std::vector<unsigned int>* indLocal;
1095  std::vector<double>* derLocal;
1096  std::vector<int>* labGlobal;
1097  std::vector<double>* derGlobal;
1098 
1099  if (not constructOK)
1100  return;
1101 
1102 // data: measurements, kinks and external seed
1103  std::vector<GblData>::iterator itData;
1104  for (itData = theData.begin(); itData != theData.end(); ++itData) {
1105  itData->getAllData(aValue, aErr, indLocal, derLocal, labGlobal,
1106  derGlobal);
1107  aMille.addData(aValue, aErr, *indLocal, *derLocal, *labGlobal,
1108  *derGlobal);
1109  }
1110  aMille.writeRecord();
1111 }
1112 
1114 
1117 void GblTrajectory::printTrajectory(unsigned int level) {
1118  if (numInnerTrans) {
1119  std::cout << "Composed GblTrajectory, " << numInnerTrans
1120  << " subtrajectories" << std::endl;
1121  } else {
1122  std::cout << "Simple GblTrajectory" << std::endl;
1123  }
1124  if (theDimension.size() < 2) {
1125  std::cout << " 2D-trajectory" << std::endl;
1126  }
1127  std::cout << " Number of GblPoints : " << numAllPoints
1128  << std::endl;
1129  std::cout << " Number of points with offsets: " << numOffsets << std::endl;
1130  std::cout << " Number of fit parameters : " << numParameters
1131  << std::endl;
1132  std::cout << " Number of measurements : " << numMeasurements
1133  << std::endl;
1134  if (externalMeasurements.GetNrows()) {
1135  std::cout << " Number of ext. measurements : "
1136  << externalMeasurements.GetNrows() << std::endl;
1137  }
1138  if (externalPoint) {
1139  std::cout << " Label of point with ext. seed: " << externalPoint
1140  << std::endl;
1141  }
1142  if (constructOK) {
1143  std::cout << " Constructed OK " << std::endl;
1144  }
1145  if (fitOK) {
1146  std::cout << " Fitted OK " << std::endl;
1147  }
1148  if (level > 0) {
1149  if (numInnerTrans) {
1150  std::cout << " Inner transformations" << std::endl;
1151  for (unsigned int i = 0; i < numInnerTrans; ++i) {
1152  innerTransformations[i].Print();
1153  }
1154  }
1155  if (externalMeasurements.GetNrows()) {
1156  std::cout << " External measurements" << std::endl;
1157  std::cout << " Measurements:" << std::endl;
1158  externalMeasurements.Print();
1159  std::cout << " Precisions:" << std::endl;
1160  externalPrecisions.Print();
1161  std::cout << " Derivatives:" << std::endl;
1162  externalDerivatives.Print();
1163  }
1164  if (externalPoint) {
1165  std::cout << " External seed:" << std::endl;
1166  externalSeed.Print();
1167  }
1168  if (fitOK) {
1169  std::cout << " Fit results" << std::endl;
1170  std::cout << " Parameters:" << std::endl;
1171  theVector.print();
1172  std::cout << " Covariance matrix (bordered band part):"
1173  << std::endl;
1174  theMatrix.printMatrix();
1175  }
1176  }
1177 }
1178 
1180 
1183 void GblTrajectory::printPoints(unsigned int level) {
1184  std::cout << "GblPoints " << std::endl;
1185  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
1186  std::vector<GblPoint>::iterator itPoint;
1187  for (itPoint = thePoints[iTraj].begin();
1188  itPoint < thePoints[iTraj].end(); ++itPoint) {
1189  itPoint->printPoint(level);
1190  }
1191  }
1192 }
1193 
1195 void GblTrajectory::printData() {
1196  std::cout << "GblData blocks " << std::endl;
1197  std::vector<GblData>::iterator itData;
1198  for (itData = theData.begin(); itData < theData.end(); ++itData) {
1199  itData->printData();
1200  }
1201 }
1202 
1203 }
gbl::GblTrajectory::getMeasResults
unsigned int getMeasResults(unsigned int aLabel, unsigned int &numRes, TVectorD &aResiduals, TVectorD &aMeasErrors, TVectorD &aResErrors, TVectorD &aDownWeights)
Get residuals from fit at point for measurement.
Definition: GblTrajectory.cc:661
gbl::GblTrajectory::externalPoint
unsigned int externalPoint
Label of external point (or 0)
Definition: GblTrajectory.h:96
gbl::GblTrajectory::getJacobian
std::pair< std::vector< unsigned int >, TMatrixD > getJacobian(int aSignedLabel) const
Get jacobian for transformation from fit to track parameters at point.
Definition: GblTrajectory.cc:412
gbl::GblTrajectory::getResults
unsigned int getResults(int aSignedLabel, TVectorD &localPar, TMatrixDSym &localCov) const
Get fit results at point.
Definition: GblTrajectory.cc:631
gbl::GblTrajectory::milleOut
void milleOut(MilleBinary &aMille)
Write valid trajectory to Millepede-II binary file.
Definition: GblTrajectory.cc:1091
gbl::GblTrajectory::getFitToLocalJacobian
void getFitToLocalJacobian(std::vector< unsigned int > &anIndex, SMatrix55 &aJacobian, const GblPoint &aPoint, unsigned int measDim, unsigned int nJacobian=1) const
Get (part of) jacobian for transformation from (trajectory) fit to track parameters at point.
Definition: GblTrajectory.cc:488
gbl::GblTrajectory::getResAndErr
void getResAndErr(unsigned int aData, double &aResidual, double &aMeadsError, double &aResError, double &aDownWeight)
Get residual and errors from data block.
Definition: GblTrajectory.cc:756
gbl::GblTrajectory::printPoints
void printPoints(unsigned int level=0)
Print GblPoints on trajectory.
Definition: GblTrajectory.cc:1183
gbl::GblPoint::addNextJacobian
void addNextJacobian(const SMatrix55 &aJac)
Define jacobian to next scatterer (by GBLTrajectory constructor)
Definition: GblPoint.cc:401
gbl::GblTrajectory::theMatrix
BorderedBandMatrix theMatrix
(Bordered band) matrix of linear equation system
Definition: GblTrajectory.h:111
gbl::GblTrajectory::numOffsets
unsigned int numOffsets
Number of (points with) offsets on trajectory.
Definition: GblTrajectory.h:90
gbl::GblTrajectory::printData
void printData()
Print GblData blocks for trajectory.
Definition: GblTrajectory.cc:1195
gbl::GblTrajectory::thePoints
std::vector< std::vector< GblPoint > > thePoints
(list of) List of points on trajectory
Definition: GblTrajectory.h:100
gbl::GblTrajectory::fitOK
bool fitOK
Trajectory has been successfully fitted (results are valid)
Definition: GblTrajectory.h:98
gbl::GblTrajectory::getScatResults
unsigned int getScatResults(unsigned int aLabel, unsigned int &numRes, TVectorD &aResiduals, TVectorD &aMeasErrors, TVectorD &aResErrors, TVectorD &aDownWeights)
Get (kink) residuals from fit at point for scatterer.
Definition: GblTrajectory.cc:690
gbl::VVector::resize
void resize(const unsigned int nRows)
Resize vector.
Definition: VMatrix.cc:262
gbl::GblPoint::getDerivatives
void getDerivatives(int aDirection, SMatrix22 &matW, SMatrix22 &matWJ, SVector2 &vecWd) const
Retrieve derivatives of local track model.
Definition: GblPoint.cc:415
gbl
Namespace for the general broken lines package.
Definition: BorderedBandMatrix.h:44
gbl::GblTrajectory::downWeight
double downWeight(unsigned int aMethod)
Down-weight all points.
Definition: GblTrajectory.cc:1022
gbl::GblData::addDerivatives
void addDerivatives(unsigned int iRow, const std::vector< unsigned int > &labDer, const SMatrix55 &matDer, unsigned int iOff, const TMatrixD &derLocal, const std::vector< int > &labGlobal, const TMatrixD &derGlobal, unsigned int nLocal, const TMatrixD &derTrans)
Add derivatives from measurement.
Definition: GblData.cc:63
gbl::GblTrajectory::numPoints
std::vector< unsigned int > numPoints
Number of points on (sub)trajectory.
Definition: GblTrajectory.h:88
gbl::GblTrajectory::numLocals
unsigned int numLocals
Total number of (additional) local parameters.
Definition: GblTrajectory.h:94
gbl::GblTrajectory::measDataIndex
std::vector< unsigned int > measDataIndex
mapping points to data blocks from measurements
Definition: GblTrajectory.h:102
gbl::VVector::print
void print() const
Print vector.
Definition: VMatrix.cc:298
gbl::GblTrajectory::printTrajectory
void printTrajectory(unsigned int level=0)
Print GblTrajectory.
Definition: GblTrajectory.cc:1117
gbl::GblTrajectory::scatDataIndex
std::vector< unsigned int > scatDataIndex
mapping points to data blocks from scatterers
Definition: GblTrajectory.h:103
gbl::BorderedBandMatrix::printMatrix
void printMatrix() const
Print bordered band matrix.
Definition: BorderedBandMatrix.cc:180
gbl::GblTrajectory::numTrajectories
unsigned int numTrajectories
Number of trajectories (in composed trajectory)
Definition: GblTrajectory.h:89
gbl::GblData
Data (block) for independent scalar measurement.
Definition: GblData.h:55
gbl::GblTrajectory::construct
void construct()
Construct trajectory from list of points.
Definition: GblTrajectory.cc:302
gbl::BorderedBandMatrix::addBlockMatrix
void addBlockMatrix(double aWeight, const std::vector< unsigned int > *anIndex, const std::vector< double > *aVector)
Add symmetric block matrix.
Definition: BorderedBandMatrix.cc:68
gbl::GblTrajectory::externalSeed
TMatrixDSym externalSeed
Precision (inverse covariance matrix) of external seed.
Definition: GblTrajectory.h:104
gbl::GblTrajectory::numMeasurements
unsigned int numMeasurements
Total number of measurements.
Definition: GblTrajectory.h:95
gbl::MilleBinary::writeRecord
void writeRecord()
Write record to file.
Definition: MilleBinary.cc:111
gbl::BorderedBandMatrix::solveAndInvertBorderedBand
void solveAndInvertBorderedBand(const VVector &aRightHandSide, VVector &aSolution)
Solve linear equation system, partially calculate inverse.
Definition: BorderedBandMatrix.cc:147
gbl::GblTrajectory::getNumPoints
unsigned int getNumPoints() const
Retrieve number of point from trajectory.
Definition: GblTrajectory.cc:294
gbl::GblTrajectory::theDimension
std::vector< unsigned int > theDimension
List of active dimensions (0=u1, 1=u2) in fit.
Definition: GblTrajectory.h:99
gbl::GblTrajectory::numInnerTrans
unsigned int numInnerTrans
Number of inner transformations to external parameters.
Definition: GblTrajectory.h:91
gbl::MilleBinary::addData
void addData(double aMeas, double aPrec, const std::vector< unsigned int > &indLocal, const std::vector< double > &derLocal, const std::vector< int > &labGlobal, const std::vector< double > &derGlobal)
Add data block to (end of) record.
Definition: MilleBinary.cc:70
gbl::GblTrajectory::predict
void predict()
Calculate predictions for all points.
Definition: GblTrajectory.cc:1011
gbl::GblTrajectory::prepare
void prepare()
Prepare fit for simple or composed trajectory.
Definition: GblTrajectory.cc:797
gbl::GblPoint::getOffset
int getOffset() const
Retrieve offset for point.
Definition: GblPoint.cc:371
gbl::GblTrajectory::theVector
VVector theVector
Vector of linear equation system.
Definition: GblTrajectory.h:110
gbl::BorderedBandMatrix::resize
void resize(unsigned int nSize, unsigned int nBorder=1, unsigned int nBand=5)
Resize bordered band matrix.
Definition: BorderedBandMatrix.cc:48
gbl::GblTrajectory::getLabels
unsigned int getLabels(std::vector< unsigned int > &aLabelList)
Get (list of) labels of points on (simple) valid trajectory.
Definition: GblTrajectory.cc:711
gbl::GblTrajectory::isValid
bool isValid() const
Retrieve validity of trajectory.
Definition: GblTrajectory.cc:289
gbl::BorderedBandMatrix::getBlockMatrix
TMatrixDSym getBlockMatrix(const std::vector< unsigned int > anIndex) const
Retrieve symmetric block matrix.
Definition: BorderedBandMatrix.cc:97
gbl::GblTrajectory::numCurvature
unsigned int numCurvature
Number of curvature parameters (0 or 1) or external parameters.
Definition: GblTrajectory.h:92
GblTrajectory.h
gbl::GblTrajectory::buildLinearEquationSystem
void buildLinearEquationSystem()
Build linear equation system from data (blocks).
Definition: GblTrajectory.cc:776
gbl::GblTrajectory::innerTransformations
std::vector< TMatrixD > innerTransformations
Transformations at innermost points of.
Definition: GblTrajectory.h:105
gbl::GblTrajectory::theData
std::vector< GblData > theData
List of data blocks.
Definition: GblTrajectory.h:101
gbl::GblTrajectory::numParameters
unsigned int numParameters
Number of fit parameters.
Definition: GblTrajectory.h:93
gbl::GblTrajectory::GblTrajectory
GblTrajectory(const std::vector< GblPoint > &aPointList, bool flagCurv=true, bool flagU1dir=true, bool flagU2dir=true)
Create new (simple) trajectory from list of points.
Definition: GblTrajectory.cc:147
gbl::GblTrajectory::constructOK
bool constructOK
Trajectory has been successfully constructed (ready for fit/output)
Definition: GblTrajectory.h:97
gbl::MilleBinary
Millepede-II (binary) record.
Definition: MilleBinary.h:68
gbl::GblPoint
Point on trajectory.
Definition: GblPoint.h:68
gbl::GblTrajectory::calcJacobians
void calcJacobians()
Calculate Jacobians to previous/next scatterer from point to point ones.
Definition: GblTrajectory.cc:366
gbl::GblTrajectory::getFitToKinkJacobian
void getFitToKinkJacobian(std::vector< unsigned int > &anIndex, SMatrix27 &aJacobian, const GblPoint &aPoint) const
Get jacobian for transformation from (trajectory) fit to kink parameters at point.
Definition: GblTrajectory.cc:583
gbl::GblTrajectory::defineOffsets
void defineOffsets()
Define offsets from list of points.
Definition: GblTrajectory.cc:344
gbl::GblTrajectory::fit
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight, std::string optionList="")
Perform fit of (valid) trajectory.
Definition: GblTrajectory.cc:1043
gbl::GblTrajectory::numAllPoints
unsigned int numAllPoints
Number of all points on trajectory.
Definition: GblTrajectory.h:87