117 bool flagCurv,
bool flagU1dir,
bool flagU2dir) :
118 numAllPoints(aPointList.size()), numPoints(), numOffsets(0), numInnerTrans(
119 0), numCurvature(flagCurv ? 1 : 0), numParameters(0), numLocals(
120 0), numMeasurements(0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
145 unsigned int aLabel,
const TMatrixDSym &aSeed,
bool flagCurv,
146 bool flagU1dir,
bool flagU2dir) :
147 numAllPoints(aPointList.size()), numPoints(), numOffsets(0), numInnerTrans(
148 0), numCurvature(flagCurv ? 1 : 0), numParameters(0), numLocals(
149 0), numMeasurements(0), externalPoint(aLabel), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(
150 aSeed), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
168 const std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > &aPointsAndTransList) :
169 numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
170 aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
171 0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
173 for (
unsigned int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
174 thePoints.push_back(aPointsAndTransList[iTraj].first);
194 const std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > &aPointsAndTransList,
195 const TMatrixD &extDerivatives,
const TVectorD &extMeasurements,
196 const TVectorD &extPrecisions) :
197 numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
198 aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
199 0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(
200 extDerivatives), externalMeasurements(extMeasurements), externalPrecisions(
203 for (
unsigned int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
204 thePoints.push_back(aPointsAndTransList[iTraj].first);
236 unsigned int aLabel = 0;
238 std::cout <<
" GblTrajectory construction failed: too few GblPoints "
246 std::vector<GblPoint>::iterator itPoint;
248 itPoint <
thePoints[iTraj].end(); ++itPoint) {
251 itPoint->setLabel(++aLabel);
258 }
catch (std::overflow_error &e) {
259 std::cout <<
" GblTrajectory construction failed: " << e.what()
281 std::vector<GblPoint>::iterator itPoint;
282 for (itPoint =
thePoints[iTraj].begin() + 1;
283 itPoint <
thePoints[iTraj].end() - 1; ++itPoint) {
284 if (itPoint->hasScatterer()) {
303 unsigned int numStep = 0;
304 std::vector<GblPoint>::iterator itPoint;
305 for (itPoint =
thePoints[iTraj].begin() + 1;
306 itPoint <
thePoints[iTraj].end(); ++itPoint) {
308 scatJacobian = itPoint->getP2pJacobian();
310 scatJacobian = itPoint->getP2pJacobian() * scatJacobian;
313 itPoint->addPrevJacobian(scatJacobian);
314 if (itPoint->getOffset() >= 0) {
317 previousPoint = &(*itPoint);
321 for (itPoint =
thePoints[iTraj].end() - 1;
322 itPoint >
thePoints[iTraj].begin(); --itPoint) {
323 if (itPoint->getOffset() >= 0) {
327 itPoint->addNextJacobian(scatJacobian);
328 scatJacobian = scatJacobian * itPoint->getP2pJacobian();
343 int aSignedLabel)
const {
348 unsigned int nBorder = nCurv + nLocals;
349 unsigned int nParBRL = nBorder + 2 * nDim;
350 unsigned int nParLoc = nLocals + 5;
351 std::vector<unsigned int> anIndex;
352 anIndex.reserve(nParBRL);
353 TMatrixD aJacobian(nParLoc, nParBRL);
356 unsigned int aLabel = abs(aSignedLabel);
357 unsigned int firstLabel = 1;
358 unsigned int lastLabel = 0;
359 unsigned int aTrajectory = 0;
364 if (aLabel <= lastLabel)
371 if (aSignedLabel > 0) {
373 if (aLabel >= lastLabel) {
379 if (aLabel <= firstLabel) {
385 std::vector<unsigned int> labDer(5);
390 for (
unsigned int i = 0; i < nLocals; ++i) {
391 aJacobian(i + 5, i) = 1.0;
392 anIndex.push_back(i + 1);
395 unsigned int iCol = nLocals;
396 for (
unsigned int i = 0; i < 5; ++i) {
398 anIndex.push_back(labDer[i]);
399 for (
unsigned int j = 0; j < 5; ++j) {
400 aJacobian(j, iCol) = matDer(j, i);
405 return std::make_pair(anIndex, aJacobian);
420 unsigned int nJacobian)
const {
430 SMatrix22 prevW, prevWJ, nextW, nextWJ, matN;
436 matN = sumWJ.Inverse(ierr);
440 const SVector2 prevNd(matN * prevWd);
441 const SVector2 nextNd(matN * nextWd);
443 unsigned int iOff = nDim * (-nOffset - 1) + nLocals + nCurv + 1;
447 aJacobian.Place_in_col(-prevNd - nextNd, 3, 0);
448 anIndex[0] = nLocals + 1;
450 aJacobian.Place_at(prevNW, 3, 1);
451 aJacobian.Place_at(nextNW, 3, 3);
452 for (
unsigned int i = 0; i < nDim; ++i) {
460 const SMatrix22 prevWPN(nextWJ * prevNW);
461 const SMatrix22 nextWPN(prevWJ * nextNW);
462 const SVector2 prevWNd(nextWJ * prevNd);
463 const SVector2 nextWNd(prevWJ * nextNd);
465 aJacobian(0, 0) = 1.0;
466 aJacobian.Place_in_col(prevWNd - nextWNd, 1, 0);
468 aJacobian.Place_at(-prevWPN, 1, 1);
469 aJacobian.Place_at(nextWPN, 1, 3);
475 unsigned int iOff1 = nDim * nOffset + nCurv + nLocals + 1;
476 unsigned int index1 = 3 - 2 * nJacobian;
477 unsigned int iOff2 = iOff1 + nDim * (nJacobian * 2 - 1);
478 unsigned int index2 = 1 + 2 * nJacobian;
480 aJacobian(3, index1) = 1.0;
481 aJacobian(4, index1 + 1) = 1.0;
482 for (
unsigned int i = 0; i < nDim; ++i) {
491 double sign = (nJacobian > 0) ? 1. : -1.;
493 aJacobian(0, 0) = 1.0;
494 aJacobian.Place_in_col(-sign * vecWd, 1, 0);
495 anIndex[0] = nLocals + 1;
497 aJacobian.Place_at(-sign * matWJ, 1, index1);
498 aJacobian.Place_at(sign * matW, 1, index2);
499 for (
unsigned int i = 0; i < nDim; ++i) {
527 const SVector2 sumWd(prevWd + nextWd);
529 unsigned int iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1;
533 aJacobian.Place_in_col(-sumWd, 0, 0);
534 anIndex[0] = nLocals + 1;
536 aJacobian.Place_at(prevW, 0, 1);
537 aJacobian.Place_at(-sumWJ, 0, 3);
538 aJacobian.Place_at(nextW, 0, 5);
539 for (
unsigned int i = 0; i < nDim; ++i) {
557 TMatrixDSym &localCov)
const {
560 std::pair<std::vector<unsigned int>, TMatrixD> indexAndJacobian =
562 unsigned int nParBrl = indexAndJacobian.first.size();
563 TVectorD aVec(nParBrl);
564 for (
unsigned int i = 0; i < nParBrl; ++i) {
565 aVec[i] =
theVector(indexAndJacobian.first[i] - 1);
568 localPar = indexAndJacobian.second * aVec;
569 localCov = aMat.Similarity(indexAndJacobian.second);
587 unsigned int &numData, TVectorD &aResiduals, TVectorD &aMeasErrors,
588 TVectorD &aResErrors, TVectorD &aDownWeights) {
595 for (
unsigned int i = 0; i < numData; ++i) {
596 getResAndErr(firstData + i, aResiduals[i], aMeasErrors[i],
597 aResErrors[i], aDownWeights[i]);
616 unsigned int &numData, TVectorD &aResiduals, TVectorD &aMeasErrors,
617 TVectorD &aResErrors, TVectorD &aDownWeights) {
624 for (
unsigned int i = 0; i < numData; ++i) {
625 getResAndErr(firstData + i, aResiduals[i], aMeasErrors[i],
626 aResErrors[i], aDownWeights[i]);
636 unsigned int aLabel = 0;
637 unsigned int nPoint =
thePoints[0].size();
638 aLabelList.resize(nPoint);
639 for (
unsigned i = 0; i < nPoint; ++i) {
640 aLabelList[i] = ++aLabel;
649 std::vector<std::vector<unsigned int> > &aLabelList) {
650 unsigned int aLabel = 0;
653 unsigned int nPoint =
thePoints[iTraj].size();
654 aLabelList[iTraj].resize(nPoint);
655 for (
unsigned i = 0; i < nPoint; ++i) {
656 aLabelList[iTraj][i] = ++aLabel;
672 double &aMeasError,
double &aResError,
double &aDownWeight) {
675 std::vector<unsigned int>* indLocal;
676 std::vector<double>* derLocal;
677 theData[aData].getResidual(aResidual, aMeasVar, aDownWeight, indLocal,
679 unsigned int nParBrl = (*indLocal).size();
680 TVectorD aVec(nParBrl);
681 for (
unsigned int j = 0; j < nParBrl; ++j) {
682 aVec[j] = (*derLocal)[j];
685 double aFitVar = aMat.Similarity(aVec);
686 aMeasError = sqrt(aMeasVar);
687 aResError = (aFitVar < aMeasVar ? sqrt(aMeasVar - aFitVar) : 0.);
695 double aValue, aWeight;
696 std::vector<unsigned int>* indLocal;
697 std::vector<double>* derLocal;
698 std::vector<GblData>::iterator itData;
699 for (itData =
theData.begin(); itData <
theData.end(); ++itData) {
700 itData->getLocalData(aValue, aWeight, indLocal, derLocal);
701 for (
unsigned int j = 0; j < indLocal->size(); ++j) {
702 theVector((*indLocal)[j] - 1) += (*derLocal)[j] * aWeight * aValue;
720 unsigned int nData = 0;
721 std::vector<TMatrixD> innerTransDer;
722 std::vector<std::vector<unsigned int> > innerTransLab;
730 std::vector<unsigned int> firstLabels(5);
735 matLocalToFit = matFitToLocal.Inverse(ierr);
736 TMatrixD localToFit(5, 5);
737 for (
unsigned int i = 0; i < 5; ++i) {
738 for (
unsigned int j = 0; j < 5; ++j) {
739 localToFit(i, j) = matLocalToFit(i, j);
744 innerTransLab.push_back(firstLabels);
750 std::vector<GblPoint>::iterator itPoint;
753 itPoint <
thePoints[iTraj].end(); ++itPoint) {
755 unsigned int nLabel = itPoint->getLabel();
756 unsigned int measDim = itPoint->hasMeasurement();
758 const TMatrixD localDer = itPoint->getLocalDerivatives();
759 const std::vector<int> globalLab = itPoint->getGlobalLabels();
760 const TMatrixD globalDer = itPoint->getGlobalDerivatives();
762 itPoint->getMeasurement(matP, aMeas, aPrec);
763 unsigned int iOff = 5 - measDim;
764 std::vector<unsigned int> labDer(5);
766 unsigned int nJacobian =
767 (itPoint <
thePoints[iTraj].end() - 1) ? 1 : 0;
771 matPDer = matP * matDer;
780 TMatrixD proDer(measDim, 5);
782 unsigned int ifirst = 0;
783 unsigned int ilabel = 0;
785 if (labDer[ilabel] > 0) {
786 while (innerTransLab[iTraj][ifirst]
787 != labDer[ilabel] and ifirst < 5) {
791 labDer[ilabel] -= 2 * nDim * (iTraj + 1);
795 for (
unsigned int k = iOff; k < 5; ++k) {
796 proDer(k - iOff, ifirst) = matPDer(k,
804 transDer = proDer * innerTransDer[iTraj];
806 for (
unsigned int i = iOff; i < 5; ++i) {
808 GblData aData(nLabel, aMeas(i), aPrec(i));
810 globalLab, globalDer,
numLocals, transDer);
827 for (itPoint =
thePoints[iTraj].begin() + 1;
828 itPoint <
thePoints[iTraj].end() - 1; ++itPoint) {
830 unsigned int nLabel = itPoint->getLabel();
831 if (itPoint->hasScatterer()) {
832 itPoint->getScatterer(matT, aMeas, aPrec);
834 std::vector<unsigned int> labDer(7);
837 matTDer = matT * matDer;
840 TMatrixD proDer(nDim, 5);
842 unsigned int ifirst = 0;
843 unsigned int ilabel = 0;
845 if (labDer[ilabel] > 0) {
846 while (innerTransLab[iTraj][ifirst]
847 != labDer[ilabel] and ifirst < 5) {
851 labDer[ilabel] -= 2 * nDim * (iTraj + 1);
855 for (
unsigned int k = 0; k < nDim; ++k) {
856 proDer(k, ifirst) = matTDer(k, ilabel);
863 transDer = proDer * innerTransDer[iTraj];
865 for (
unsigned int i = 0; i < nDim; ++i) {
867 if (aPrec(iDim) > 0.) {
868 GblData aData(nLabel, aMeas(iDim), aPrec(iDim));
883 std::pair<std::vector<unsigned int>, TMatrixD> indexAndJacobian =
885 std::vector<unsigned int> externalSeedIndex = indexAndJacobian.first;
886 std::vector<double> externalSeedDerivatives(externalSeedIndex.size());
888 const TVectorD valEigen = externalSeedEigen.GetEigenValues();
889 TMatrixD vecEigen = externalSeedEigen.GetEigenVectors();
890 vecEigen = vecEigen.T() * indexAndJacobian.second;
892 if (valEigen(i) > 0.) {
894 externalSeedDerivatives[j] = vecEigen(i, j);
909 for (
unsigned int iExt = 0; iExt < nExt; ++iExt) {
910 for (
unsigned int iCol = 0; iCol <
numCurvature; ++iCol) {
911 index[iCol] = iCol + 1;
926 std::vector<GblData>::iterator itData;
927 for (itData =
theData.begin(); itData <
theData.end(); ++itData) {
938 std::vector<GblData>::iterator itData;
939 for (itData =
theData.begin(); itData <
theData.end(); ++itData) {
940 aLoss += (1. - itData->setDownWeighting(aMethod));
956 std::string optionList) {
957 const double normChi2[4] = { 1.0, 0.8737, 0.9326, 0.8228 };
958 const std::string methodList =
"TtHhCc";
966 unsigned int aMethod = 0;
970 unsigned int ierr = 0;
976 for (
unsigned int i = 0; i < optionList.size(); ++i)
978 size_t aPosition = methodList.find(optionList[i]);
979 if (aPosition != std::string::npos) {
980 aMethod = aPosition / 2 + 1;
989 for (
unsigned int i = 0; i <
theData.size(); ++i) {
992 Chi2 /= normChi2[aMethod];
996 std::cout <<
" GblTrajectory::fit exception " << e << std::endl;
1006 std::vector<unsigned int>* indLocal;
1007 std::vector<double>* derLocal;
1008 std::vector<int>* labGlobal;
1009 std::vector<double>* derGlobal;
1015 std::vector<GblData>::iterator itData;
1016 for (itData =
theData.begin(); itData !=
theData.end(); ++itData) {
1017 itData->getAllData(aValue, aErr, indLocal, derLocal, labGlobal,
1019 aMille.
addData(aValue, aErr, *indLocal, *derLocal, *labGlobal,
1032 <<
" subtrajectories" << std::endl;
1034 std::cout <<
"Simple GblTrajectory" << std::endl;
1037 std::cout <<
" 2D-trajectory" << std::endl;
1041 std::cout <<
" Number of points with offsets: " <<
numOffsets << std::endl;
1042 std::cout <<
" Number of fit parameters : " <<
numParameters
1047 std::cout <<
" Number of ext. measurements : "
1051 std::cout <<
" Label of point with ext. seed: " <<
externalPoint
1055 std::cout <<
" Constructed OK " << std::endl;
1058 std::cout <<
" Fitted OK " << std::endl;
1062 std::cout <<
" Inner transformations" << std::endl;
1068 std::cout <<
" External measurements" << std::endl;
1069 std::cout <<
" Measurements:" << std::endl;
1071 std::cout <<
" Precisions:" << std::endl;
1073 std::cout <<
" Derivatives:" << std::endl;
1077 std::cout <<
" External seed:" << std::endl;
1081 std::cout <<
" Fit results" << std::endl;
1082 std::cout <<
" Parameters:" << std::endl;
1084 std::cout <<
" Covariance matrix (bordered band part):"
1096 std::cout <<
"GblPoints " << std::endl;
1098 std::vector<GblPoint>::iterator itPoint;
1099 for (itPoint =
thePoints[iTraj].begin();
1100 itPoint <
thePoints[iTraj].end(); ++itPoint) {
1101 itPoint->printPoint(level);
1108 std::cout <<
"GblData blocks " << std::endl;
1109 std::vector<GblData>::iterator itData;
1110 for (itData =
theData.begin(); itData <
theData.end(); ++itData) {
1111 itData->printData();