Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GblTrajectory.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GblTrajectory.cc
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 }