GENFIT  Rev:NoNumberAvailable
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 
102 #include "GblTrajectory.h"
103 
105 namespace gbl {
106 
108 
116 GblTrajectory::GblTrajectory(const std::vector<GblPoint> &aPointList,
117  bool flagCurv, bool flagU1dir, bool flagU2dir) :
118  numAllPoints(aPointList.size()), numPoints(), numOffsets(0), numInnerTrans(
119  0), numCurvature(flagCurv ? 1 : 0), numParameters(0), numLocals(
120  0), numMeasurements(0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
121 
122  if (flagU1dir)
123  theDimension.push_back(0);
124  if (flagU2dir)
125  theDimension.push_back(1);
126  // simple (single) trajectory
127  thePoints.push_back(aPointList);
128  numPoints.push_back(numAllPoints);
129  construct(); // construct trajectory
130 }
131 
133 
144 GblTrajectory::GblTrajectory(const std::vector<GblPoint> &aPointList,
145  unsigned int aLabel, const TMatrixDSym &aSeed, bool flagCurv,
146  bool flagU1dir, bool flagU2dir) :
147  numAllPoints(aPointList.size()), numPoints(), numOffsets(0), numInnerTrans(
148  0), numCurvature(flagCurv ? 1 : 0), numParameters(0), numLocals(
149  0), numMeasurements(0), externalPoint(aLabel), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(
150  aSeed), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
151 
152  if (flagU1dir)
153  theDimension.push_back(0);
154  if (flagU2dir)
155  theDimension.push_back(1);
156  // simple (single) trajectory
157  thePoints.push_back(aPointList);
158  numPoints.push_back(numAllPoints);
159  construct(); // construct trajectory
160 }
161 
163 
168  const std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > &aPointsAndTransList) :
169  numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
170  aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
171  0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
172 
173  for (unsigned int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
174  thePoints.push_back(aPointsAndTransList[iTraj].first);
175  numPoints.push_back(thePoints.back().size());
176  numAllPoints += numPoints.back();
177  innerTransformations.push_back(aPointsAndTransList[iTraj].second);
178  }
179  theDimension.push_back(0);
180  theDimension.push_back(1);
181  numCurvature = innerTransformations[0].GetNcols();
182  construct(); // construct (composed) trajectory
183 }
184 
186 
194  const std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > &aPointsAndTransList,
195  const TMatrixD &extDerivatives, const TVectorD &extMeasurements,
196  const TVectorD &extPrecisions) :
197  numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
198  aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
199  0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(
200  extDerivatives), externalMeasurements(extMeasurements), externalPrecisions(
201  extPrecisions) {
202 
203  for (unsigned int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
204  thePoints.push_back(aPointsAndTransList[iTraj].first);
205  numPoints.push_back(thePoints.back().size());
206  numAllPoints += numPoints.back();
207  innerTransformations.push_back(aPointsAndTransList[iTraj].second);
208  }
209  theDimension.push_back(0);
210  theDimension.push_back(1);
211  numCurvature = innerTransformations[0].GetNcols();
212  construct(); // construct (composed) trajectory
213 }
214 
216 }
217 
220  return constructOK;
221 }
222 
224 unsigned int GblTrajectory::getNumPoints() const {
225  return numAllPoints;
226 }
227 
229 
233 
234  constructOK = false;
235  fitOK = false;
236  unsigned int aLabel = 0;
237  if (numAllPoints < 2) {
238  std::cout << " GblTrajectory construction failed: too few GblPoints "
239  << std::endl;
240  return;
241  }
242  // loop over trajectories
243  numTrajectories = thePoints.size();
244  //std::cout << " numTrajectories: " << numTrajectories << ", " << innerTransformations.size() << std::endl;
245  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
246  std::vector<GblPoint>::iterator itPoint;
247  for (itPoint = thePoints[iTraj].begin();
248  itPoint < thePoints[iTraj].end(); ++itPoint) {
249  numLocals = std::max(numLocals, itPoint->getNumLocals());
250  numMeasurements += itPoint->hasMeasurement();
251  itPoint->setLabel(++aLabel);
252  }
253  }
254  defineOffsets();
255  calcJacobians();
256  try {
257  prepare();
258  } catch (std::overflow_error &e) {
259  std::cout << " GblTrajectory construction failed: " << e.what()
260  << std::endl;
261  return;
262  }
263  constructOK = true;
264  // number of fit parameters
267 }
268 
270 
275 
276  // loop over trajectories
277  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
278  // first point is offset
279  thePoints[iTraj].front().setOffset(numOffsets++);
280  // intermediate scatterers are offsets
281  std::vector<GblPoint>::iterator itPoint;
282  for (itPoint = thePoints[iTraj].begin() + 1;
283  itPoint < thePoints[iTraj].end() - 1; ++itPoint) {
284  if (itPoint->hasScatterer()) {
285  itPoint->setOffset(numOffsets++);
286  } else {
287  itPoint->setOffset(-numOffsets);
288  }
289  }
290  // last point is offset
291  thePoints[iTraj].back().setOffset(numOffsets++);
292  }
293 }
294 
297 
298  SMatrix55 scatJacobian;
299  // loop over trajectories
300  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
301  // forward propagation (all)
302  GblPoint* previousPoint = &thePoints[iTraj].front();
303  unsigned int numStep = 0;
304  std::vector<GblPoint>::iterator itPoint;
305  for (itPoint = thePoints[iTraj].begin() + 1;
306  itPoint < thePoints[iTraj].end(); ++itPoint) {
307  if (numStep == 0) {
308  scatJacobian = itPoint->getP2pJacobian();
309  } else {
310  scatJacobian = itPoint->getP2pJacobian() * scatJacobian;
311  }
312  numStep++;
313  itPoint->addPrevJacobian(scatJacobian); // iPoint -> previous scatterer
314  if (itPoint->getOffset() >= 0) {
315  previousPoint->addNextJacobian(scatJacobian); // lastPoint -> next scatterer
316  numStep = 0;
317  previousPoint = &(*itPoint);
318  }
319  }
320  // backward propagation (without scatterers)
321  for (itPoint = thePoints[iTraj].end() - 1;
322  itPoint > thePoints[iTraj].begin(); --itPoint) {
323  if (itPoint->getOffset() >= 0) {
324  scatJacobian = itPoint->getP2pJacobian();
325  continue; // skip offsets
326  }
327  itPoint->addNextJacobian(scatJacobian); // iPoint -> next scatterer
328  scatJacobian = scatJacobian * itPoint->getP2pJacobian();
329  }
330  }
331 }
332 
334 
342 std::pair<std::vector<unsigned int>, TMatrixD> GblTrajectory::getJacobian(
343  int aSignedLabel) const {
344 
345  unsigned int nDim = theDimension.size();
346  unsigned int nCurv = numCurvature;
347  unsigned int nLocals = numLocals;
348  unsigned int nBorder = nCurv + nLocals;
349  unsigned int nParBRL = nBorder + 2 * nDim;
350  unsigned int nParLoc = nLocals + 5;
351  std::vector<unsigned int> anIndex;
352  anIndex.reserve(nParBRL);
353  TMatrixD aJacobian(nParLoc, nParBRL);
354  aJacobian.Zero();
355 
356  unsigned int aLabel = abs(aSignedLabel);
357  unsigned int firstLabel = 1;
358  unsigned int lastLabel = 0;
359  unsigned int aTrajectory = 0;
360  // loop over trajectories
361  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
362  aTrajectory = iTraj;
363  lastLabel += numPoints[iTraj];
364  if (aLabel <= lastLabel)
365  break;
366  if (iTraj < numTrajectories - 1)
367  firstLabel += numPoints[iTraj];
368  }
369  int nJacobian; // 0: prev, 1: next
370  // check consistency of (index, direction)
371  if (aSignedLabel > 0) {
372  nJacobian = 1;
373  if (aLabel >= lastLabel) {
374  aLabel = lastLabel;
375  nJacobian = 0;
376  }
377  } else {
378  nJacobian = 0;
379  if (aLabel <= firstLabel) {
380  aLabel = firstLabel;
381  nJacobian = 1;
382  }
383  }
384  const GblPoint aPoint = thePoints[aTrajectory][aLabel - firstLabel];
385  std::vector<unsigned int> labDer(5);
386  SMatrix55 matDer;
387  getFitToLocalJacobian(labDer, matDer, aPoint, 5, nJacobian);
388 
389  // from local parameters
390  for (unsigned int i = 0; i < nLocals; ++i) {
391  aJacobian(i + 5, i) = 1.0;
392  anIndex.push_back(i + 1);
393  }
394  // from trajectory parameters
395  unsigned int iCol = nLocals;
396  for (unsigned int i = 0; i < 5; ++i) {
397  if (labDer[i] > 0) {
398  anIndex.push_back(labDer[i]);
399  for (unsigned int j = 0; j < 5; ++j) {
400  aJacobian(j, iCol) = matDer(j, i);
401  }
402  ++iCol;
403  }
404  }
405  return std::make_pair(anIndex, aJacobian);
406 }
407 
409 
418 void GblTrajectory::getFitToLocalJacobian(std::vector<unsigned int> &anIndex,
419  SMatrix55 &aJacobian, const GblPoint &aPoint, unsigned int measDim,
420  unsigned int nJacobian) const {
421 
422  unsigned int nDim = theDimension.size();
423  unsigned int nCurv = numCurvature;
424  unsigned int nLocals = numLocals;
425 
426  int nOffset = aPoint.getOffset();
427 
428  if (nOffset < 0) // need interpolation
429  {
430  SMatrix22 prevW, prevWJ, nextW, nextWJ, matN;
431  SVector2 prevWd, nextWd;
432  int ierr;
433  aPoint.getDerivatives(0, prevW, prevWJ, prevWd); // W-, W- * J-, W- * d-
434  aPoint.getDerivatives(1, nextW, nextWJ, nextWd); // W-, W- * J-, W- * d-
435  const SMatrix22 sumWJ(prevWJ + nextWJ);
436  matN = sumWJ.Inverse(ierr); // N = (W- * J- + W+ * J+)^-1
437  // derivatives for u_int
438  const SMatrix22 prevNW(matN * prevW); // N * W-
439  const SMatrix22 nextNW(matN * nextW); // N * W+
440  const SVector2 prevNd(matN * prevWd); // N * W- * d-
441  const SVector2 nextNd(matN * nextWd); // N * W+ * d+
442 
443  unsigned int iOff = nDim * (-nOffset - 1) + nLocals + nCurv + 1; // first offset ('i' in u_i)
444 
445  // local offset
446  if (nCurv > 0) {
447  aJacobian.Place_in_col(-prevNd - nextNd, 3, 0); // from curvature
448  anIndex[0] = nLocals + 1;
449  }
450  aJacobian.Place_at(prevNW, 3, 1); // from 1st Offset
451  aJacobian.Place_at(nextNW, 3, 3); // from 2nd Offset
452  for (unsigned int i = 0; i < nDim; ++i) {
453  anIndex[1 + theDimension[i]] = iOff + i;
454  anIndex[3 + theDimension[i]] = iOff + nDim + i;
455  }
456 
457  // local slope and curvature
458  if (measDim > 2) {
459  // derivatives for u'_int
460  const SMatrix22 prevWPN(nextWJ * prevNW); // W+ * J+ * N * W-
461  const SMatrix22 nextWPN(prevWJ * nextNW); // W- * J- * N * W+
462  const SVector2 prevWNd(nextWJ * prevNd); // W+ * J+ * N * W- * d-
463  const SVector2 nextWNd(prevWJ * nextNd); // W- * J- * N * W+ * d+
464  if (nCurv > 0) {
465  aJacobian(0, 0) = 1.0;
466  aJacobian.Place_in_col(prevWNd - nextWNd, 1, 0); // from curvature
467  }
468  aJacobian.Place_at(-prevWPN, 1, 1); // from 1st Offset
469  aJacobian.Place_at(nextWPN, 1, 3); // from 2nd Offset
470  }
471  } else { // at point
472  // anIndex must be sorted
473  // forward : iOff2 = iOff1 + nDim, index1 = 1, index2 = 3
474  // backward: iOff2 = iOff1 - nDim, index1 = 3, index2 = 1
475  unsigned int iOff1 = nDim * nOffset + nCurv + nLocals + 1; // first offset ('i' in u_i)
476  unsigned int index1 = 3 - 2 * nJacobian; // index of first offset
477  unsigned int iOff2 = iOff1 + nDim * (nJacobian * 2 - 1); // second offset ('i' in u_i)
478  unsigned int index2 = 1 + 2 * nJacobian; // index of second offset
479  // local offset
480  aJacobian(3, index1) = 1.0; // from 1st Offset
481  aJacobian(4, index1 + 1) = 1.0;
482  for (unsigned int i = 0; i < nDim; ++i) {
483  anIndex[index1 + theDimension[i]] = iOff1 + i;
484  }
485 
486  // local slope and curvature
487  if (measDim > 2) {
488  SMatrix22 matW, matWJ;
489  SVector2 vecWd;
490  aPoint.getDerivatives(nJacobian, matW, matWJ, vecWd); // W, W * J, W * d
491  double sign = (nJacobian > 0) ? 1. : -1.;
492  if (nCurv > 0) {
493  aJacobian(0, 0) = 1.0;
494  aJacobian.Place_in_col(-sign * vecWd, 1, 0); // from curvature
495  anIndex[0] = nLocals + 1;
496  }
497  aJacobian.Place_at(-sign * matWJ, 1, index1); // from 1st Offset
498  aJacobian.Place_at(sign * matW, 1, index2); // from 2nd Offset
499  for (unsigned int i = 0; i < nDim; ++i) {
500  anIndex[index2 + theDimension[i]] = iOff2 + i;
501  }
502  }
503  }
504 }
505 
507 
513 void GblTrajectory::getFitToKinkJacobian(std::vector<unsigned int> &anIndex,
514  SMatrix27 &aJacobian, const GblPoint &aPoint) const {
515 
516  unsigned int nDim = theDimension.size();
517  unsigned int nCurv = numCurvature;
518  unsigned int nLocals = numLocals;
519 
520  int nOffset = aPoint.getOffset();
521 
522  SMatrix22 prevW, prevWJ, nextW, nextWJ;
523  SVector2 prevWd, nextWd;
524  aPoint.getDerivatives(0, prevW, prevWJ, prevWd); // W-, W- * J-, W- * d-
525  aPoint.getDerivatives(1, nextW, nextWJ, nextWd); // W-, W- * J-, W- * d-
526  const SMatrix22 sumWJ(prevWJ + nextWJ); // W- * J- + W+ * J+
527  const SVector2 sumWd(prevWd + nextWd); // W+ * d+ + W- * d-
528 
529  unsigned int iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1; // first offset ('i' in u_i)
530 
531  // local offset
532  if (nCurv > 0) {
533  aJacobian.Place_in_col(-sumWd, 0, 0); // from curvature
534  anIndex[0] = nLocals + 1;
535  }
536  aJacobian.Place_at(prevW, 0, 1); // from 1st Offset
537  aJacobian.Place_at(-sumWJ, 0, 3); // from 2nd Offset
538  aJacobian.Place_at(nextW, 0, 5); // from 1st Offset
539  for (unsigned int i = 0; i < nDim; ++i) {
540  anIndex[1 + theDimension[i]] = iOff + i;
541  anIndex[3 + theDimension[i]] = iOff + nDim + i;
542  anIndex[5 + theDimension[i]] = iOff + nDim * 2 + i;
543  }
544 }
545 
547 
556 unsigned int GblTrajectory::getResults(int aSignedLabel, TVectorD &localPar,
557  TMatrixDSym &localCov) const {
558  if (not fitOK)
559  return 1;
560  std::pair<std::vector<unsigned int>, TMatrixD> indexAndJacobian =
561  getJacobian(aSignedLabel);
562  unsigned int nParBrl = indexAndJacobian.first.size();
563  TVectorD aVec(nParBrl); // compressed vector
564  for (unsigned int i = 0; i < nParBrl; ++i) {
565  aVec[i] = theVector(indexAndJacobian.first[i] - 1);
566  }
567  TMatrixDSym aMat = theMatrix.getBlockMatrix(indexAndJacobian.first); // compressed matrix
568  localPar = indexAndJacobian.second * aVec;
569  localCov = aMat.Similarity(indexAndJacobian.second);
570  return 0;
571 }
572 
574 
586 unsigned int GblTrajectory::getMeasResults(unsigned int aLabel,
587  unsigned int &numData, TVectorD &aResiduals, TVectorD &aMeasErrors,
588  TVectorD &aResErrors, TVectorD &aDownWeights) {
589  numData = 0;
590  if (not fitOK)
591  return 1;
592 
593  unsigned int firstData = measDataIndex[aLabel - 1]; // first data block with measurement
594  numData = measDataIndex[aLabel] - firstData; // number of data blocks
595  for (unsigned int i = 0; i < numData; ++i) {
596  getResAndErr(firstData + i, aResiduals[i], aMeasErrors[i],
597  aResErrors[i], aDownWeights[i]);
598  }
599  return 0;
600 }
601 
603 
615 unsigned int GblTrajectory::getScatResults(unsigned int aLabel,
616  unsigned int &numData, TVectorD &aResiduals, TVectorD &aMeasErrors,
617  TVectorD &aResErrors, TVectorD &aDownWeights) {
618  numData = 0;
619  if (not fitOK)
620  return 1;
621 
622  unsigned int firstData = scatDataIndex[aLabel - 1]; // first data block with scatterer
623  numData = scatDataIndex[aLabel] - firstData; // number of data blocks
624  for (unsigned int i = 0; i < numData; ++i) {
625  getResAndErr(firstData + i, aResiduals[i], aMeasErrors[i],
626  aResErrors[i], aDownWeights[i]);
627  }
628  return 0;
629 }
630 
632 
635 void GblTrajectory::getLabels(std::vector<unsigned int> &aLabelList) {
636  unsigned int aLabel = 0;
637  unsigned int nPoint = thePoints[0].size();
638  aLabelList.resize(nPoint);
639  for (unsigned i = 0; i < nPoint; ++i) {
640  aLabelList[i] = ++aLabel;
641  }
642 }
643 
645 
649  std::vector<std::vector<unsigned int> > &aLabelList) {
650  unsigned int aLabel = 0;
651  aLabelList.resize(numTrajectories);
652  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
653  unsigned int nPoint = thePoints[iTraj].size();
654  aLabelList[iTraj].resize(nPoint);
655  for (unsigned i = 0; i < nPoint; ++i) {
656  aLabelList[iTraj][i] = ++aLabel;
657  }
658  }
659 }
660 
662 
671 void GblTrajectory::getResAndErr(unsigned int aData, double &aResidual,
672  double &aMeasError, double &aResError, double &aDownWeight) {
673 
674  double aMeasVar;
675  std::vector<unsigned int>* indLocal;
676  std::vector<double>* derLocal;
677  theData[aData].getResidual(aResidual, aMeasVar, aDownWeight, indLocal,
678  derLocal);
679  unsigned int nParBrl = (*indLocal).size();
680  TVectorD aVec(nParBrl); // compressed vector of derivatives
681  for (unsigned int j = 0; j < nParBrl; ++j) {
682  aVec[j] = (*derLocal)[j];
683  }
684  TMatrixDSym aMat = theMatrix.getBlockMatrix(*indLocal); // compressed (covariance) matrix
685  double aFitVar = aMat.Similarity(aVec); // variance from track fit
686  aMeasError = sqrt(aMeasVar); // error of measurement
687  aResError = (aFitVar < aMeasVar ? sqrt(aMeasVar - aFitVar) : 0.); // error of residual
688 }
689 
692  unsigned int nBorder = numCurvature + numLocals;
694  theMatrix.resize(numParameters, nBorder);
695  double aValue, aWeight;
696  std::vector<unsigned int>* indLocal;
697  std::vector<double>* derLocal;
698  std::vector<GblData>::iterator itData;
699  for (itData = theData.begin(); itData < theData.end(); ++itData) {
700  itData->getLocalData(aValue, aWeight, indLocal, derLocal);
701  for (unsigned int j = 0; j < indLocal->size(); ++j) {
702  theVector((*indLocal)[j] - 1) += (*derLocal)[j] * aWeight * aValue;
703  }
704  theMatrix.addBlockMatrix(aWeight, indLocal, derLocal);
705  }
706 }
707 
709 
713  unsigned int nDim = theDimension.size();
714  // upper limit
715  unsigned int maxData = numMeasurements + nDim * (numOffsets - 2)
716  + externalSeed.GetNrows();
717  theData.reserve(maxData);
718  measDataIndex.resize(numAllPoints + 3); // include external seed and measurements
719  scatDataIndex.resize(numAllPoints + 1);
720  unsigned int nData = 0;
721  std::vector<TMatrixD> innerTransDer;
722  std::vector<std::vector<unsigned int> > innerTransLab;
723  // composed trajectory ?
724  if (numInnerTrans > 0) {
725  //std::cout << "composed trajectory" << std::endl;
726  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
727  // innermost point
728  GblPoint* innerPoint = &thePoints[iTraj].front();
729  // transformation fit to local track parameters
730  std::vector<unsigned int> firstLabels(5);
731  SMatrix55 matFitToLocal, matLocalToFit;
732  getFitToLocalJacobian(firstLabels, matFitToLocal, *innerPoint, 5);
733  // transformation local track to fit parameters
734  int ierr;
735  matLocalToFit = matFitToLocal.Inverse(ierr);
736  TMatrixD localToFit(5, 5);
737  for (unsigned int i = 0; i < 5; ++i) {
738  for (unsigned int j = 0; j < 5; ++j) {
739  localToFit(i, j) = matLocalToFit(i, j);
740  }
741  }
742  // transformation external to fit parameters at inner (first) point
743  innerTransDer.push_back(localToFit * innerTransformations[iTraj]);
744  innerTransLab.push_back(firstLabels);
745  }
746  }
747  // measurements
748  SMatrix55 matP;
749  // loop over trajectories
750  std::vector<GblPoint>::iterator itPoint;
751  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
752  for (itPoint = thePoints[iTraj].begin();
753  itPoint < thePoints[iTraj].end(); ++itPoint) {
754  SVector5 aMeas, aPrec;
755  unsigned int nLabel = itPoint->getLabel();
756  unsigned int measDim = itPoint->hasMeasurement();
757  if (measDim) {
758  const TMatrixD localDer = itPoint->getLocalDerivatives();
759  const std::vector<int> globalLab = itPoint->getGlobalLabels();
760  const TMatrixD globalDer = itPoint->getGlobalDerivatives();
761  TMatrixD transDer;
762  itPoint->getMeasurement(matP, aMeas, aPrec);
763  unsigned int iOff = 5 - measDim; // first active component
764  std::vector<unsigned int> labDer(5);
765  SMatrix55 matDer, matPDer;
766  unsigned int nJacobian =
767  (itPoint < thePoints[iTraj].end() - 1) ? 1 : 0; // last point needs backward propagation
768  getFitToLocalJacobian(labDer, matDer, *itPoint, measDim,
769  nJacobian);
770  if (measDim > 2) {
771  matPDer = matP * matDer;
772  } else { // 'shortcut' for position measurements
773  matPDer.Place_at(
774  matP.Sub<SMatrix22>(3, 3)
775  * matDer.Sub<SMatrix25>(3, 0), 3, 0);
776  }
777 
778  if (numInnerTrans > 0) {
779  // transform for external parameters
780  TMatrixD proDer(measDim, 5);
781  // match parameters
782  unsigned int ifirst = 0;
783  unsigned int ilabel = 0;
784  while (ilabel < 5) {
785  if (labDer[ilabel] > 0) {
786  while (innerTransLab[iTraj][ifirst]
787  != labDer[ilabel] and ifirst < 5) {
788  ++ifirst;
789  }
790  if (ifirst >= 5) {
791  labDer[ilabel] -= 2 * nDim * (iTraj + 1); // adjust label
792  } else {
793  // match
794  labDer[ilabel] = 0; // mark as related to external parameters
795  for (unsigned int k = iOff; k < 5; ++k) {
796  proDer(k - iOff, ifirst) = matPDer(k,
797  ilabel);
798  }
799  }
800  }
801  ++ilabel;
802  }
803  transDer.ResizeTo(measDim, numCurvature);
804  transDer = proDer * innerTransDer[iTraj];
805  }
806  for (unsigned int i = iOff; i < 5; ++i) {
807  if (aPrec(i) > 0.) {
808  GblData aData(nLabel, aMeas(i), aPrec(i));
809  aData.addDerivatives(i, labDer, matPDer, iOff, localDer,
810  globalLab, globalDer, numLocals, transDer);
811  theData.push_back(aData);
812  nData++;
813  }
814  }
815 
816  }
817  measDataIndex[nLabel] = nData;
818  }
819  }
820 
821  // pseudo measurements from kinks
822  SMatrix22 matT;
823  scatDataIndex[0] = nData;
824  scatDataIndex[1] = nData;
825  // loop over trajectories
826  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
827  for (itPoint = thePoints[iTraj].begin() + 1;
828  itPoint < thePoints[iTraj].end() - 1; ++itPoint) {
829  SVector2 aMeas, aPrec;
830  unsigned int nLabel = itPoint->getLabel();
831  if (itPoint->hasScatterer()) {
832  itPoint->getScatterer(matT, aMeas, aPrec);
833  TMatrixD transDer;
834  std::vector<unsigned int> labDer(7);
835  SMatrix27 matDer, matTDer;
836  getFitToKinkJacobian(labDer, matDer, *itPoint);
837  matTDer = matT * matDer;
838  if (numInnerTrans > 0) {
839  // transform for external parameters
840  TMatrixD proDer(nDim, 5);
841  // match parameters
842  unsigned int ifirst = 0;
843  unsigned int ilabel = 0;
844  while (ilabel < 7) {
845  if (labDer[ilabel] > 0) {
846  while (innerTransLab[iTraj][ifirst]
847  != labDer[ilabel] and ifirst < 5) {
848  ++ifirst;
849  }
850  if (ifirst >= 5) {
851  labDer[ilabel] -= 2 * nDim * (iTraj + 1); // adjust label
852  } else {
853  // match
854  labDer[ilabel] = 0; // mark as related to external parameters
855  for (unsigned int k = 0; k < nDim; ++k) {
856  proDer(k, ifirst) = matTDer(k, ilabel);
857  }
858  }
859  }
860  ++ilabel;
861  }
862  transDer.ResizeTo(nDim, numCurvature);
863  transDer = proDer * innerTransDer[iTraj];
864  }
865  for (unsigned int i = 0; i < nDim; ++i) {
866  unsigned int iDim = theDimension[i];
867  if (aPrec(iDim) > 0.) {
868  GblData aData(nLabel, aMeas(iDim), aPrec(iDim));
869  aData.addDerivatives(iDim, labDer, matTDer, numLocals,
870  transDer);
871  theData.push_back(aData);
872  nData++;
873  }
874  }
875  }
876  scatDataIndex[nLabel] = nData;
877  }
878  scatDataIndex[thePoints[iTraj].back().getLabel()] = nData;
879  }
880 
881  // external seed
882  if (externalPoint > 0) {
883  std::pair<std::vector<unsigned int>, TMatrixD> indexAndJacobian =
885  std::vector<unsigned int> externalSeedIndex = indexAndJacobian.first;
886  std::vector<double> externalSeedDerivatives(externalSeedIndex.size());
887  const TMatrixDSymEigen externalSeedEigen(externalSeed);
888  const TVectorD valEigen = externalSeedEigen.GetEigenValues();
889  TMatrixD vecEigen = externalSeedEigen.GetEigenVectors();
890  vecEigen = vecEigen.T() * indexAndJacobian.second;
891  for (int i = 0; i < externalSeed.GetNrows(); ++i) {
892  if (valEigen(i) > 0.) {
893  for (int j = 0; j < externalSeed.GetNcols(); ++j) {
894  externalSeedDerivatives[j] = vecEigen(i, j);
895  }
896  GblData aData(externalPoint, 0., valEigen(i));
897  aData.addDerivatives(externalSeedIndex, externalSeedDerivatives);
898  theData.push_back(aData);
899  nData++;
900  }
901  }
902  }
903  measDataIndex[numAllPoints + 1] = nData;
904  // external measurements
905  unsigned int nExt = externalMeasurements.GetNrows();
906  if (nExt > 0) {
907  std::vector<unsigned int> index(numCurvature);
908  std::vector<double> derivatives(numCurvature);
909  for (unsigned int iExt = 0; iExt < nExt; ++iExt) {
910  for (unsigned int iCol = 0; iCol < numCurvature; ++iCol) {
911  index[iCol] = iCol + 1;
912  derivatives[iCol] = externalDerivatives(iExt, iCol);
913  }
914  GblData aData(1U, externalMeasurements(iExt),
915  externalPrecisions(iExt));
916  aData.addDerivatives(index, derivatives);
917  theData.push_back(aData);
918  nData++;
919  }
920  }
921  measDataIndex[numAllPoints + 2] = nData;
922 }
923 
926  std::vector<GblData>::iterator itData;
927  for (itData = theData.begin(); itData < theData.end(); ++itData) {
928  itData->setPrediction(theVector);
929  }
930 }
931 
933 
936 double GblTrajectory::downWeight(unsigned int aMethod) {
937  double aLoss = 0.;
938  std::vector<GblData>::iterator itData;
939  for (itData = theData.begin(); itData < theData.end(); ++itData) {
940  aLoss += (1. - itData->setDownWeighting(aMethod));
941  }
942  return aLoss;
943 }
944 
946 
955 unsigned int GblTrajectory::fit(double &Chi2, int &Ndf, double &lostWeight,
956  std::string optionList) {
957  const double normChi2[4] = { 1.0, 0.8737, 0.9326, 0.8228 };
958  const std::string methodList = "TtHhCc";
959 
960  Chi2 = 0.;
961  Ndf = -1;
962  lostWeight = 0.;
963  if (not constructOK)
964  return 10;
965 
966  unsigned int aMethod = 0;
967 
969  lostWeight = 0.;
970  unsigned int ierr = 0;
971  try {
972 
974  predict();
975 
976  for (unsigned int i = 0; i < optionList.size(); ++i) // down weighting iterations
977  {
978  size_t aPosition = methodList.find(optionList[i]);
979  if (aPosition != std::string::npos) {
980  aMethod = aPosition / 2 + 1;
981  lostWeight = downWeight(aMethod);
984  predict();
985  }
986  }
987  Ndf = theData.size() - numParameters;
988  Chi2 = 0.;
989  for (unsigned int i = 0; i < theData.size(); ++i) {
990  Chi2 += theData[i].getChi2();
991  }
992  Chi2 /= normChi2[aMethod];
993  fitOK = true;
994 
995  } catch (int e) {
996  std::cout << " GblTrajectory::fit exception " << e << std::endl;
997  ierr = e;
998  }
999  return ierr;
1000 }
1001 
1004  double aValue;
1005  double aErr;
1006  std::vector<unsigned int>* indLocal;
1007  std::vector<double>* derLocal;
1008  std::vector<int>* labGlobal;
1009  std::vector<double>* derGlobal;
1010 
1011  if (not constructOK)
1012  return;
1013 
1014 // data: measurements, kinks and external seed
1015  std::vector<GblData>::iterator itData;
1016  for (itData = theData.begin(); itData != theData.end(); ++itData) {
1017  itData->getAllData(aValue, aErr, indLocal, derLocal, labGlobal,
1018  derGlobal);
1019  aMille.addData(aValue, aErr, *indLocal, *derLocal, *labGlobal,
1020  *derGlobal);
1021  }
1022  aMille.writeRecord();
1023 }
1024 
1026 
1029 void GblTrajectory::printTrajectory(unsigned int level) {
1030  if (numInnerTrans) {
1031  std::cout << "Composed GblTrajectory, " << numInnerTrans
1032  << " subtrajectories" << std::endl;
1033  } else {
1034  std::cout << "Simple GblTrajectory" << std::endl;
1035  }
1036  if (theDimension.size() < 2) {
1037  std::cout << " 2D-trajectory" << std::endl;
1038  }
1039  std::cout << " Number of GblPoints : " << numAllPoints
1040  << std::endl;
1041  std::cout << " Number of points with offsets: " << numOffsets << std::endl;
1042  std::cout << " Number of fit parameters : " << numParameters
1043  << std::endl;
1044  std::cout << " Number of measurements : " << numMeasurements
1045  << std::endl;
1046  if (externalMeasurements.GetNrows()) {
1047  std::cout << " Number of ext. measurements : "
1048  << externalMeasurements.GetNrows() << std::endl;
1049  }
1050  if (externalPoint) {
1051  std::cout << " Label of point with ext. seed: " << externalPoint
1052  << std::endl;
1053  }
1054  if (constructOK) {
1055  std::cout << " Constructed OK " << std::endl;
1056  }
1057  if (fitOK) {
1058  std::cout << " Fitted OK " << std::endl;
1059  }
1060  if (level > 0) {
1061  if (numInnerTrans) {
1062  std::cout << " Inner transformations" << std::endl;
1063  for (unsigned int i = 0; i < numInnerTrans; ++i) {
1064  innerTransformations[i].Print();
1065  }
1066  }
1067  if (externalMeasurements.GetNrows()) {
1068  std::cout << " External measurements" << std::endl;
1069  std::cout << " Measurements:" << std::endl;
1070  externalMeasurements.Print();
1071  std::cout << " Precisions:" << std::endl;
1072  externalPrecisions.Print();
1073  std::cout << " Derivatives:" << std::endl;
1074  externalDerivatives.Print();
1075  }
1076  if (externalPoint) {
1077  std::cout << " External seed:" << std::endl;
1078  externalSeed.Print();
1079  }
1080  if (fitOK) {
1081  std::cout << " Fit results" << std::endl;
1082  std::cout << " Parameters:" << std::endl;
1083  theVector.print();
1084  std::cout << " Covariance matrix (bordered band part):"
1085  << std::endl;
1087  }
1088  }
1089 }
1090 
1092 
1095 void GblTrajectory::printPoints(unsigned int level) {
1096  std::cout << "GblPoints " << std::endl;
1097  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
1098  std::vector<GblPoint>::iterator itPoint;
1099  for (itPoint = thePoints[iTraj].begin();
1100  itPoint < thePoints[iTraj].end(); ++itPoint) {
1101  itPoint->printPoint(level);
1102  }
1103  }
1104 }
1105 
1108  std::cout << "GblData blocks " << std::endl;
1109  std::vector<GblData>::iterator itData;
1110  for (itData = theData.begin(); itData < theData.end(); ++itData) {
1111  itData->printData();
1112  }
1113 }
1114 
1115 }
gbl::GblTrajectory::externalPoint
unsigned int externalPoint
Label of external point (or 0)
Definition: GblTrajectory.h:69
gbl::VVector::resize
void resize(const unsigned int nRows)
Resize vector.
Definition: VMatrix.cc:240
gbl::GblPoint::addNextJacobian
void addNextJacobian(const SMatrix55 &aJac)
Define jacobian to next scatterer (by GBLTrajectory constructor)
Definition: GblPoint.cc:379
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:116
gbl::GblTrajectory::getLabels
void getLabels(std::vector< unsigned int > &aLabelList)
Get (list of) labels of points on (simple) trajectory.
Definition: GblTrajectory.cc:635
SVector5
ROOT::Math::SVector< double, 5 > SVector5
Definition: GblPoint.h:31
gbl::MilleBinary
Millepede-II (binary) record.
Definition: MilleBinary.h:46
gbl::GblTrajectory::calcJacobians
void calcJacobians()
Calculate Jacobians to previous/next scatterer from point to point ones.
Definition: GblTrajectory.cc:296
gbl::MilleBinary::writeRecord
void writeRecord()
Write record to file.
Definition: MilleBinary.cc:89
gbl::GblTrajectory::scatDataIndex
std::vector< unsigned int > scatDataIndex
mapping points to data blocks from scatterers
Definition: GblTrajectory.h:76
gbl::GblTrajectory::numMeasurements
unsigned int numMeasurements
Total number of measurements.
Definition: GblTrajectory.h:68
gbl::GblTrajectory::printTrajectory
void printTrajectory(unsigned int level=0)
Print GblTrajectory.
Definition: GblTrajectory.cc:1029
gbl::GblPoint::getOffset
int getOffset() const
Retrieve offset for point.
Definition: GblPoint.cc:349
gbl::GblTrajectory::numParameters
unsigned int numParameters
Number of fit parameters.
Definition: GblTrajectory.h:66
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:342
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:671
SVector2
ROOT::Math::SVector< double, 2 > SVector2
Definition: GblPoint.h:30
gbl::GblTrajectory::numOffsets
unsigned int numOffsets
Number of (points with) offsets on trajectory.
Definition: GblTrajectory.h:63
gbl::GblTrajectory::theVector
VVector theVector
Vector of linear equation system.
Definition: GblTrajectory.h:83
gbl::GblTrajectory::innerTransformations
std::vector< TMatrixD > innerTransformations
Transformations at innermost points of.
Definition: GblTrajectory.h:78
SMatrix22
ROOT::Math::SMatrix< double, 2 > SMatrix22
Definition: GblPoint.h:23
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:41
gbl::GblTrajectory::milleOut
void milleOut(MilleBinary &aMille)
Write valid trajectory to Millepede-II binary file.
Definition: GblTrajectory.cc:1003
gbl::GblTrajectory::theData
std::vector< GblData > theData
List of data blocks.
Definition: GblTrajectory.h:74
gbl::GblTrajectory::buildLinearEquationSystem
void buildLinearEquationSystem()
Build linear equation system from data (blocks).
Definition: GblTrajectory.cc:691
gbl::GblTrajectory::isValid
bool isValid() const
Retrieve validity of trajectory.
Definition: GblTrajectory.cc:219
gbl::GblPoint
Point on trajectory.
Definition: GblPoint.h:48
SMatrix27
ROOT::Math::SMatrix< double, 2, 7 > SMatrix27
Definition: GblData.h:22
gbl
Namespace for the general broken lines package.
Definition: BorderedBandMatrix.cc:11
gbl::BorderedBandMatrix::getBlockMatrix
TMatrixDSym getBlockMatrix(const std::vector< unsigned int > anIndex) const
Retrieve symmetric block matrix.
Definition: BorderedBandMatrix.cc:75
gbl::GblTrajectory::externalMeasurements
TVectorD externalMeasurements
Definition: GblTrajectory.h:81
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:46
gbl::GblTrajectory::numInnerTrans
unsigned int numInnerTrans
Number of inner transformations to external parameters.
Definition: GblTrajectory.h:64
gbl::GblTrajectory::getNumPoints
unsigned int getNumPoints() const
Retrieve number of point from trajectory.
Definition: GblTrajectory.cc:224
gbl::GblTrajectory::numCurvature
unsigned int numCurvature
Number of curvature parameters (0 or 1) or external parameters.
Definition: GblTrajectory.h:65
gbl::VVector::print
void print() const
Print vector.
Definition: VMatrix.cc:276
gbl::GblTrajectory::getMeasResults
unsigned int getMeasResults(unsigned int aLabel, unsigned int &numRes, TVectorD &aResiduals, TVectorD &aMeasErrors, TVectorD &aResErrors, TVectorD &aDownWeights)
Get residuals at point from measurement.
Definition: GblTrajectory.cc:586
gbl::GblTrajectory::thePoints
std::vector< std::vector< GblPoint > > thePoints
(list of) List of points on trajectory
Definition: GblTrajectory.h:73
gbl::GblTrajectory::externalPrecisions
TVectorD externalPrecisions
Definition: GblTrajectory.h:82
gbl::GblTrajectory::externalDerivatives
TMatrixD externalDerivatives
Definition: GblTrajectory.h:80
gbl::GblTrajectory::construct
void construct()
Construct trajectory from list of points.
Definition: GblTrajectory.cc:232
gbl::GblTrajectory::theDimension
std::vector< unsigned int > theDimension
List of active dimensions (0=u1, 1=u2) in fit.
Definition: GblTrajectory.h:72
gbl::GblTrajectory::getResults
unsigned int getResults(int aSignedLabel, TVectorD &localPar, TMatrixDSym &localCov) const
Get fit results at point.
Definition: GblTrajectory.cc:556
gbl::GblTrajectory::prepare
void prepare()
Prepare fit for simple or composed trajectory.
Definition: GblTrajectory.cc:712
gbl::GblPoint::getP2pJacobian
const SMatrix55 & getP2pJacobian() const
Retrieve point-to-(previous)point jacobian.
Definition: GblPoint.cc:354
SMatrix55
ROOT::Math::SMatrix< double, 5, 5 > SMatrix55
Definition: GblData.h:23
gbl::GblTrajectory::predict
void predict()
Calculate predictions for all points.
Definition: GblTrajectory.cc:925
gbl::GblTrajectory::~GblTrajectory
virtual ~GblTrajectory()
Definition: GblTrajectory.cc:215
gbl::BorderedBandMatrix::printMatrix
void printMatrix() const
Print bordered band matrix.
Definition: BorderedBandMatrix.cc:158
gbl::GblTrajectory::numPoints
std::vector< unsigned int > numPoints
Number of points on (sub)trajectory.
Definition: GblTrajectory.h:61
gbl::GblTrajectory::downWeight
double downWeight(unsigned int aMethod)
Down-weight all points.
Definition: GblTrajectory.cc:936
gbl::GblTrajectory::externalSeed
TMatrixDSym externalSeed
Precision (inverse covariance matrix) of external seed.
Definition: GblTrajectory.h:77
gbl::GblPoint::getDerivatives
void getDerivatives(int aDirection, SMatrix22 &matW, SMatrix22 &matWJ, SVector2 &vecWd) const
Retrieve derivatives of local track model.
Definition: GblPoint.cc:393
gbl::GblTrajectory::defineOffsets
void defineOffsets()
Define offsets from list of points.
Definition: GblTrajectory.cc:274
gbl::GblTrajectory::numLocals
unsigned int numLocals
Total number of (additional) local parameters.
Definition: GblTrajectory.h:67
gbl::GblTrajectory::numAllPoints
unsigned int numAllPoints
Number of all points on trajectory.
Definition: GblTrajectory.h:60
gbl::GblTrajectory::measDataIndex
std::vector< unsigned int > measDataIndex
mapping points to data blocks from measurements
Definition: GblTrajectory.h:75
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:418
gbl::BorderedBandMatrix::solveAndInvertBorderedBand
void solveAndInvertBorderedBand(const VVector &aRightHandSide, VVector &aSolution)
Solve linear equation system, partially calculate inverse.
Definition: BorderedBandMatrix.cc:125
gbl::GblTrajectory::fit
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight, std::string optionList="")
Perform fit of trajectory.
Definition: GblTrajectory.cc:955
gbl::GblTrajectory::printData
void printData()
Print GblData blocks for trajectory.
Definition: GblTrajectory.cc:1107
GblTrajectory.h
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:513
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:48
gbl::GblTrajectory::numTrajectories
unsigned int numTrajectories
Number of trajectories (in composed trajectory)
Definition: GblTrajectory.h:62
gbl::GblTrajectory::theMatrix
BorderedBandMatrix theMatrix
(Bordered band) matrix of linear equation system
Definition: GblTrajectory.h:84
gbl::GblTrajectory::getScatResults
unsigned int getScatResults(unsigned int aLabel, unsigned int &numRes, TVectorD &aResiduals, TVectorD &aMeasErrors, TVectorD &aResErrors, TVectorD &aDownWeights)
Get (kink) residuals at point from scatterer.
Definition: GblTrajectory.cc:615
gbl::BorderedBandMatrix::resize
void resize(unsigned int nSize, unsigned int nBorder=1, unsigned int nBand=5)
Resize bordered band matrix.
Definition: BorderedBandMatrix.cc:26
gbl::GblTrajectory::fitOK
bool fitOK
Trajectory has been successfully fitted (results are valid)
Definition: GblTrajectory.h:71
gbl::GblTrajectory::printPoints
void printPoints(unsigned int level=0)
Print GblPoints on trajectory.
Definition: GblTrajectory.cc:1095
gbl::GblData
Data (block) for independent scalar measurement.
Definition: GblData.h:33
gbl::GblTrajectory::constructOK
bool constructOK
Trajectory has been successfully constructed (ready for fit/output)
Definition: GblTrajectory.h:70
SMatrix25
ROOT::Math::SMatrix< double, 2, 5 > SMatrix25
Definition: GblData.h:21