GeneralBrokenLines  trunk_rev46
 All Classes Files Functions Variables Typedefs Pages
GblTrajectory.cpp
Go to the documentation of this file.
1 /*
2  * GblTrajectory.cpp
3  *
4  * Created on: Aug 18, 2011
5  * Author: kleinwrt
6  */
7 
91 #include "GblTrajectory.h"
92 
94 
102 GblTrajectory::GblTrajectory(const std::vector<GblPoint> &aPointList,
103  bool flagCurv, bool flagU1dir, bool flagU2dir) :
104  numAllPoints(aPointList.size()), numPoints(), numOffsets(0), numInnerTrans(
105  0), numCurvature(flagCurv ? 1 : 0), numParameters(0), numLocals(
106  0), numMeasurements(0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
107 
108  if (flagU1dir)
109  theDimension.push_back(0);
110  if (flagU2dir)
111  theDimension.push_back(1);
112  // simple (single) trajectory
113  thePoints.push_back(aPointList);
114  numPoints.push_back(numAllPoints);
115  construct(); // construct trajectory
116 }
117 
119 
130 GblTrajectory::GblTrajectory(const std::vector<GblPoint> &aPointList,
131  unsigned int aLabel, const TMatrixDSym &aSeed, bool flagCurv,
132  bool flagU1dir, bool flagU2dir) :
133  numAllPoints(aPointList.size()), numPoints(), numOffsets(0), numInnerTrans(
134  0), numCurvature(flagCurv ? 1 : 0), numParameters(0), numLocals(
135  0), numMeasurements(0), externalPoint(aLabel), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalIndex(), externalSeed(
136  aSeed), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
137 
138  if (flagU1dir)
139  theDimension.push_back(0);
140  if (flagU2dir)
141  theDimension.push_back(1);
142  // simple (single) trajectory
143  thePoints.push_back(aPointList);
144  numPoints.push_back(numAllPoints);
145  construct(); // construct trajectory
146 }
147 
149 
153 GblTrajectory::GblTrajectory(
154  const std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > &aPointsAndTransList) :
155  numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
156  aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
157  0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
158 
159  for (unsigned int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
160  thePoints.push_back(aPointsAndTransList[iTraj].first);
161  numPoints.push_back(thePoints.back().size());
162  numAllPoints += numPoints.back();
163  innerTransformations.push_back(aPointsAndTransList[iTraj].second);
164  }
165  theDimension.push_back(0);
166  theDimension.push_back(1);
167  numCurvature = innerTransformations[0].GetNcols();
168  construct(); // construct (composed) trajectory
169 }
170 
172 
179 GblTrajectory::GblTrajectory(
180  const std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > &aPointsAndTransList,
181  const TMatrixD &extDerivatives, const TVectorD &extMeasurements,
182  const TVectorD &extPrecisions) :
183  numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
184  aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
185  0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalIndex(), externalSeed(), innerTransformations(), externalDerivatives(
186  extDerivatives), externalMeasurements(extMeasurements), externalPrecisions(
187  extPrecisions) {
188 
189  for (unsigned int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
190  thePoints.push_back(aPointsAndTransList[iTraj].first);
191  numPoints.push_back(thePoints.back().size());
192  numAllPoints += numPoints.back();
193  innerTransformations.push_back(aPointsAndTransList[iTraj].second);
194  }
195  theDimension.push_back(0);
196  theDimension.push_back(1);
197  numCurvature = innerTransformations[0].GetNcols();
198  construct(); // construct (composed) trajectory
199 }
200 
201 GblTrajectory::~GblTrajectory() {
202 }
203 
205 unsigned int GblTrajectory::getNumPoints() const {
206  return numAllPoints;
207 }
208 
210 
213 void GblTrajectory::construct() {
214 
215  fitOK = false;
216  unsigned int aLabel = 0;
217  // loop over trajectories
218  numTrajectories = thePoints.size();
219  //std::cout << " numTrajectories: " << numTrajectories << ", "<< innerTransformations.size() << std::endl;
220  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
221  std::vector<GblPoint>::iterator itPoint;
222  for (itPoint = thePoints[iTraj].begin();
223  itPoint < thePoints[iTraj].end(); ++itPoint) {
224  numLocals = std::max(numLocals, itPoint->getNumLocals());
225  numMeasurements += itPoint->hasMeasurement();
226  itPoint->setLabel(++aLabel);
227  }
228  }
229  defineOffsets();
230  calcJacobians();
231  prepare();
232  // number of fit parameters
233  numParameters = (numOffsets - 2 * numInnerTrans) * theDimension.size()
234  + numCurvature + numLocals;
235 }
236 
238 
242 void GblTrajectory::defineOffsets() {
243 
244  // loop over trajectories
245  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
246  // first point is offset
247  thePoints[iTraj].front().setOffset(numOffsets++);
248  // intermediate scatterers are offsets
249  std::vector<GblPoint>::iterator itPoint;
250  for (itPoint = thePoints[iTraj].begin() + 1;
251  itPoint < thePoints[iTraj].end() - 1; ++itPoint) {
252  if (itPoint->hasScatterer()) {
253  itPoint->setOffset(numOffsets++);
254  } else {
255  itPoint->setOffset(-numOffsets);
256  }
257  }
258  // last point is offset
259  thePoints[iTraj].back().setOffset(numOffsets++);
260  }
261 }
262 
264 void GblTrajectory::calcJacobians() {
265 
266  SMatrix55 scatJacobian;
267  // loop over trajectories
268  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
269  // forward propagation (all)
270  GblPoint* previousPoint = &thePoints[iTraj].front();
271  unsigned int numStep = 0;
272  std::vector<GblPoint>::iterator itPoint;
273  for (itPoint = thePoints[iTraj].begin() + 1;
274  itPoint < thePoints[iTraj].end(); ++itPoint) {
275  if (numStep == 0) {
276  scatJacobian = itPoint->getP2pJacobian();
277  } else {
278  scatJacobian = itPoint->getP2pJacobian() * scatJacobian;
279  }
280  numStep++;
281  itPoint->addPrevJacobian(scatJacobian); // iPoint -> previous scatterer
282  if (itPoint->getOffset() >= 0) {
283  previousPoint->addNextJacobian(scatJacobian); // lastPoint -> next scatterer
284  numStep = 0;
285  previousPoint = &(*itPoint);
286  }
287  }
288  // backward propagation (without scatterers)
289  numStep = 0;
290  for (itPoint = thePoints[iTraj].end() - 1;
291  itPoint > thePoints[iTraj].begin() + 1; --itPoint) {
292  if (itPoint->getOffset() >= 0) {
293  numStep = 0;
294  continue; // skip offsets
295  }
296  if (numStep == 0) {
297  scatJacobian = itPoint->getP2pJacobian();
298  } else {
299  scatJacobian = scatJacobian * itPoint->getP2pJacobian();
300  }
301  numStep++;
302  itPoint->addNextJacobian(scatJacobian); // iPoint -> next scatterer
303  }
304  }
305 }
306 
308 
316 std::pair<std::vector<unsigned int>, TMatrixD> GblTrajectory::getJacobian(
317  int aSignedLabel) const {
318 
319  unsigned int nDim = theDimension.size();
320  unsigned int nCurv = numCurvature;
321  unsigned int nLocals = numLocals;
322  unsigned int nBorder = nCurv + nLocals;
323  unsigned int nParBRL = nBorder + 2 * nDim;
324  unsigned int nParLoc = nLocals + 5;
325  std::vector<unsigned int> anIndex;
326  anIndex.reserve(nParBRL);
327  TMatrixD aJacobian(nParLoc, nParBRL);
328  aJacobian.Zero();
329 
330  unsigned int aLabel = abs(aSignedLabel);
331  unsigned int firstLabel = 1;
332  unsigned int lastLabel = 0;
333  unsigned int aTrajectory = 0;
334  // loop over trajectories
335  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
336  aTrajectory = iTraj;
337  lastLabel += numPoints[iTraj];
338  if (aLabel <= lastLabel)
339  break;
340  if (iTraj < numTrajectories - 1)
341  firstLabel += numPoints[iTraj];
342  }
343  int nJacobian; // 0: prev, 1: next
344  // check consistency of (index, direction)
345  if (aSignedLabel > 0) {
346  nJacobian = 1;
347  if (aLabel >= lastLabel) {
348  aLabel = lastLabel;
349  nJacobian = 0;
350  }
351  } else {
352  nJacobian = 0;
353  if (aLabel <= firstLabel) {
354  aLabel = firstLabel;
355  nJacobian = 1;
356  }
357  }
358  const GblPoint aPoint = thePoints[aTrajectory][aLabel - firstLabel];
359  std::vector<unsigned int> labDer(5);
360  SMatrix55 matDer;
361  getFitToLocalJacobian(labDer, matDer, aPoint, 5, nJacobian);
362 
363  // from local parameters
364  for (unsigned int i = 0; i < nLocals; ++i) {
365  aJacobian(i + 5, i) = 1.0;
366  anIndex.push_back(i + 1);
367  }
368  // from trajectory parameters
369  unsigned int iCol = nLocals;
370  for (unsigned int i = 0; i < 5; ++i) {
371  if (labDer[i] > 0) {
372  anIndex.push_back(labDer[i]);
373  for (unsigned int j = 0; j < 5; ++j) {
374  aJacobian(j, iCol) = matDer(j, i);
375  }
376  ++iCol;
377  }
378  }
379  return std::make_pair(anIndex, aJacobian);
380 }
381 
383 
392 void GblTrajectory::getFitToLocalJacobian(std::vector<unsigned int> &anIndex,
393  SMatrix55 &aJacobian, const GblPoint &aPoint, unsigned int measDim,
394  unsigned int nJacobian) const {
395 
396  unsigned int nDim = theDimension.size();
397  unsigned int nCurv = numCurvature;
398  unsigned int nLocals = numLocals;
399 
400  int nOffset = aPoint.getOffset();
401 
402  if (nOffset < 0) // need interpolation
403  {
404  SMatrix22 prevW, prevWJ, nextW, nextWJ, matN;
405  SVector2 prevWd, nextWd;
406  int ierr;
407  aPoint.getDerivatives(0, prevW, prevWJ, prevWd); // W-, W- * J-, W- * d-
408  aPoint.getDerivatives(1, nextW, nextWJ, nextWd); // W-, W- * J-, W- * d-
409  const SMatrix22 sumWJ(prevWJ + nextWJ);
410  matN = sumWJ.Inverse(ierr); // N = (W- * J- + W+ * J+)^-1
411  // derivatives for u_int
412  const SMatrix22 prevNW(matN * prevW); // N * W-
413  const SMatrix22 nextNW(matN * nextW); // N * W+
414  const SVector2 prevNd(matN * prevWd); // N * W- * d-
415  const SVector2 nextNd(matN * nextWd); // N * W+ * d+
416 
417  unsigned int iOff = nDim * (-nOffset - 1) + nLocals + nCurv + 1; // first offset ('i' in u_i)
418 
419  // local offset
420  if (nCurv > 0) {
421  aJacobian.Place_in_col(-prevNd - nextNd, 3, 0); // from curvature
422  anIndex[0] = nLocals + 1;
423  }
424  aJacobian.Place_at(prevNW, 3, 1); // from 1st Offset
425  aJacobian.Place_at(nextNW, 3, 3); // from 2nd Offset
426  for (unsigned int i = 0; i < nDim; ++i) {
427  anIndex[1 + theDimension[i]] = iOff + i;
428  anIndex[3 + theDimension[i]] = iOff + nDim + i;
429  }
430 
431  // local slope and curvature
432  if (measDim > 2) {
433  // derivatives for u'_int
434  const SMatrix22 prevWPN(nextWJ * prevNW); // W+ * J+ * N * W-
435  const SMatrix22 nextWPN(prevWJ * nextNW); // W- * J- * N * W+
436  const SVector2 prevWNd(nextWJ * prevNd); // W+ * J+ * N * W- * d-
437  const SVector2 nextWNd(prevWJ * nextNd); // W- * J- * N * W+ * d+
438  if (nCurv > 0) {
439  aJacobian(0, 0) = 1.0;
440  aJacobian.Place_in_col(prevWNd - nextWNd, 1, 0); // from curvature
441  }
442  aJacobian.Place_at(-prevWPN, 1, 1); // from 1st Offset
443  aJacobian.Place_at(nextWPN, 1, 3); // from 2nd Offset
444  }
445  } else { // at point
446  // anIndex must be sorted
447  // forward : iOff2 = iOff1 + nDim, index1 = 1, index2 = 3
448  // backward: iOff2 = iOff1 - nDim, index1 = 3, index2 = 1
449  unsigned int iOff1 = nDim * nOffset + nCurv + nLocals + 1; // first offset ('i' in u_i)
450  unsigned int index1 = 3 - 2 * nJacobian; // index of first offset
451  unsigned int iOff2 = iOff1 + nDim * (nJacobian * 2 - 1); // second offset ('i' in u_i)
452  unsigned int index2 = 1 + 2 * nJacobian; // index of second offset
453  // local offset
454  aJacobian(3, index1) = 1.0; // from 1st Offset
455  aJacobian(4, index1 + 1) = 1.0;
456  for (unsigned int i = 0; i < nDim; ++i) {
457  anIndex[index1 + theDimension[i]] = iOff1 + i;
458  }
459 
460  // local slope and curvature
461  if (measDim > 2) {
462  SMatrix22 matW, matWJ;
463  SVector2 vecWd;
464  aPoint.getDerivatives(nJacobian, matW, matWJ, vecWd); // W, W * J, W * d
465  if (nCurv > 0) {
466  aJacobian(0, 0) = 1.0;
467  aJacobian.Place_in_col(-vecWd, 1, 0); // from curvature
468  anIndex[0] = nLocals + 1;
469  }
470  aJacobian.Place_at(-matWJ, 1, index1); // from 1st Offset
471  aJacobian.Place_at(matW, 1, index2); // from 2nd Offset
472  for (unsigned int i = 0; i < nDim; ++i) {
473  anIndex[index2 + theDimension[i]] = iOff2 + i;
474  }
475  }
476  }
477 }
478 
480 
486 void GblTrajectory::getFitToKinkJacobian(std::vector<unsigned int> &anIndex,
487  SMatrix27 &aJacobian, const GblPoint &aPoint) const {
488 
489  unsigned int nDim = theDimension.size();
490  unsigned int nCurv = numCurvature;
491  unsigned int nLocals = numLocals;
492 
493  int nOffset = aPoint.getOffset();
494 
495  SMatrix22 prevW, prevWJ, nextW, nextWJ;
496  SVector2 prevWd, nextWd;
497  aPoint.getDerivatives(0, prevW, prevWJ, prevWd); // W-, W- * J-, W- * d-
498  aPoint.getDerivatives(1, nextW, nextWJ, nextWd); // W-, W- * J-, W- * d-
499  const SMatrix22 sumWJ(prevWJ + nextWJ); // W- * J- + W+ * J+
500  const SVector2 sumWd(prevWd + nextWd); // W+ * d+ + W- * d-
501 
502  unsigned int iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1; // first offset ('i' in u_i)
503 
504  // local offset
505  if (nCurv > 0) {
506  aJacobian.Place_in_col(-sumWd, 0, 0); // from curvature
507  anIndex[0] = nLocals + 1;
508  }
509  aJacobian.Place_at(prevW, 0, 1); // from 1st Offset
510  aJacobian.Place_at(-sumWJ, 0, 3); // from 2nd Offset
511  aJacobian.Place_at(nextW, 0, 5); // from 1st Offset
512  for (unsigned int i = 0; i < nDim; ++i) {
513  anIndex[1 + theDimension[i]] = iOff + i;
514  anIndex[3 + theDimension[i]] = iOff + nDim + i;
515  anIndex[5 + theDimension[i]] = iOff + nDim * 2 + i;
516  }
517 }
518 
520 
529 unsigned int GblTrajectory::getResults(int aSignedLabel, TVectorD &localPar,
530  TMatrixDSym &localCov) const {
531  if (not fitOK)
532  return 1;
533  std::pair<std::vector<unsigned int>, TMatrixD> indexAndJacobian =
534  getJacobian(aSignedLabel);
535  unsigned int nParBrl = indexAndJacobian.first.size();
536  TVectorD aVec(nParBrl); // compressed vector
537  for (unsigned int i = 0; i < nParBrl; ++i) {
538  aVec[i] = theVector(indexAndJacobian.first[i] - 1);
539  }
540  TMatrixDSym aMat = theMatrix.getBlockMatrix(indexAndJacobian.first); // compressed matrix
541  localPar = indexAndJacobian.second * aVec;
542  localCov = aMat.Similarity(indexAndJacobian.second);
543  return 0;
544 }
545 
547 
558 unsigned int GblTrajectory::getMeasResults(unsigned int aLabel,
559  unsigned int &numData, TVectorD &aResiduals, TVectorD &aMeasErrors,
560  TVectorD &aResErrors, TVectorD &aDownWeights) {
561  numData = 0;
562  if (not fitOK)
563  return 1;
564 
565  unsigned int firstData = measDataIndex[aLabel - 1]; // first data block with measurement
566  numData = measDataIndex[aLabel] - firstData; // number of data blocks
567  for (unsigned int i = 0; i < numData; ++i) {
568  getResAndErr(firstData + i, aResiduals[i], aMeasErrors[i],
569  aResErrors[i], aDownWeights[i]);
570  }
571  return 0;
572 }
573 
575 
586 unsigned int GblTrajectory::getScatResults(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 = scatDataIndex[aLabel - 1]; // first data block with scatterer
594  numData = scatDataIndex[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 
612 void GblTrajectory::getResAndErr(unsigned int aData, double &aResidual,
613  double &aMeasError, double &aResError, double &aDownWeight) {
614 
615  double aMeasVar;
616  std::vector<unsigned int>* indLocal;
617  std::vector<double>* derLocal;
618  theData[aData].getResidual(aResidual, aMeasVar, aDownWeight, indLocal,
619  derLocal);
620  unsigned int nParBrl = (*indLocal).size();
621  TVectorD aVec(nParBrl); // compressed vector of derivatives
622  for (unsigned int j = 0; j < nParBrl; ++j) {
623  aVec[j] = (*derLocal)[j];
624  }
625  TMatrixDSym aMat = theMatrix.getBlockMatrix(*indLocal); // compressed (covariance) matrix
626  double aFitVar = aMat.Similarity(aVec); // variance from track fit
627  aMeasError = sqrt(aMeasVar); // error of measurement
628  aResError = (aFitVar < aMeasVar ? sqrt(aMeasVar - aFitVar) : 0.); // error of residual
629 }
630 
632 void GblTrajectory::buildLinearEquationSystem() {
633  unsigned int nBorder = numCurvature + numLocals;
634  theVector.resize(numParameters);
635  theMatrix.resize(numParameters, nBorder);
636  double aValue, aWeight;
637  std::vector<unsigned int>* indLocal;
638  std::vector<double>* derLocal;
639  std::vector<GblData>::iterator itData;
640  for (itData = theData.begin(); itData < theData.end(); ++itData) {
641  itData->getLocalData(aValue, aWeight, indLocal, derLocal);
642  for (unsigned int j = 0; j < indLocal->size(); ++j) {
643  theVector((*indLocal)[j] - 1) += (*derLocal)[j] * aWeight * aValue;
644  }
645  theMatrix.addBlockMatrix(aWeight, indLocal, derLocal);
646  }
647 }
648 
650 
653 void GblTrajectory::prepare() {
654  unsigned int nDim = theDimension.size();
655  // upper limit
656  unsigned int maxData = numMeasurements + nDim * (numOffsets - 2)
657  + externalSeed.GetNrows();
658  theData.reserve(maxData);
659  measDataIndex.resize(numAllPoints + 3); // include external seed and measurements
660  scatDataIndex.resize(numAllPoints + 1);
661  unsigned int nData = 0;
662  std::vector<TMatrixD> innerTransDer;
663  std::vector<std::vector<unsigned int> > innerTransLab;
664  // composed trajectory ?
665  if (numInnerTrans > 0) {
666  //std::cout << "composed trajectory" << std::endl;
667  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
668  // innermost point
669  GblPoint* innerPoint = &thePoints[iTraj].front();
670  // transformation fit to local track parameters
671  std::vector<unsigned int> firstLabels(5);
672  SMatrix55 matFitToLocal, matLocalToFit;
673  getFitToLocalJacobian(firstLabels, matFitToLocal, *innerPoint, 5);
674  // transformation local track to fit parameters
675  int ierr;
676  matLocalToFit = matFitToLocal.Inverse(ierr);
677  TMatrixD localToFit(5, 5);
678  for (unsigned int i = 0; i < 5; ++i) {
679  for (unsigned int j = 0; j < 5; ++j) {
680  localToFit(i, j) = matLocalToFit(i, j);
681  }
682  }
683  // transformation external to fit parameters at inner (first) point
684  innerTransDer.push_back(localToFit * innerTransformations[iTraj]);
685  innerTransLab.push_back(firstLabels);
686  }
687  }
688  // measurements
689  SMatrix55 matP;
690  // loop over trajectories
691  std::vector<GblPoint>::iterator itPoint;
692  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
693  for (itPoint = thePoints[iTraj].begin();
694  itPoint < thePoints[iTraj].end(); ++itPoint) {
695  SVector5 aMeas, aPrec;
696  unsigned int nLabel = itPoint->getLabel();
697  unsigned int measDim = itPoint->hasMeasurement();
698  if (measDim) {
699  const TMatrixD localDer = itPoint->getLocalDerivatives();
700  const std::vector<int> globalLab = itPoint->getGlobalLabels();
701  const TMatrixD globalDer = itPoint->getGlobalDerivatives();
702  TMatrixD transDer;
703  itPoint->getMeasurement(matP, aMeas, aPrec);
704  unsigned int iOff = 5 - measDim; // first active component
705  std::vector<unsigned int> labDer(5);
706  SMatrix55 matDer, matPDer;
707  getFitToLocalJacobian(labDer, matDer, *itPoint, measDim);
708  if (measDim > 2) {
709  matPDer = matP * matDer;
710  } else { // 'shortcut' for position measurements
711  matPDer.Place_at(
712  matP.Sub<SMatrix22>(3, 3)
713  * matDer.Sub<SMatrix25>(3, 0), 3, 0);
714  }
715 
716  if (numInnerTrans > 0) {
717  // transform for external parameters
718  TMatrixD proDer(measDim, 5);
719  // match parameters
720  unsigned int ifirst = 0;
721  unsigned int ilabel = 0;
722  while (ilabel < 5) {
723  if (labDer[ilabel] > 0) {
724  while (innerTransLab[iTraj][ifirst]
725  != labDer[ilabel] and ifirst < 5) {
726  ++ifirst;
727  }
728  if (ifirst >= 5) {
729  labDer[ilabel] -= 2 * nDim * (iTraj + 1); // adjust label
730  } else {
731  // match
732  labDer[ilabel] = 0; // mark as related to external parameters
733  for (unsigned int k = iOff; k < 5; ++k) {
734  proDer(k - iOff, ifirst) = matPDer(k,
735  ilabel);
736  }
737  }
738  }
739  ++ilabel;
740  }
741  transDer.ResizeTo(measDim, numCurvature);
742  transDer = proDer * innerTransDer[iTraj];
743  }
744  for (unsigned int i = iOff; i < 5; ++i) {
745  if (aPrec(i) > 0.) {
746  GblData aData(nLabel, aMeas(i), aPrec(i));
747  aData.addDerivatives(i, labDer, matPDer, iOff, localDer,
748  globalLab, globalDer, numLocals, transDer);
749  theData.push_back(aData);
750  nData++;
751  }
752  }
753 
754  }
755  measDataIndex[nLabel] = nData;
756  }
757  }
758 
759  // pseudo measurements from kinks
760  scatDataIndex[0] = nData;
761  scatDataIndex[1] = nData;
762  // loop over trajectories
763  for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
764  for (itPoint = thePoints[iTraj].begin() + 1;
765  itPoint < thePoints[iTraj].end() - 1; ++itPoint) {
766  SVector2 aMeas, aPrec;
767  unsigned int nLabel = itPoint->getLabel();
768  if (itPoint->hasScatterer()) {
769  itPoint->getScatterer(aMeas, aPrec);
770  TMatrixD transDer;
771  std::vector<unsigned int> labDer(7);
772  SMatrix27 matDer;
773  getFitToKinkJacobian(labDer, matDer, *itPoint);
774 
775  if (numInnerTrans > 0) {
776  // transform for external parameters
777  TMatrixD proDer(nDim, 5);
778  // match parameters
779  unsigned int ifirst = 0;
780  unsigned int ilabel = 0;
781  while (ilabel < 7) {
782  if (labDer[ilabel] > 0) {
783  while (innerTransLab[iTraj][ifirst]
784  != labDer[ilabel] and ifirst < 5) {
785  ++ifirst;
786  }
787  if (ifirst >= 5) {
788  labDer[ilabel] -= 2 * nDim * (iTraj + 1); // adjust label
789  } else {
790  // match
791  labDer[ilabel] = 0; // mark as related to external parameters
792  for (unsigned int k = 0; k < nDim; ++k) {
793  proDer(k, ifirst) = matDer(k, ilabel);
794  }
795  }
796  }
797  ++ilabel;
798  }
799  transDer.ResizeTo(nDim, numCurvature);
800  transDer = proDer * innerTransDer[iTraj];
801  }
802  for (unsigned int i = 0; i < nDim; ++i) {
803  unsigned int iDim = theDimension[i];
804  if (aPrec(iDim) > 0.) {
805  GblData aData(nLabel, aMeas(iDim), aPrec(iDim));
806  aData.addDerivatives(iDim, labDer, matDer, numLocals,
807  transDer);
808  theData.push_back(aData);
809  nData++;
810  }
811  }
812  }
813  scatDataIndex[nLabel] = nData;
814  }
815  scatDataIndex[thePoints[iTraj].back().getLabel()] = nData;
816  }
817 
818  // external seed
819  if (externalPoint > 0) {
820  std::pair<std::vector<unsigned int>, TMatrixD> indexAndJacobian =
821  getJacobian(externalPoint);
822  externalIndex = indexAndJacobian.first;
823  std::vector<double> externalDerivatives(externalIndex.size());
824  const TMatrixDSymEigen externalEigen(externalSeed);
825  const TVectorD valEigen = externalEigen.GetEigenValues();
826  TMatrixD vecEigen = externalEigen.GetEigenVectors();
827  vecEigen = vecEigen.T() * indexAndJacobian.second;
828  for (int i = 0; i < externalSeed.GetNrows(); ++i) {
829  if (valEigen(i) > 0.) {
830  for (int j = 0; j < externalSeed.GetNcols(); ++j) {
831  externalDerivatives[j] = vecEigen(i, j);
832  }
833  GblData aData(externalPoint, 0., valEigen(i));
834  aData.addDerivatives(externalIndex, externalDerivatives);
835  theData.push_back(aData);
836  nData++;
837  }
838  }
839  }
840  measDataIndex[numAllPoints + 2] = nData;
841  // external measurements
842  unsigned int nExt = externalMeasurements.GetNrows();
843  if (nExt > 0) {
844  std::vector<unsigned int> index(numCurvature);
845  std::vector<double> derivatives(numCurvature);
846  for (unsigned int iExt = 0; iExt < nExt; ++iExt) {
847  for (unsigned int iCol = 0; iCol < numCurvature; ++iCol) {
848  index[iCol] = iCol + 1;
849  derivatives[iCol] = externalDerivatives(iExt, iCol);
850  }
851  GblData aData(1U, externalMeasurements(iExt),
852  externalPrecisions(iExt));
853  aData.addDerivatives(index, derivatives);
854  theData.push_back(aData);
855  nData++;
856  }
857  }
858  measDataIndex[numAllPoints + 3] = nData;
859 }
860 
862 void GblTrajectory::predict() {
863  std::vector<GblData>::iterator itData;
864  for (itData = theData.begin(); itData < theData.end(); ++itData) {
865  itData->setPrediction(theVector);
866  }
867 }
868 
870 
873 double GblTrajectory::downWeight(unsigned int aMethod) {
874  double aLoss = 0.;
875  std::vector<GblData>::iterator itData;
876  for (itData = theData.begin(); itData < theData.end(); ++itData) {
877  aLoss += (1. - itData->setDownWeighting(aMethod));
878  }
879  return aLoss;
880 }
881 
883 
892 unsigned int GblTrajectory::fit(double &Chi2, int &Ndf, double &lostWeight,
893  std::string optionList) {
894  const double normChi2[4] = { 1.0, 0.8737, 0.9326, 0.8228 };
895  const std::string methodList = "TtHhCc";
896 
897  unsigned int aMethod = 0;
898 
899  buildLinearEquationSystem();
900  lostWeight = 0.;
901  unsigned int ierr = 0;
902  try {
903 
904  theMatrix.solveAndInvertBorderedBand(theVector, theVector);
905  predict();
906 
907  for (unsigned int i = 0; i < optionList.size(); ++i) // down weighting iterations
908  {
909  size_t aPosition = methodList.find(optionList[i]);
910  if (aPosition != std::string::npos) {
911  aMethod = aPosition / 2 + 1;
912  lostWeight = downWeight(aMethod);
913  buildLinearEquationSystem();
914  theMatrix.solveAndInvertBorderedBand(theVector, theVector);
915  predict();
916  }
917  }
918  Ndf = theData.size() - numParameters;
919  Chi2 = 0.;
920  for (unsigned int i = 0; i < theData.size(); ++i) {
921  Chi2 += theData[i].getChi2();
922  }
923  Chi2 /= normChi2[aMethod];
924  fitOK = true;
925 
926  } catch (int e) {
927  std::cout << " GblTrajectory::fit exception " << e << std::endl;
928  Chi2 = 0.;
929  Ndf = -1;
930  lostWeight = 0.;
931  ierr = e;
932  }
933  return ierr;
934 }
935 
937 void GblTrajectory::milleOut(MilleBinary &aMille) {
938  float fValue;
939  float fErr;
940  std::vector<unsigned int>* indLocal;
941  std::vector<double>* derLocal;
942  std::vector<int>* labGlobal;
943  std::vector<double>* derGlobal;
944 
945 // data: measurements, kinks and external seed
946  std::vector<GblData>::iterator itData;
947  for (itData = theData.begin(); itData != theData.end(); ++itData) {
948  itData->getAllData(fValue, fErr, indLocal, derLocal, labGlobal,
949  derGlobal);
950  aMille.addData(fValue, fErr, *indLocal, *derLocal, *labGlobal,
951  *derGlobal);
952  }
953  aMille.writeRecord();
954 }
Generated on Thu Sep 20 2012 15:45:38 for GeneralBrokenLines by   doxygen 1.8.1