GENFIT Rev: NoNumberAvailable
Loading...
Searching...
No Matches
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
132
133#include "GblTrajectory.h"
134
136namespace gbl {
137
139
147GblTrajectory::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
175GblTrajectory::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
198GblTrajectory::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
224GblTrajectory::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
254GblTrajectory::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
287
289bool GblTrajectory::isValid() const {
290 return constructOK;
291}
292
294unsigned int GblTrajectory::getNumPoints() const {
295 return numAllPoints;
296}
297
299
302void 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
344void 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
366void 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
412std::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
488void 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
583void 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
631unsigned 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
661unsigned 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
690unsigned 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
711unsigned 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
729unsigned 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
756void 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
776void 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
797void 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
1011void 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
1022double 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
1043unsigned 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
1091void 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
1117void 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
1183void 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
1195void 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}
SMatrix25
ROOT::Math::SMatrix< double, 2, 5 > SMatrix25
Definition GblData.h:43
SMatrix55
ROOT::Math::SMatrix< double, 5, 5 > SMatrix55
Definition GblData.h:45
SMatrix27
ROOT::Math::SMatrix< double, 2, 7 > SMatrix27
Definition GblData.h:44
SVector5
ROOT::Math::SVector< double, 5 > SVector5
Definition GblPoint.h:52
SVector2
ROOT::Math::SVector< double, 2 > SVector2
Definition GblPoint.h:51
SMatrix22
ROOT::Math::SMatrix< double, 2 > SMatrix22
Definition GblPoint.h:44
gbl::GblData
Data (block) for independent scalar measurement.
Definition GblData.h:55
gbl::GblData::addDerivatives
void addDerivatives(unsigned int iRow, const std::vector< unsigned int > &labDer, const SMatrix55 &matDer, unsigned int iOff, const TMatrixD &derLocal, const std::vector< int > &labGlobal, const TMatrixD &derGlobal, unsigned int nLocal, const TMatrixD &derTrans)
Add derivatives from measurement.
Definition GblData.cc:63
gbl::GblPoint
Point on trajectory.
Definition GblPoint.h:68
gbl::GblPoint::getDerivatives
void getDerivatives(int aDirection, SMatrix22 &matW, SMatrix22 &matWJ, SVector2 &vecWd) const
Retrieve derivatives of local track model.
Definition GblPoint.cc:415
gbl::GblPoint::addNextJacobian
void addNextJacobian(const SMatrix55 &aJac)
Define jacobian to next scatterer (by GBLTrajectory constructor)
Definition GblPoint.cc:401
gbl::GblPoint::getOffset
int getOffset() const
Retrieve offset for point.
Definition GblPoint.cc:371
gbl::GblTrajectory::printData
void printData()
Print GblData blocks for trajectory.
Definition GblTrajectory.cc:1195
gbl::GblTrajectory::numPoints
std::vector< unsigned int > numPoints
Number of points on (sub)trajectory.
Definition GblTrajectory.h:88
gbl::GblTrajectory::scatDataIndex
std::vector< unsigned int > scatDataIndex
mapping points to data blocks from scatterers
Definition GblTrajectory.h:103
gbl::GblTrajectory::theMatrix
BorderedBandMatrix theMatrix
(Bordered band) matrix of linear equation system
Definition GblTrajectory.h:111
gbl::GblTrajectory::thePoints
std::vector< std::vector< GblPoint > > thePoints
(list of) List of points on trajectory
Definition GblTrajectory.h:100
gbl::GblTrajectory::getFitToLocalJacobian
void getFitToLocalJacobian(std::vector< unsigned int > &anIndex, SMatrix55 &aJacobian, const GblPoint &aPoint, unsigned int measDim, unsigned int nJacobian=1) const
Get (part of) jacobian for transformation from (trajectory) fit to track parameters at point.
Definition GblTrajectory.cc:488
gbl::GblTrajectory::getLabels
unsigned int getLabels(std::vector< unsigned int > &aLabelList)
Get (list of) labels of points on (simple) valid trajectory.
Definition GblTrajectory.cc:711
gbl::GblTrajectory::numTrajectories
unsigned int numTrajectories
Number of trajectories (in composed trajectory)
Definition GblTrajectory.h:89
gbl::GblTrajectory::fit
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight, std::string optionList="")
Perform fit of (valid) trajectory.
Definition GblTrajectory.cc:1043
gbl::GblTrajectory::construct
void construct()
Construct trajectory from list of points.
Definition GblTrajectory.cc:302
gbl::GblTrajectory::externalPoint
unsigned int externalPoint
Label of external point (or 0)
Definition GblTrajectory.h:96
gbl::GblTrajectory::getMeasResults
unsigned int getMeasResults(unsigned int aLabel, unsigned int &numRes, TVectorD &aResiduals, TVectorD &aMeasErrors, TVectorD &aResErrors, TVectorD &aDownWeights)
Get residuals from fit at point for measurement.
Definition GblTrajectory.cc:661
gbl::GblTrajectory::getResults
unsigned int getResults(int aSignedLabel, TVectorD &localPar, TMatrixDSym &localCov) const
Get fit results at point.
Definition GblTrajectory.cc:631
gbl::GblTrajectory::numAllPoints
unsigned int numAllPoints
Number of all points on trajectory.
Definition GblTrajectory.h:87
gbl::GblTrajectory::downWeight
double downWeight(unsigned int aMethod)
Down-weight all points.
Definition GblTrajectory.cc:1022
gbl::GblTrajectory::theVector
VVector theVector
Vector of linear equation system.
Definition GblTrajectory.h:110
gbl::GblTrajectory::calcJacobians
void calcJacobians()
Calculate Jacobians to previous/next scatterer from point to point ones.
Definition GblTrajectory.cc:366
gbl::GblTrajectory::innerTransformations
std::vector< TMatrixD > innerTransformations
Transformations at innermost points of.
Definition GblTrajectory.h:105
gbl::GblTrajectory::getJacobian
std::pair< std::vector< unsigned int >, TMatrixD > getJacobian(int aSignedLabel) const
Get jacobian for transformation from fit to track parameters at point.
Definition GblTrajectory.cc:412
gbl::GblTrajectory::numCurvature
unsigned int numCurvature
Number of curvature parameters (0 or 1) or external parameters.
Definition GblTrajectory.h:92
gbl::GblTrajectory::predict
void predict()
Calculate predictions for all points.
Definition GblTrajectory.cc:1011
gbl::GblTrajectory::printPoints
void printPoints(unsigned int level=0)
Print GblPoints on trajectory.
Definition GblTrajectory.cc:1183
gbl::GblTrajectory::defineOffsets
void defineOffsets()
Define offsets from list of points.
Definition GblTrajectory.cc:344
gbl::GblTrajectory::numInnerTrans
unsigned int numInnerTrans
Number of inner transformations to external parameters.
Definition GblTrajectory.h:91
gbl::GblTrajectory::numLocals
unsigned int numLocals
Total number of (additional) local parameters.
Definition GblTrajectory.h:94
gbl::GblTrajectory::getFitToKinkJacobian
void getFitToKinkJacobian(std::vector< unsigned int > &anIndex, SMatrix27 &aJacobian, const GblPoint &aPoint) const
Get jacobian for transformation from (trajectory) fit to kink parameters at point.
Definition GblTrajectory.cc:583
gbl::GblTrajectory::milleOut
void milleOut(MilleBinary &aMille)
Write valid trajectory to Millepede-II binary file.
Definition GblTrajectory.cc:1091
gbl::GblTrajectory::isValid
bool isValid() const
Retrieve validity of trajectory.
Definition GblTrajectory.cc:289
gbl::GblTrajectory::GblTrajectory
GblTrajectory(const std::vector< GblPoint > &aPointList, bool flagCurv=true, bool flagU1dir=true, bool flagU2dir=true)
Create new (simple) trajectory from list of points.
Definition GblTrajectory.cc:147
gbl::GblTrajectory::prepare
void prepare()
Prepare fit for simple or composed trajectory.
Definition GblTrajectory.cc:797
gbl::GblTrajectory::theData
std::vector< GblData > theData
List of data blocks.
Definition GblTrajectory.h:101
gbl::GblTrajectory::constructOK
bool constructOK
Trajectory has been successfully constructed (ready for fit/output)
Definition GblTrajectory.h:97
gbl::GblTrajectory::externalMeasurements
TVectorD externalMeasurements
Definition GblTrajectory.h:108
gbl::GblTrajectory::measDataIndex
std::vector< unsigned int > measDataIndex
mapping points to data blocks from measurements
Definition GblTrajectory.h:102
gbl::GblTrajectory::buildLinearEquationSystem
void buildLinearEquationSystem()
Build linear equation system from data (blocks).
Definition GblTrajectory.cc:776
gbl::GblTrajectory::getScatResults
unsigned int getScatResults(unsigned int aLabel, unsigned int &numRes, TVectorD &aResiduals, TVectorD &aMeasErrors, TVectorD &aResErrors, TVectorD &aDownWeights)
Get (kink) residuals from fit at point for scatterer.
Definition GblTrajectory.cc:690
gbl::GblTrajectory::numParameters
unsigned int numParameters
Number of fit parameters.
Definition GblTrajectory.h:93
gbl::GblTrajectory::numMeasurements
unsigned int numMeasurements
Total number of measurements.
Definition GblTrajectory.h:95
gbl::GblTrajectory::externalSeed
TMatrixDSym externalSeed
Precision (inverse covariance matrix) of external seed.
Definition GblTrajectory.h:104
gbl::GblTrajectory::externalDerivatives
TMatrixD externalDerivatives
Definition GblTrajectory.h:107
gbl::GblTrajectory::printTrajectory
void printTrajectory(unsigned int level=0)
Print GblTrajectory.
Definition GblTrajectory.cc:1117
gbl::GblTrajectory::fitOK
bool fitOK
Trajectory has been successfully fitted (results are valid)
Definition GblTrajectory.h:98
gbl::GblTrajectory::numOffsets
unsigned int numOffsets
Number of (points with) offsets on trajectory.
Definition GblTrajectory.h:90
gbl::GblTrajectory::getResAndErr
void getResAndErr(unsigned int aData, double &aResidual, double &aMeadsError, double &aResError, double &aDownWeight)
Get residual and errors from data block.
Definition GblTrajectory.cc:756
gbl::GblTrajectory::theDimension
std::vector< unsigned int > theDimension
List of active dimensions (0=u1, 1=u2) in fit.
Definition GblTrajectory.h:99
gbl::GblTrajectory::externalPrecisions
TVectorD externalPrecisions
Definition GblTrajectory.h:109
gbl::GblTrajectory::getNumPoints
unsigned int getNumPoints() const
Retrieve number of point from trajectory.
Definition GblTrajectory.cc:294
gbl::MilleBinary
Millepede-II (binary) record.
Definition MilleBinary.h:68
gbl::MilleBinary::writeRecord
void writeRecord()
Write record to file.
Definition MilleBinary.cc:111
gbl::MilleBinary::addData
void addData(double aMeas, double aPrec, const std::vector< unsigned int > &indLocal, const std::vector< double > &derLocal, const std::vector< int > &labGlobal, const std::vector< double > &derGlobal)
Add data block to (end of) record.
Definition MilleBinary.cc:70
gbl
Namespace for the general broken lines package.
Definition BorderedBandMatrix.cc:33
Generated by doxygen 1.13.2