gbl::GblTrajectory Class Reference

GBL trajectory. More...

#include <GblTrajectory.h>

Public Member Functions
unsigned int	fit (double &Chi2, int &Ndf, double &lostWeight, const std::string &optionList="", unsigned int aLabel=0)
	Perform fit of (valid) trajectory. More...
	GblTrajectory (const std::vector< GblPoint > &aPointList, bool flagCurv=true, bool flagU1dir=true, bool flagU2dir=true)
	Create new (simple) trajectory from list of points. More...
template<typename Seed >
	GblTrajectory (const std::vector< GblPoint > &aPointList, unsigned int aLabel, const Eigen::MatrixBase< Seed > &aSeed, bool flagCurv=true, bool flagU1dir=true, bool flagU2dir=true)
	Create new (simple) trajectory from list of points with external seed. More...
	GblTrajectory (const std::vector< std::pair< std::vector< GblPoint >, Eigen::MatrixXd > > &aPointsAndTransList)
template<typename Derivatives , typename Measurements , typename Precision , typename std::enable_if<(Measurements::ColsAtCompileTime==1)>::type * = nullptr, typename std::enable_if<(Precision::ColsAtCompileTime!=1)>::type * = nullptr>
	GblTrajectory (const std::vector< std::pair< std::vector< GblPoint >, Eigen::MatrixXd > > &aPointsAndTransList, const Eigen::MatrixBase< Derivatives > &extDerivatives, const Eigen::MatrixBase< Measurements > &extMeasurements, const Eigen::MatrixBase< Precision > &extPrecisions)
	Create new composed trajectory from list of points and transformations with (independent or correlated) external measurements. More...
template<typename Derivatives , typename Measurements , typename Precision , typename std::enable_if<(Measurements::ColsAtCompileTime==1)>::type * = nullptr, typename std::enable_if<(Precision::ColsAtCompileTime==1)>::type * = nullptr>
	GblTrajectory (const std::vector< std::pair< std::vector< GblPoint >, Eigen::MatrixXd > > &aPointsAndTransList, const Eigen::MatrixBase< Derivatives > &extDerivatives, const Eigen::MatrixBase< Measurements > &extMeasurements, const Eigen::MatrixBase< Precision > &extPrecisions)
unsigned int	getLabels (std::vector< unsigned int > &aLabelList)
	Get (list of) labels of points on (simple) valid trajectory. More...
unsigned int	getLabels (std::vector< std::vector< unsigned int > > &aLabelList)
	Get (list of lists of) labels of points on (composed) valid trajectory. More...
unsigned int	getMeasResults (unsigned int aLabel, unsigned int &numRes, Eigen::VectorXd &aResiduals, Eigen::VectorXd &aMeasErrors, Eigen::VectorXd &aResErrors, Eigen::VectorXd &aDownWeights)
	Get residuals from fit at point for measurement. More...
unsigned int	getNumPoints () const
	Retrieve number of point from trajectory. More...
unsigned int	getResults (int aSignedLabel, Eigen::VectorXd &localPar, Eigen::MatrixXd &localCov) const
	Get fit results at point. More...
unsigned int	getScatResults (unsigned int aLabel, unsigned int &numRes, Eigen::VectorXd &aResiduals, Eigen::VectorXd &aMeasErrors, Eigen::VectorXd &aResErrors, Eigen::VectorXd &aDownWeights)
	Get (kink) residuals from fit at point for scatterer. More...
bool	isValid () const
	Retrieve validity of trajectory. More...
void	milleOut (MilleBinary &aMille)
	Write valid trajectory to Millepede-II binary file. More...
void	printData ()
	Print GblData blocks for trajectory. More...
void	printPoints (unsigned int level=0)
	Print GblPoints on trajectory. More...
void	printTrajectory (unsigned int level=0)
	Print GblTrajectory. More...
virtual	~GblTrajectory ()

Private Member Functions
void	buildLinearEquationSystem ()
	Build linear equation system from data (blocks). More...
void	calcJacobians ()
	Calculate Jacobians to previous/next scatterer from point to point ones. More...
void	construct ()
	Construct trajectory from list of points. More...
void	defineOffsets ()
	Define offsets from list of points. More...
double	downWeight (unsigned int aMethod)
	Down-weight all points. More...
void	getFitToKinkJacobian (std::array< unsigned int, 7 > &anIndex, Matrix27d &aJacobian, const GblPoint &aPoint) const
	Get jacobian for transformation from (trajectory) fit to kink parameters at point. More...
void	getFitToLocalJacobian (std::array< unsigned int, 5 > &anIndex, Matrix5d &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. More...
std::pair< std::vector< unsigned int >, Eigen::MatrixXd >	getJacobian (int aSignedLabel) const
	Get jacobian for transformation from fit to track parameters at point. More...
void	getResAndErr (unsigned int aData, bool used, double &aResidual, double &aMeadsError, double &aResError, double &aDownWeight)
	Get residual and errors from data block. More...
void	predict ()
	Calculate predictions for all points. More...
void	prepare ()
	Prepare fit for simple or composed trajectory. More...

Private Attributes
bool	constructOK
	Trajectory has been successfully constructed (ready for fit/output) More...
Eigen::MatrixXd	externalDerivatives
Eigen::VectorXd	externalMeasurements
unsigned int	externalPoint
	Label of external point (or 0) More...
Eigen::VectorXd	externalPrecisions
Eigen::MatrixXd	externalSeed
	Precision (inverse covariance matrix) of external seed. More...
bool	fitOK
	Trajectory has been successfully fitted (results are valid) More...
std::vector< Eigen::MatrixXd >	innerTransformations
	Transformations at innermost points of. More...
unsigned int	maxNumGlobals
	Max. number of global labels/derivatives per point. More...
std::vector< unsigned int >	measDataIndex
	mapping points to data blocks from measurements More...
unsigned int	numAllPoints
	Number of all points on trajectory. More...
unsigned int	numCurvature
	Number of curvature parameters (0 or 1) or external parameters. More...
unsigned int	numInnerTrans
	Number of inner transformations to external parameters. More...
unsigned int	numLocals
	Total number of (additional) local parameters. More...
unsigned int	numMeasurements
	Total number of measurements. More...
unsigned int	numOffsets
	Number of (points with) offsets on trajectory. More...
unsigned int	numParameters
	Number of fit parameters. More...
std::vector< unsigned int >	numPoints
	Number of points on (sub)trajectory. More...
unsigned int	numTrajectories
	Number of trajectories (in composed trajectory) More...
std::vector< unsigned int >	scatDataIndex
	mapping points to data blocks from scatterers More...
unsigned int	skippedMeasLabel
	Label of point with measurements skipped in fit (for unbiased residuals) (or 0) More...
std::vector< GblData >	theData
	List of data blocks. More...
std::vector< unsigned int >	theDimension
	List of active dimensions (0=u1, 1=u2) in fit. More...
BorderedBandMatrix	theMatrix
	(Bordered band) matrix of linear equation system More...
std::vector< std::vector< GblPoint > >	thePoints
	(list of) List of points on trajectory More...
VVector	theVector
	Vector of linear equation system. More...

Detailed Description

GBL trajectory.

List of GblPoints ordered by arc length. Can be fitted and optionally written to MP-II binary file.

Definition at line 47 of file GblTrajectory.h.

Constructor & Destructor Documentation

gbl::GblTrajectory::GblTrajectory	(	const std::vector< GblPoint > &	aPointList,
		bool	flagCurv = true,
		bool	flagU1dir = true,
		bool	flagU2dir = true
	)

Create new (simple) trajectory from list of points.

Curved trajectory in space (default) or without curvature (q/p) or in one plane (u-direction) only.

Parameters

[in]	aPointList	List of points
[in]	flagCurv	Use q/p
[in]	flagU1dir	Use in u1 direction
[in]	flagU2dir	Use in u2 direction

Definition at line 153 of file GblTrajectory.cc.

References construct(), externalDerivatives, externalMeasurements, externalPoint, externalPrecisions, externalSeed, plotBeamSpotDB::first, mps_fire::i, innerTransformations, maxNumGlobals, measDataIndex, eostools::move(), numAllPoints, numCurvature, numInnerTrans, numLocals, numMeasurements, numOffsets, numParameters, numPoints, scatDataIndex, edm::second(), skippedMeasLabel, theData, theDimension, and thePoints.

Referenced by GblTrajectory().

                                                                              :
    numAllPoints(aPointList.size()), numPoints(), numOffsets(0),
    numInnerTrans(0), numCurvature(flagCurv ? 1 : 0), numParameters(0),
    numLocals(0), numMeasurements(0), externalPoint(0), skippedMeasLabel(0),
    maxNumGlobals(0), theDimension(0), thePoints(), theData(), measDataIndex(),
    scatDataIndex(), externalSeed(), innerTransformations(),
    externalDerivatives(), externalMeasurements(), externalPrecisions()
  {

    if (flagU1dir)
      theDimension.push_back(0);
    if (flagU2dir)
      theDimension.push_back(1);
    // simple (single) trajectory
    thePoints.emplace_back(std::move(aPointList));
    numPoints.push_back(numAllPoints);
    construct(); // construct trajectory
  }

template<typename Seed >

gbl::GblTrajectory::GblTrajectory	(	const std::vector< GblPoint > &	aPointList,
		unsigned int	aLabel,
		const Eigen::MatrixBase< Seed > &	aSeed,
		bool	flagCurv = true,
		bool	flagU1dir = true,
		bool	flagU2dir = true
	)

Create new (simple) trajectory from list of points with external seed.

Curved trajectory in space (default) or without curvature (q/p) or in one plane (u-direction) only.

Parameters

[in]	aPointList	List of points
[in]	aLabel	(Signed) label of point for external seed (<0: in front, >0: after point, slope changes at scatterer!)
[in]	aSeed	Precision matrix of external seed
[in]	flagCurv	Use q/p
[in]	flagU1dir	Use in u1 direction
[in]	flagU2dir	Use in u2 direction

Definition at line 175 of file GblTrajectory.h.

References construct(), numAllPoints, numPoints, theDimension, and thePoints.

                                                               :
    numAllPoints(aPointList.size()), numPoints(), numOffsets(0),
    numInnerTrans(0), numCurvature(flagCurv ? 1 : 0), numParameters(0),
    numLocals(0), numMeasurements(0), externalPoint(aLabel),
    skippedMeasLabel(0), maxNumGlobals(0), theDimension(0), thePoints(),
    theData(), measDataIndex(), scatDataIndex(), externalSeed(aSeed),
    innerTransformations(), externalDerivatives(), externalMeasurements(),
    externalPrecisions()
  {

    if (flagU1dir)
      theDimension.push_back(0);
    if (flagU2dir)
      theDimension.push_back(1);
    // simple (single) trajectory
    thePoints.push_back(aPointList);
    numPoints.push_back(numAllPoints);
    construct(); // construct trajectory
  }

gbl::GblTrajectory::GblTrajectory

(

const std::vector< std::pair< std::vector< GblPoint >, Eigen::MatrixXd > > &

aPointsAndTransList

)

template<typename Derivatives , typename Measurements , typename Precision , typename std::enable_if<(Measurements::ColsAtCompileTime==1)>::type * , typename std::enable_if<(Precision::ColsAtCompileTime==1)>::type * >

gbl::GblTrajectory::GblTrajectory	(	const std::vector< std::pair< std::vector< GblPoint >, Eigen::MatrixXd > > &	aPointsAndTransList,
		const Eigen::MatrixBase< Derivatives > &	extDerivatives,
		const Eigen::MatrixBase< Measurements > &	extMeasurements,
		const Eigen::MatrixBase< Precision > &	extPrecisions
	)

Create new composed trajectory from list of points and transformations with (independent or correlated) external measurements.

Composed of curved trajectories in space. The precision matrix for the external measurements can specified as a vector for independent measurements or as arbitrary matrix which will be diagonalized.

Parameters

[in]	aPointsAndTransList	List containing pairs with list of points and transformation (at inner (first) point)
[in]	extDerivatives	Derivatives of external measurements vs external parameters
[in]	extMeasurements	External measurements (residuals)
[in]	extPrecisions	Precision of external measurements (matrix)

Composed of curved trajectories in space. The precision matrix for the external measurements can specified as a vector for independent measurements or as arbitrary matrix which will be diagonalized.

Parameters

[in]	aPointsAndTransList	List containing pairs with list of points and transformation (at inner (first) point)
[in]	extDerivatives	Derivatives of external measurements vs external parameters
[in]	extMeasurements	External measurements (residuals)
[in]	extPrecisions	Precision of external measurements (diagonal)

Definition at line 212 of file GblTrajectory.h.

References construct(), externalDerivatives, externalMeasurements, externalPoint, externalPrecisions, externalSeed, plotBeamSpotDB::first, GblTrajectory(), innerTransformations, maxNumGlobals, measDataIndex, numAllPoints, numCurvature, numInnerTrans, numLocals, numMeasurements, numOffsets, numParameters, numPoints, scatDataIndex, edm::second(), skippedMeasLabel, theData, theDimension, and thePoints.

                                                                                :
    numAllPoints(), numPoints(), numOffsets(0),
    numInnerTrans(aPointsAndTransList.size()), numParameters(0), numLocals(0),
    numMeasurements(0), externalPoint(0), skippedMeasLabel(0), maxNumGlobals(0),
    theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(),
    externalSeed(), innerTransformations()
  {
    // diagonalize external measurement
    Eigen::SelfAdjointEigenSolver<typename Precision::PlainObject> extEigen{extPrecisions};
    // @TODO   if (extEigen.info() != Success) abort();
    auto extTransformation = extEigen.eigenvectors().transpose();
    externalDerivatives.resize(extDerivatives.rows(),
                               extDerivatives.cols());
    externalDerivatives = extTransformation * extDerivatives;
    externalMeasurements.resize(extMeasurements.size());
    externalMeasurements = extTransformation * extMeasurements;
    externalPrecisions.resize(extMeasurements.size());
    externalPrecisions = extEigen.eigenvalues();

    for (unsigned int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
      thePoints.push_back(aPointsAndTransList[iTraj].first);
      numPoints.push_back(thePoints.back().size());
      numAllPoints += numPoints.back();
      innerTransformations.push_back(aPointsAndTransList[iTraj].second);
    }
    theDimension.push_back(0);
    theDimension.push_back(1);
    numCurvature = innerTransformations[0].cols();
    construct(); // construct (composed) trajectory
  }

template<typename Derivatives , typename Measurements , typename Precision , typename std::enable_if<(Measurements::ColsAtCompileTime==1)>::type * = nullptr, typename std::enable_if<(Precision::ColsAtCompileTime==1)>::type * = nullptr>

gbl::GblTrajectory::GblTrajectory	(	const std::vector< std::pair< std::vector< GblPoint >, Eigen::MatrixXd > > &	aPointsAndTransList,
		const Eigen::MatrixBase< Derivatives > &	extDerivatives,
		const Eigen::MatrixBase< Measurements > &	extMeasurements,
		const Eigen::MatrixBase< Precision > &	extPrecisions
	)

gbl::GblTrajectory::~GblTrajectory ( )

virtual

Definition at line 392 of file GblTrajectory.cc.

                                {
  }

Member Function Documentation

void gbl::GblTrajectory::buildLinearEquationSystem ( )

private

Build linear equation system from data (blocks).

Definition at line 994 of file GblTrajectory.cc.

References gbl::BorderedBandMatrix::addBlockMatrix(), gbl::InternalMeasurement, numCurvature, numLocals, numParameters, gbl::VVector::resize(), gbl::BorderedBandMatrix::resize(), skippedMeasLabel, theData, theMatrix, and theVector.

Referenced by fit().

                                                {
    unsigned int nBorder = numCurvature + numLocals;
    theVector.resize(numParameters);
    theMatrix.resize(numParameters, nBorder);
    double aValue, aWeight;
    unsigned int* indLocal;
    double* derLocal;
    unsigned int numLocal;

    std::vector<GblData>::iterator itData;
    for (itData = theData.begin(); itData < theData.end(); ++itData) {
      // skipped (internal) measurement ?
      if (itData->getLabel() == skippedMeasLabel
          && itData->getType() == InternalMeasurement)
        continue;
      itData->getLocalData(aValue, aWeight, numLocal, indLocal, derLocal);
      for (unsigned int j = 0; j < numLocal; ++j) {
        theVector(indLocal[j] - 1) += derLocal[j] * aWeight * aValue;
      }
      theMatrix.addBlockMatrix(aWeight, numLocal, indLocal, derLocal);
    }
  }

void gbl::GblTrajectory::calcJacobians ( )

private

Calculate Jacobians to previous/next scatterer from point to point ones.

Definition at line 474 of file GblTrajectory.cc.

References gbl::GblPoint::addNextJacobian(), begin, end, gbl::GblPoint::getP2pJacobian(), numTrajectories, and thePoints.

Referenced by construct().

                                    {

    Matrix5d scatJacobian;
    // loop over trajectories
    for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
      // forward propagation (all)
      GblPoint* previousPoint = &thePoints[iTraj].front();
      unsigned int numStep = 0;
      std::vector<GblPoint>::iterator itPoint;
      for (itPoint = thePoints[iTraj].begin() + 1;
           itPoint < thePoints[iTraj].end(); ++itPoint) {
        if (numStep == 0) {
          scatJacobian = itPoint->getP2pJacobian();
        } else {
          scatJacobian = itPoint->getP2pJacobian() * scatJacobian;
        }
        numStep++;
        itPoint->addPrevJacobian(scatJacobian); // iPoint -> previous scatterer
        if (itPoint->getOffset() >= 0) {
          previousPoint->addNextJacobian(scatJacobian); // lastPoint -> next scatterer
          numStep = 0;
          previousPoint = &(*itPoint);
        }
      }
      // backward propagation (without scatterers)
      for (itPoint = thePoints[iTraj].end() - 1;
           itPoint > thePoints[iTraj].begin(); --itPoint) {
        if (itPoint->getOffset() >= 0) {
          scatJacobian = itPoint->getP2pJacobian();
          continue; // skip offsets
        }
        itPoint->addNextJacobian(scatJacobian); // iPoint -> next scatterer
        scatJacobian = scatJacobian * itPoint->getP2pJacobian();
      }
    }
  }

void gbl::GblTrajectory::construct ( )

private

Construct trajectory from list of points.

Trajectory is prepared for fit or output to binary file, may consists of sub-trajectories.

Definition at line 409 of file GblTrajectory.cc.

References begin, calcJacobians(), constructOK, gather_cfg::cout, defineOffsets(), MillePedeFileConverter_cfg::e, fitOK, hpstanc_transforms::max, numAllPoints, numCurvature, numInnerTrans, numLocals, numMeasurements, numOffsets, numParameters, numTrajectories, prepare(), theDimension, and thePoints.

Referenced by GblTrajectory().

                                {

    constructOK = false;
    fitOK = false;
    unsigned int aLabel = 0;
    if (numAllPoints < 2) {
      std::cout << " GblTrajectory construction failed: too few GblPoints "
                << std::endl;
      return;
    }
    // loop over trajectories
    numTrajectories = thePoints.size();
    //std::cout << " numTrajectories: " << numTrajectories << ", " << innerTransformations.size() << std::endl;
    for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
      std::vector<GblPoint>::iterator itPoint;
      for (itPoint = thePoints[iTraj].begin();
           itPoint < thePoints[iTraj].end(); ++itPoint) {
        numLocals = std::max(numLocals, itPoint->getNumLocals());
        numMeasurements += itPoint->hasMeasurement();
        itPoint->setLabel(++aLabel);
      }
    }
    defineOffsets();
    calcJacobians();
    try {
      prepare();
    } catch (std::overflow_error &e) {
      std::cout << " GblTrajectory construction failed: " << e.what()
                << std::endl;
      return;
    }

    constructOK = true;
    // number of fit parameters
    numParameters = (numOffsets - 2 * numInnerTrans) * theDimension.size()
      + numCurvature + numLocals;
  }

void gbl::GblTrajectory::defineOffsets ( )

private

Define offsets from list of points.

Define offsets at points with scatterers and first and last point. All other points need interpolation from adjacent points with offsets.

Definition at line 452 of file GblTrajectory.cc.

References begin, numOffsets, numTrajectories, and thePoints.

Referenced by construct().

                                    {

    // loop over trajectories
    for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
      // first point is offset
      thePoints[iTraj].front().setOffset(numOffsets++);
      // intermediate scatterers are offsets
      std::vector<GblPoint>::iterator itPoint;
      for (itPoint = thePoints[iTraj].begin() + 1;
           itPoint < thePoints[iTraj].end() - 1; ++itPoint) {
        if (itPoint->hasScatterer()) {
          itPoint->setOffset(numOffsets++);
        } else {
          itPoint->setOffset(-numOffsets);
        }
      }
      // last point is offset
      thePoints[iTraj].back().setOffset(numOffsets++);
    }
  }

double gbl::GblTrajectory::downWeight ( unsigned int  aMethod )

private

Down-weight all points.

Parameters

[in]

aMethod

M-estimator (1: Tukey, 2:Huber, 3:Cauchy)

Definition at line 1253 of file GblTrajectory.cc.

References theData.

Referenced by fit().

                                                       {
    double aLoss = 0.;
    std::vector<GblData>::iterator itData;
    for (itData = theData.begin(); itData < theData.end(); ++itData) {
      aLoss += (1. - itData->setDownWeighting(aMethod));
    }
    return aLoss;
  }

unsigned int gbl::GblTrajectory::fit	(	double &	Chi2,
		int &	Ndf,
		double &	lostWeight,
		const std::string &	optionList = "",
		unsigned int	aLabel = 0
	)

Perform fit of (valid) trajectory.

Optionally iterate for outlier down-weighting. Fit may fail due to singular or not positive definite matrices (internal exceptions 1-3).

Parameters

[out]	Chi2	Chi2 sum (corrected for down-weighting)
[out]	Ndf	Number of degrees of freedom
[out]	lostWeight	Sum of weights lost due to down-weighting
[in]	optionList	Iterations for down-weighting (One character per iteration: t,h,c (or T,H,C) for Tukey, Huber or Cauchy function)
[in]	aLabel	Label of point where to skip measurements (for unbiased residuals)

Returns

Error code (non zero value indicates failure of fit)

Definition at line 1275 of file GblTrajectory.cc.

References buildLinearEquationSystem(), constructOK, gather_cfg::cout, downWeight(), MillePedeFileConverter_cfg::e, fitOK, mps_fire::i, gbl::InternalMeasurement, numParameters, predict(), skippedMeasLabel, gbl::BorderedBandMatrix::solveAndInvertBorderedBand(), AlCaHLTBitMon_QueryRunRegistry::string, theData, theMatrix, and theVector.

Referenced by trackingPlots.Iteration::modules().

                                                                                    {
    const double normChi2[4] = { 1.0, 0.8737, 0.9326, 0.8228 };
    const std::string methodList = "TtHhCc";

    Chi2 = 0.;
    Ndf = -1;
    lostWeight = 0.;
    if (not constructOK)
      return 10;

    unsigned int aMethod = 0;
    skippedMeasLabel = aLabel;

    buildLinearEquationSystem();
    lostWeight = 0.;
    unsigned int ierr = 0;
    try {

      theMatrix.solveAndInvertBorderedBand(theVector, theVector);
      predict();

      for (unsigned int i = 0; i < optionList.size(); ++i) // down weighting iterations
        {
          size_t aPosition = methodList.find(optionList[i]);
          if (aPosition != std::string::npos) {
            aMethod = aPosition / 2 + 1;
            lostWeight = downWeight(aMethod);
            buildLinearEquationSystem();
            theMatrix.solveAndInvertBorderedBand(theVector, theVector);
            predict();
          }
        }
      Ndf = -numParameters;
      Chi2 = 0.;
      for (unsigned int i = 0; i < theData.size(); ++i) {
        // skipped (internal) measurement ?
        if (theData[i].getLabel() == skippedMeasLabel
            && theData[i].getType() == InternalMeasurement)
          continue;
        Chi2 += theData[i].getChi2();
        Ndf++;
      }
      Chi2 /= normChi2[aMethod];
      fitOK = true;

    } catch (int e) {
      std::cout << " GblTrajectory::fit exception " << e << std::endl;
      ierr = e;
    }
    return ierr;
  }

void gbl::GblTrajectory::getFitToKinkJacobian ( std::array< unsigned int, 7 > &  anIndex, Matrix27d &  aJacobian, const GblPoint &  aPoint  ) const

private

Get jacobian for transformation from (trajectory) fit to kink parameters at point.

Jacobian broken lines (q/p,..,u_i-1,u_i,u_i+1..) to kink (du') parameters.

Parameters

[out]	anIndex	List of fit parameters with non zero derivatives
[out]	aJacobian	Corresponding transformation matrix
[in]	aPoint	Point to use

Definition at line 691 of file GblTrajectory.cc.

References gbl::GblPoint::getDerivatives(), gbl::GblPoint::getOffset(), mps_fire::i, numCurvature, numLocals, and theDimension.

Referenced by prepare().

                                                                                               {

    unsigned int nDim = theDimension.size();
    unsigned int nCurv = numCurvature;
    unsigned int nLocals = numLocals;

    int nOffset = aPoint.getOffset();

    aJacobian.setZero();

    Matrix2d prevW, prevWJ, nextW, nextWJ;
    Vector2d prevWd, nextWd;
    aPoint.getDerivatives(0, prevW, prevWJ, prevWd); // W-, W- * J-, W- * d-
    aPoint.getDerivatives(1, nextW, nextWJ, nextWd); // W-, W- * J-, W- * d-
    const Matrix2d sumWJ(prevWJ + nextWJ); // W- * J- + W+ * J+
    const Vector2d sumWd(prevWd + nextWd); // W+ * d+ + W- * d-

    unsigned int iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1; // first offset ('i' in u_i)

    // local offset
    if (nCurv > 0) {
      aJacobian.block<2, 1>(0, 0) = -sumWd; // from curvature
      anIndex[0] = nLocals + 1;
    }
    aJacobian.block<2, 2>(0, 1) = prevW; // from 1st Offset
    aJacobian.block<2, 2>(0, 3) = -sumWJ; // from 2nd Offset
    aJacobian.block<2, 2>(0, 5) = nextW; // from 1st Offset
    for (unsigned int i = 0; i < nDim; ++i) {
      anIndex[1 + theDimension[i]] = iOff + i;
      anIndex[3 + theDimension[i]] = iOff + nDim + i;
      anIndex[5 + theDimension[i]] = iOff + nDim * 2 + i;
    }
  }

void gbl::GblTrajectory::getFitToLocalJacobian ( std::array< unsigned int, 5 > &  anIndex, Matrix5d &  aJacobian, const GblPoint &  aPoint, unsigned int  measDim, unsigned int  nJacobian = 1  ) const

private

Get (part of) jacobian for transformation from (trajectory) fit to track parameters at point.

Jacobian broken lines (q/p,..,u_i,u_i+1..) to local (q/p,u',u) parameters.

Parameters

[out]	anIndex	List of fit parameters with non zero derivatives
[out]	aJacobian	Corresponding transformation matrix
[in]	aPoint	Point to use
[in]	measDim	Dimension of 'measurement' (<=2: calculate only offset part, >2: complete matrix)
[in]	nJacobian	Direction (0: to previous offset, 1: to next offset)

Definition at line 596 of file GblTrajectory.cc.

References gbl::GblPoint::getDerivatives(), gbl::GblPoint::getOffset(), mps_fire::i, numCurvature, numLocals, Validation_hcalonly_cfi::sign, and theDimension.

Referenced by getJacobian(), and prepare().

                                                                          {

    unsigned int nDim = theDimension.size();
    unsigned int nCurv = numCurvature;
    unsigned int nLocals = numLocals;

    int nOffset = aPoint.getOffset();

    aJacobian.setZero();
    if (nOffset < 0) // need interpolation
      {
        Matrix2d prevW, prevWJ, nextW, nextWJ, matN;
        Vector2d prevWd, nextWd;
        aPoint.getDerivatives(0, prevW, prevWJ, prevWd); // W-, W- * J-, W- * d-
        aPoint.getDerivatives(1, nextW, nextWJ, nextWd); // W+, W+ * J+, W+ * d+
        const Matrix2d sumWJ(prevWJ + nextWJ);
        matN = sumWJ.inverse(); // N = (W- * J- + W+ * J+)^-1
        // derivatives for u_int
        const Matrix2d prevNW(matN * prevW); // N * W-
        const Matrix2d nextNW(matN * nextW); // N * W+
        const Vector2d prevNd(matN * prevWd); // N * W- * d-
        const Vector2d nextNd(matN * nextWd); // N * W+ * d+

        unsigned int iOff = nDim * (-nOffset - 1) + nLocals + nCurv + 1; // first offset ('i' in u_i)

        // local offset
        if (nCurv > 0) {
          aJacobian.block<2, 1>(3, 0) = -prevNd - nextNd; // from curvature
          anIndex[0] = nLocals + 1;
        }
        aJacobian.block<2, 2>(3, 1) = prevNW; // from 1st Offset
        aJacobian.block<2, 2>(3, 3) = nextNW; // from 2nd Offset
        for (unsigned int i = 0; i < nDim; ++i) {
          anIndex[1 + theDimension[i]] = iOff + i;
          anIndex[3 + theDimension[i]] = iOff + nDim + i;
        }

        // local slope and curvature
        if (measDim > 2) {
          // derivatives for u'_int
          const Matrix2d prevWPN(nextWJ * prevNW); // W+ * J+ * N * W-
          const Matrix2d nextWPN(prevWJ * nextNW); // W- * J- * N * W+
          const Vector2d prevWNd(nextWJ * prevNd); // W+ * J+ * N * W- * d-
          const Vector2d nextWNd(prevWJ * nextNd); // W- * J- * N * W+ * d+
          if (nCurv > 0) {
            aJacobian(0, 0) = 1.0;
            aJacobian.block<2, 1>(1, 0) = prevWNd - nextWNd; // from curvature
          }
          aJacobian.block<2, 2>(1, 1) = -prevWPN; // from 1st Offset
          aJacobian.block<2, 2>(1, 3) = nextWPN; // from 2nd Offset
        }
      } else { // at point
      // anIndex must be sorted
      // forward : iOff2 = iOff1 + nDim, index1 = 1, index2 = 3
      // backward: iOff2 = iOff1 - nDim, index1 = 3, index2 = 1
      unsigned int iOff1 = nDim * nOffset + nCurv + nLocals + 1; // first offset ('i' in u_i)
      unsigned int index1 = 3 - 2 * nJacobian; // index of first offset
      unsigned int iOff2 = iOff1 + nDim * (nJacobian * 2 - 1); // second offset ('i' in u_i)
      unsigned int index2 = 1 + 2 * nJacobian; // index of second offset
      // local offset
      aJacobian(3, index1) = 1.0; // from 1st Offset
      aJacobian(4, index1 + 1) = 1.0;
      for (unsigned int i = 0; i < nDim; ++i) {
        anIndex[index1 + theDimension[i]] = iOff1 + i;
      }

      // local slope and curvature
      if (measDim > 2) {
        Matrix2d matW, matWJ;
        Vector2d vecWd;
        aPoint.getDerivatives(nJacobian, matW, matWJ, vecWd); // W, W * J, W * d
        double sign = (nJacobian > 0) ? 1. : -1.;
        if (nCurv > 0) {
          aJacobian(0, 0) = 1.0;
          aJacobian.block<2, 1>(1, 0) = -sign * vecWd; // from curvature
          anIndex[0] = nLocals + 1;
        }
        aJacobian.block<2, 2>(1, index1) = -sign * matWJ; // from 1st Offset
        aJacobian.block<2, 2>(1, index2) = sign * matW; // from 2nd Offset
        for (unsigned int i = 0; i < nDim; ++i) {
          anIndex[index2 + theDimension[i]] = iOff2 + i;
        }
      }
    }
  }

std::pair< std::vector< unsigned int >, MatrixXd > gbl::GblTrajectory::getJacobian ( int  aSignedLabel ) const

private

Get jacobian for transformation from fit to track parameters at point.

Jacobian broken lines (q/p,..,u_i,u_i+1..) to track (q/p,u',u) parameters including additional local parameters.

Parameters

[in]

aSignedLabel

(Signed) label of point for external seed (<0: in front, >0: after point, slope changes at scatterer!)

Returns

List of fit parameters with non zero derivatives and corresponding transformation matrix

Definition at line 520 of file GblTrajectory.cc.

References funct::abs(), getFitToLocalJacobian(), mps_fire::i, numCurvature, numLocals, numPoints, numTrajectories, theDimension, and thePoints.

Referenced by getResults(), getScatResults(), and prepare().

                                                                                                    {

    unsigned int nDim = theDimension.size();
    unsigned int nCurv = numCurvature;
    unsigned int nLocals = numLocals;
    unsigned int nBorder = nCurv + nLocals;
    unsigned int nParBRL = nBorder + 2 * nDim;
    unsigned int nParLoc = nLocals + 5;
    std::vector<unsigned int> anIndex;
    anIndex.reserve(nParBRL);
    MatrixXd aJacobian(nParLoc, nParBRL);
    aJacobian.setZero();

    unsigned int aLabel = abs(aSignedLabel);
    unsigned int firstLabel = 1;
    unsigned int lastLabel = 0;
    unsigned int aTrajectory = 0;
    // loop over trajectories
    for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
      aTrajectory = iTraj;
      lastLabel += numPoints[iTraj];
      if (aLabel <= lastLabel)
        break;
      if (iTraj < numTrajectories - 1)
        firstLabel += numPoints[iTraj];
    }
    int nJacobian; // 0: prev, 1: next
    // check consistency of (index, direction)
    if (aSignedLabel > 0) {
      nJacobian = 1;
      if (aLabel >= lastLabel) {
        aLabel = lastLabel;
        nJacobian = 0;
      }
    } else {
      nJacobian = 0;
      if (aLabel <= firstLabel) {
        aLabel = firstLabel;
        nJacobian = 1;
      }
    }
    const GblPoint aPoint = thePoints[aTrajectory][aLabel - firstLabel];
    std::array<unsigned int, 5> labDer;
    Matrix5d matDer;
    getFitToLocalJacobian(labDer, matDer, aPoint, 5, nJacobian);

    // from local parameters
    for (unsigned int i = 0; i < nLocals; ++i) {
      aJacobian(i + 5, i) = 1.0;
      anIndex.push_back(i + 1);
    }
    // from trajectory parameters
    unsigned int iCol = nLocals;
    for (unsigned int i = 0; i < 5; ++i) {
      if (labDer[i] > 0) {
        anIndex.push_back(labDer[i]);
        for (unsigned int j = 0; j < 5; ++j) {
          aJacobian(j, iCol) = matDer(j, i);
        }
        ++iCol;
      }
    }
    return std::make_pair(anIndex, aJacobian);
  }

unsigned int gbl::GblTrajectory::getLabels

(

std::vector< unsigned int > &

aLabelList

)

Get (list of) labels of points on (simple) valid trajectory.

Parameters

[out]

aLabelList

List of labels (aLabelList[i] = i+1)

Returns

error code (non-zero if trajectory not valid (constructed successfully))

Definition at line 924 of file GblTrajectory.cc.

References constructOK, mps_fire::i, and thePoints.

                                                                           {
    if (not constructOK)
      return 1;

    unsigned int aLabel = 0;
    unsigned int nPoint = thePoints[0].size();
    aLabelList.resize(nPoint);
    for (unsigned i = 0; i < nPoint; ++i) {
      aLabelList[i] = ++aLabel;
    }
    return 0;
  }

unsigned int gbl::GblTrajectory::getLabels

(

std::vector< std::vector< unsigned int > > &

aLabelList

)

Get (list of lists of) labels of points on (composed) valid trajectory.

Parameters

[out]

aLabelList

List of of lists of labels

Returns

error code (non-zero if trajectory not valid (constructed successfully))

Definition at line 942 of file GblTrajectory.cc.

References constructOK, mps_fire::i, numTrajectories, and thePoints.

                                                                                         {
    if (not constructOK)
      return 1;

    unsigned int aLabel = 0;
    aLabelList.resize(numTrajectories);
    for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
      unsigned int nPoint = thePoints[iTraj].size();
      aLabelList[iTraj].resize(nPoint);
      for (unsigned i = 0; i < nPoint; ++i) {
        aLabelList[iTraj][i] = ++aLabel;
      }
    }
    return 0;
  }

unsigned int gbl::GblTrajectory::getMeasResults	(	unsigned int	aLabel,
		unsigned int &	numRes,
		Eigen::VectorXd &	aResiduals,
		Eigen::VectorXd &	aMeasErrors,
		Eigen::VectorXd &	aResErrors,
		Eigen::VectorXd &	aDownWeights
	)

Get residuals from fit at point for measurement.

Get (diagonalized) residual, error of measurement and residual and down-weighting factor for measurement at point

Parameters

[in]	aLabel	Label of point on trajectory
[out]	numData	Number of data blocks from measurement at point
[out]	aResiduals	Measurements-Predictions
[out]	aMeasErrors	Errors of Measurements
[out]	aResErrors	Errors of Residuals (including correlations from track fit)
[out]	aDownWeights	Down-Weighting factors

Returns

error code (non-zero if trajectory not fitted successfully)

Definition at line 772 of file GblTrajectory.cc.

References fitOK, getResAndErr(), mps_fire::i, measDataIndex, and skippedMeasLabel.

Referenced by getScatResults().

                                                                                           {
    numData = 0;
    if (not fitOK)
      return 1;

    unsigned int firstData = measDataIndex[aLabel - 1]; // first data block with measurement
    numData = measDataIndex[aLabel] - firstData; // number of data blocks
    for (unsigned int i = 0; i < numData; ++i) {
      getResAndErr(firstData + i, (aLabel != skippedMeasLabel), aResiduals(i),
                   aMeasErrors(i), aResErrors(i), aDownWeights(i));
    }
    return 0;
  }

unsigned int gbl::GblTrajectory::getNumPoints

(

)

const

Retrieve number of point from trajectory.

Definition at line 401 of file GblTrajectory.cc.

References numAllPoints.

                                                 {
    return numAllPoints;
  }

void gbl::GblTrajectory::getResAndErr ( unsigned int  aData, bool  used, double &  aResidual, double &  aMeasError, double &  aResError, double &  aDownWeight  )

private

Get residual and errors from data block.

Get residual, error of measurement and residual and down-weighting factor for (single) data block

Parameters

[in]	aData	Label of data block
[in]	used	Flag for usage of data block in fit
[out]	aResidual	Measurement-Prediction
[out]	aMeasError	Error of Measurement
[out]	aResError	Error of Residual (including correlations from track fit)
[out]	aDownWeight	Down-Weighting factor

Definition at line 970 of file GblTrajectory.cc.

References gbl::BorderedBandMatrix::getBlockMatrix(), mathSSE::sqrt(), theData, and theMatrix.

Referenced by getMeasResults(), and getScatResults().

                                                        {

    double aMeasVar;
    unsigned int numLocal;
    unsigned int* indLocal;
    double* derLocal;
    theData[aData].getResidual(aResidual, aMeasVar, aDownWeight, numLocal,
                               indLocal, derLocal);
    VectorXd aVec(numLocal); // compressed vector of derivatives
    for (unsigned int j = 0; j < numLocal; ++j) {
      aVec[j] = derLocal[j];
    }
    MatrixXd aMat = theMatrix.getBlockMatrix(numLocal, indLocal); // compressed (covariance) matrix
    double aFitVar = aVec.transpose() * aMat * aVec; // variance from track fit
    aMeasError = sqrt(aMeasVar); // error of measurement
    if (used)
      aResError = (aFitVar < aMeasVar ? sqrt(aMeasVar - aFitVar) : 0.); // error of (biased) residual
    else
      aResError = sqrt(aMeasVar + aFitVar); // error of (unbiased) residual
  }

unsigned int gbl::GblTrajectory::getResults	(	int	aSignedLabel,
		Eigen::VectorXd &	localPar,
		Eigen::MatrixXd &	localCov
	)		const

Get fit results at point.

Get corrections and covariance matrix for local track and additional parameters in forward or backward direction.

The point is identified by its label (1..number(points)), the sign distinguishes the backward (facing previous point) and forward 'side' (facing next point). For scatterers the track direction may change in between.

Parameters

[in]	aSignedLabel	(Signed) label of point on trajectory (<0: in front, >0: after point, slope changes at scatterer!)
[out]	localPar	Corrections for local parameters
[out]	localCov	Covariance for local parameters

Returns

error code (non-zero if trajectory not fitted successfully)

Definition at line 741 of file GblTrajectory.cc.

References fitOK, gbl::BorderedBandMatrix::getBlockMatrix(), getJacobian(), mps_fire::i, theMatrix, and theVector.

Referenced by getScatResults().

                                                                   {
    if (not fitOK)
      return 1;
    std::pair<std::vector<unsigned int>, MatrixXd> indexAndJacobian =
      getJacobian(aSignedLabel);
    unsigned int nParBrl = indexAndJacobian.first.size();
    VectorXd aVec(nParBrl); // compressed vector
    for (unsigned int i = 0; i < nParBrl; ++i) {
      aVec[i] = theVector(indexAndJacobian.first[i] - 1);
    }
    MatrixXd aMat = theMatrix.getBlockMatrix(indexAndJacobian.first); // compressed matrix
    localPar = indexAndJacobian.second * aVec;
    localCov = indexAndJacobian.second * aMat
      * indexAndJacobian.second.adjoint();
    return 0;
  }

unsigned int gbl::GblTrajectory::getScatResults	(	unsigned int	aLabel,
		unsigned int &	numRes,
		Eigen::VectorXd &	aResiduals,
		Eigen::VectorXd &	aMeasErrors,
		Eigen::VectorXd &	aResErrors,
		Eigen::VectorXd &	aDownWeights
	)

Get (kink) residuals from fit at point for scatterer.

Get (diagonalized) residual, error of measurement and residual and down-weighting factor for scatterering kinks at point

Parameters

[in]	aLabel	Label of point on trajectory
[out]	numData	Number of data blocks from scatterer at point
[out]	aResiduals	(kink)Measurements-(kink)Predictions
[out]	aMeasErrors	Errors of (kink)Measurements
[out]	aResErrors	Errors of Residuals (including correlations from track fit)
[out]	aDownWeights	Down-Weighting factors

Returns

error code (non-zero if trajectory not fitted successfully)

Definition at line 801 of file GblTrajectory.cc.

References fitOK, gbl::BorderedBandMatrix::getBlockMatrix(), getJacobian(), getMeasResults(), getResAndErr(), getResults(), mps_fire::i, measDataIndex, scatDataIndex, skippedMeasLabel, theMatrix, and theVector.

                                                                                           {
    numData = 0;
    if (not fitOK)
      return 1;

    unsigned int firstData = scatDataIndex[aLabel - 1]; // first data block with scatterer
    numData = scatDataIndex[aLabel] - firstData; // number of data blocks
    for (unsigned int i = 0; i < numData; ++i) {
      getResAndErr(firstData + i, true, aResiduals(i), aMeasErrors(i),
                   aResErrors(i), aDownWeights(i));
    }
    return 0;
  }

bool gbl::GblTrajectory::isValid

(

void

)

const

Retrieve validity of trajectory.

Definition at line 396 of file GblTrajectory.cc.

References constructOK.

Referenced by ntupleDataFormat._Object::_checkIsValid(), and core.AutoHandle.AutoHandle::ReallyLoad().

                                    {
    return constructOK;
  }

void gbl::GblTrajectory::milleOut

(

MilleBinary &

aMille

)

Write valid trajectory to Millepede-II binary file.

Definition at line 1329 of file GblTrajectory.cc.

References gbl::MilleBinary::addData(), constructOK, gbl::InternalMeasurement, maxNumGlobals, theData, thePoints, and gbl::MilleBinary::writeRecord().

Referenced by MillePedeAlignmentAlgorithm::addReferenceTrajectory().

                                                  {
    double aValue;
    double aErr;
    unsigned int aTraj;
    unsigned int aPoint;
    unsigned int aRow;
    unsigned int numLocal;
    unsigned int* labLocal;
    double* derLocal;
    std::vector<int> labGlobal;
    std::vector<double> derGlobal;

    if (not constructOK)
      return;

    //   data: measurements, kinks and external seed
    labGlobal.reserve(maxNumGlobals);
    derGlobal.reserve(maxNumGlobals);
    std::vector<GblData>::iterator itData;
    for (itData = theData.begin(); itData != theData.end(); ++itData) {
      itData->getAllData(aValue, aErr, numLocal, labLocal, derLocal, aTraj,
                         aPoint, aRow);
      if (itData->getType() == InternalMeasurement)
        thePoints[aTraj][aPoint].getGlobalLabelsAndDerivatives(aRow,
                                                               labGlobal, derGlobal);
      else
        labGlobal.resize(0);
      aMille.addData(aValue, aErr, numLocal, labLocal, derLocal, labGlobal,
                     derGlobal);
    }
    aMille.writeRecord();
  }

void gbl::GblTrajectory::predict ( )

private

Calculate predictions for all points.

Definition at line 1242 of file GblTrajectory.cc.

References theData, and theVector.

Referenced by fit().

                              {
    std::vector<GblData>::iterator itData;
    for (itData = theData.begin(); itData < theData.end(); ++itData) {
      itData->setPrediction(theVector);
    }
  }

void gbl::GblTrajectory::prepare ( )

private

Prepare fit for simple or composed trajectory.

Generate data (blocks) from measurements, kinks, external seed and measurements.

Definition at line 1021 of file GblTrajectory.cc.

References gbl::GblData::addDerivatives(), begin, externalDerivatives, gbl::ExternalMeasurement, externalMeasurements, externalPoint, externalPrecisions, gbl::ExternalSeed, externalSeed, getFitToKinkJacobian(), getFitToLocalJacobian(), getJacobian(), mps_fire::i, diffTreeTool::index, innerTransformations, gbl::InternalKink, gbl::InternalMeasurement, gen::k, hpstanc_transforms::max, maxNumGlobals, measDataIndex, eostools::move(), numAllPoints, numCurvature, numInnerTrans, numLocals, numMeasurements, numOffsets, numTrajectories, scatDataIndex, theData, theDimension, thePoints, and mitigatedMETSequence_cff::U.

Referenced by construct().

                              {
    unsigned int nDim = theDimension.size();
    // upper limit
    unsigned int maxData = numMeasurements + nDim * (numOffsets - 2)
      + externalSeed.rows();
    theData.reserve(maxData);
    measDataIndex.resize(numAllPoints + 3); // include external seed and measurements
    scatDataIndex.resize(numAllPoints + 1);
    unsigned int nData = 0;
    std::vector<MatrixXd> innerTransDer;
    std::vector<std::array<unsigned int, 5> > innerTransLab;
    // composed trajectory ?
    if (numInnerTrans > 0) {
      //std::cout << "composed trajectory" << std::endl;
      for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
        // innermost point
        GblPoint* innerPoint = &thePoints[iTraj].front();
        // transformation fit to local track parameters
        std::array<unsigned int, 5> firstLabels;
        Matrix5d matFitToLocal, matLocalToFit;
        getFitToLocalJacobian(firstLabels, matFitToLocal, *innerPoint, 5);
        // transformation local track to fit parameters
        matLocalToFit = matFitToLocal.inverse();
        // transformation external to fit parameters at inner (first) point
        innerTransDer.emplace_back(matLocalToFit * innerTransformations[iTraj]);
        innerTransLab.push_back(firstLabels);
      }
    }
    Matrix5d matP;              // measurements
    std::vector<GblPoint>::iterator itPoint;
    // limit the scope of proDer:
    {
      // transform for external parameters
      Eigen::Matrix<double, Eigen::Dynamic, 5,
                    Eigen::ColMajor /* default */, 5, 5> proDer;
      // loop over trajectories
      for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
        for (itPoint = thePoints[iTraj].begin();
             itPoint < thePoints[iTraj].end(); ++itPoint) {
          Vector5d aMeas, aPrec;
          unsigned int nLabel = itPoint->getLabel();
          unsigned int measDim = itPoint->hasMeasurement();
          if (measDim) {
            const MatrixXd localDer = itPoint->getLocalDerivatives();
            maxNumGlobals = std::max(maxNumGlobals,
                                     itPoint->getNumGlobals());
            MatrixXd transDer;
            itPoint->getMeasurement(matP, aMeas, aPrec);
            double minPrecision = itPoint->getMeasPrecMin();
            unsigned int iOff = 5 - measDim; // first active component
            std::array<unsigned int, 5> labDer;
            Matrix5d matDer, matPDer;
            unsigned int nJacobian =
              (itPoint < thePoints[iTraj].end() - 1) ? 1 : 0; // last point needs backward propagation
            getFitToLocalJacobian(labDer, matDer, *itPoint, measDim,
                                  nJacobian);
            if (measDim > 2) {
              matPDer = matP * matDer;
            } else { // 'shortcut' for position measurements
              matPDer.block<2, 5>(3, 0) = matP.block<2, 2>(3, 3)
                * matDer.block<2, 5>(3, 0);
            }

            if (numInnerTrans > 0) {
              // transform for external parameters
              proDer.resize(measDim, Eigen::NoChange);
              // match parameters
              unsigned int ifirst = 0;
              unsigned int ilabel = 0;
              while (ilabel < 5) {
                if (labDer[ilabel] > 0) {
                  while (innerTransLab[iTraj][ifirst]
                         != labDer[ilabel] and ifirst < 5) {
                    ++ifirst;
                  }
                  if (ifirst >= 5) {
                    labDer[ilabel] -= 2 * nDim * (iTraj + 1); // adjust label
                  } else {
                    // match
                    labDer[ilabel] = 0; // mark as related to external parameters
                    for (unsigned int k = iOff; k < 5; ++k) {
                      proDer(k - iOff, ifirst) = matPDer(k,
                                                         ilabel);
                    }
                  }
                }
                ++ilabel;
              }
              transDer.resize(measDim, numCurvature);
              transDer = proDer * innerTransDer[iTraj];
            }
            for (unsigned int i = iOff; i < 5; ++i) {
              if (aPrec(i) > minPrecision) {
                GblData aData(nLabel, InternalMeasurement, aMeas(i),
                              aPrec(i), iTraj,
                              itPoint - thePoints[iTraj].begin());
                aData.addDerivatives(i, labDer, matPDer, iOff, localDer,
                                     numLocals, transDer);
                theData.emplace_back(aData);
                nData++;
              }
            }

          }
          measDataIndex[nLabel] = nData;
        }
      }
    } // end of scope for proDer

    Matrix2d matT;              // measurements
    // limit the scope of proDer:
    {
      // transform for external parameters
      Eigen::Matrix<double, Eigen::Dynamic, 5,
                    Eigen::ColMajor /* default */, 5, 5> proDer;
      scatDataIndex[0] = nData;
      scatDataIndex[1] = nData;
      // loop over trajectories
      for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
        for (itPoint = thePoints[iTraj].begin() + 1;
             itPoint < thePoints[iTraj].end() - 1; ++itPoint) {
          Vector2d aMeas, aPrec;
          unsigned int nLabel = itPoint->getLabel();
          if (itPoint->hasScatterer()) {
            itPoint->getScatterer(matT, aMeas, aPrec);
            MatrixXd transDer;
            std::array<unsigned int, 7> labDer;
            Matrix27d matDer, matTDer;
            getFitToKinkJacobian(labDer, matDer, *itPoint);
            matTDer = matT * matDer;
            if (numInnerTrans > 0) {
              // transform for external parameters
              proDer.resize(nDim, Eigen::NoChange);
              // match parameters
              unsigned int ifirst = 0;
              unsigned int ilabel = 0;
              while (ilabel < 7) {
                if (labDer[ilabel] > 0) {
                  while (innerTransLab[iTraj][ifirst]
                         != labDer[ilabel] and ifirst < 5) {
                    ++ifirst;
                  }
                  if (ifirst >= 5) {
                    labDer[ilabel] -= 2 * nDim * (iTraj + 1); // adjust label
                  } else {
                    // match
                    labDer[ilabel] = 0; // mark as related to external parameters
                    for (unsigned int k = 0; k < nDim; ++k) {
                      proDer(k, ifirst) = matTDer(k, ilabel);
                    }
                  }
                }
                ++ilabel;
              }
              transDer.resize(nDim, numCurvature);
              transDer = proDer * innerTransDer[iTraj];
            }
            for (unsigned int i = 0; i < nDim; ++i) {
              unsigned int iDim = theDimension[i];
              if (aPrec(iDim) > 0.) {
                GblData aData(nLabel, InternalKink, aMeas(iDim),
                              aPrec(iDim), iTraj,
                              itPoint - thePoints[iTraj].begin());
                aData.addDerivatives(iDim, labDer, matTDer, numLocals,
                                     transDer);
                theData.emplace_back(aData);
                nData++;
              }
            }
          }
          scatDataIndex[nLabel] = nData;
        }
        scatDataIndex[thePoints[iTraj].back().getLabel()] = nData;
      }
    }

    // external seed
    if (externalPoint > 0) {
      std::pair<std::vector<unsigned int>, MatrixXd> indexAndJacobian =
        getJacobian(externalPoint);
      std::vector<unsigned int> externalSeedIndex = indexAndJacobian.first;
      std::vector<double> externalSeedDerivatives(externalSeedIndex.size());
      SelfAdjointEigenSolver<MatrixXd> externalSeedEigen(externalSeed);
      VectorXd valEigen = externalSeedEigen.eigenvalues();
      MatrixXd vecEigen = externalSeedEigen.eigenvectors();
      vecEigen = vecEigen.transpose() * indexAndJacobian.second;
      for (int i = 0; i < externalSeed.rows(); ++i) {
        if (valEigen(i) > 0.) {
          for (int j = 0; j < externalSeed.cols(); ++j) {
            externalSeedDerivatives[j] = vecEigen(i, j);
          }
          GblData aData(externalPoint, ExternalSeed, 0., valEigen(i));
          aData.addDerivatives(externalSeedIndex,
                               externalSeedDerivatives);
          theData.emplace_back(std::move(aData));
          nData++;
        }
      }
    }
    measDataIndex[numAllPoints + 1] = nData;
    // external measurements
    unsigned int nExt = externalMeasurements.rows();
    if (nExt > 0) {
      std::vector<unsigned int> index(numCurvature);
      std::vector<double> derivatives(numCurvature);
      for (unsigned int iExt = 0; iExt < nExt; ++iExt) {
        for (unsigned int iCol = 0; iCol < numCurvature; ++iCol) {
          index[iCol] = iCol + 1;
          derivatives[iCol] = externalDerivatives(iExt, iCol);
        }
        GblData aData(1U, ExternalMeasurement, externalMeasurements(iExt),
                      externalPrecisions(iExt));
        aData.addDerivatives(index, derivatives);
        theData.emplace_back(std::move(aData));
        nData++;
      }
    }
    measDataIndex[numAllPoints + 2] = nData;
  }

void gbl::GblTrajectory::printData

(

)

Print GblData blocks for trajectory.

Definition at line 1446 of file GblTrajectory.cc.

References gather_cfg::cout, and theData.

                                {
    std::cout << "GblData blocks " << std::endl;
    std::vector<GblData>::iterator itData;
    for (itData = theData.begin(); itData < theData.end(); ++itData) {
      itData->printData();
    }
  }

void gbl::GblTrajectory::printPoints

(

unsigned int

level = 0

)

Print GblPoints on trajectory.

Parameters

[in]

level

print level (0: minimum, >0: more)

Definition at line 1434 of file GblTrajectory.cc.

References begin, gather_cfg::cout, numTrajectories, and thePoints.

                                                    {
    std::cout << "GblPoints " << std::endl;
    for (unsigned int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
      std::vector<GblPoint>::iterator itPoint;
      for (itPoint = thePoints[iTraj].begin();
           itPoint < thePoints[iTraj].end(); ++itPoint) {
        itPoint->printPoint(level);
      }
    }
  }

void gbl::GblTrajectory::printTrajectory

(

unsigned int

level = 0

)

Print GblTrajectory.

Parameters

[in]

level

print level (0: minimum, >0: more)

Definition at line 1366 of file GblTrajectory.cc.

References constructOK, gather_cfg::cout, externalDerivatives, externalMeasurements, externalPoint, externalPrecisions, externalSeed, fitOK, mps_fire::i, innerTransformations, numAllPoints, numInnerTrans, numMeasurements, numOffsets, numParameters, gbl::VVector::print(), gbl::BorderedBandMatrix::printMatrix(), theDimension, theMatrix, and theVector.

                                                        {
    if (numInnerTrans) {
      std::cout << "Composed GblTrajectory, " << numInnerTrans
                << " subtrajectories" << std::endl;
    } else {
      std::cout << "Simple GblTrajectory" << std::endl;
    }
    if (theDimension.size() < 2) {
      std::cout << " 2D-trajectory" << std::endl;
    }
    std::cout << " Number of GblPoints          : " << numAllPoints
              << std::endl;
    std::cout << " Number of points with offsets: " << numOffsets << std::endl;
    std::cout << " Number of fit parameters     : " << numParameters
              << std::endl;
    std::cout << " Number of measurements       : " << numMeasurements
              << std::endl;
    if (externalMeasurements.rows()) {
      std::cout << " Number of ext. measurements  : "
                << externalMeasurements.rows() << std::endl;
    }
    if (externalPoint) {
      std::cout << " Label of point with ext. seed: " << externalPoint
                << std::endl;
    }
    if (constructOK) {
      std::cout << " Constructed OK " << std::endl;
    }
    if (fitOK) {
      std::cout << " Fitted OK " << std::endl;
    }
    if (level > 0) {
      IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
      if (numInnerTrans) {
        std::cout << " Inner transformations" << std::endl;
        for (unsigned int i = 0; i < numInnerTrans; ++i) {
          std::cout << innerTransformations[i].format(CleanFmt)
                    << std::endl;
        }
      }
      if (externalMeasurements.rows()) {
        std::cout << " External measurements" << std::endl;
        std::cout << "  Measurements:" << std::endl;
        std::cout << externalMeasurements.format(CleanFmt) << std::endl;
        std::cout << "  Precisions:" << std::endl;
        std::cout << externalPrecisions.format(CleanFmt) << std::endl;
        std::cout << "  Derivatives:" << std::endl;
        std::cout << externalDerivatives.format(CleanFmt) << std::endl;
      }
      if (externalPoint) {
        std::cout << " External seed:" << std::endl;
        std::cout << externalSeed.format(CleanFmt) << std::endl;
      }
      if (fitOK) {
        std::cout << " Fit results" << std::endl;
        std::cout << "  Parameters:" << std::endl;
        theVector.print();
        std::cout << "  Covariance matrix (bordered band part):"
                  << std::endl;
        theMatrix.printMatrix();
      }
    }
  }

Member Data Documentation

bool gbl::GblTrajectory::constructOK

private

Trajectory has been successfully constructed (ready for fit/output)

Definition at line 127 of file GblTrajectory.h.

Referenced by construct(), fit(), getLabels(), isValid(), milleOut(), and printTrajectory().

Eigen::MatrixXd gbl::GblTrajectory::externalDerivatives

private

Definition at line 137 of file GblTrajectory.h.

Referenced by GblTrajectory(), prepare(), and printTrajectory().

Eigen::VectorXd gbl::GblTrajectory::externalMeasurements

private

Definition at line 138 of file GblTrajectory.h.

Referenced by GblTrajectory(), prepare(), and printTrajectory().

unsigned int gbl::GblTrajectory::externalPoint

private

Label of external point (or 0)

Definition at line 124 of file GblTrajectory.h.

Referenced by GblTrajectory(), prepare(), and printTrajectory().

Eigen::VectorXd gbl::GblTrajectory::externalPrecisions

private

Definition at line 139 of file GblTrajectory.h.

Referenced by GblTrajectory(), prepare(), and printTrajectory().

Eigen::MatrixXd gbl::GblTrajectory::externalSeed

private

Precision (inverse covariance matrix) of external seed.

Definition at line 134 of file GblTrajectory.h.

Referenced by GblTrajectory(), prepare(), and printTrajectory().

bool gbl::GblTrajectory::fitOK

private

Trajectory has been successfully fitted (results are valid)

Definition at line 128 of file GblTrajectory.h.

Referenced by construct(), fit(), getMeasResults(), getResults(), getScatResults(), and printTrajectory().

std::vector<Eigen::MatrixXd> gbl::GblTrajectory::innerTransformations

private

Transformations at innermost points of.

Definition at line 135 of file GblTrajectory.h.

Referenced by GblTrajectory(), prepare(), and printTrajectory().

unsigned int gbl::GblTrajectory::maxNumGlobals

private

Max. number of global labels/derivatives per point.

Definition at line 126 of file GblTrajectory.h.

Referenced by GblTrajectory(), milleOut(), and prepare().

std::vector<unsigned int> gbl::GblTrajectory::measDataIndex

private

mapping points to data blocks from measurements

Definition at line 132 of file GblTrajectory.h.

Referenced by GblTrajectory(), getMeasResults(), getScatResults(), and prepare().

unsigned int gbl::GblTrajectory::numAllPoints

private

Number of all points on trajectory.

Definition at line 115 of file GblTrajectory.h.

Referenced by construct(), GblTrajectory(), getNumPoints(), prepare(), and printTrajectory().

unsigned int gbl::GblTrajectory::numCurvature

private

Number of curvature parameters (0 or 1) or external parameters.

Definition at line 120 of file GblTrajectory.h.

Referenced by buildLinearEquationSystem(), construct(), GblTrajectory(), getFitToKinkJacobian(), getFitToLocalJacobian(), getJacobian(), and prepare().

unsigned int gbl::GblTrajectory::numInnerTrans

private

Number of inner transformations to external parameters.

Definition at line 119 of file GblTrajectory.h.

Referenced by construct(), GblTrajectory(), prepare(), and printTrajectory().

unsigned int gbl::GblTrajectory::numLocals

private

Total number of (additional) local parameters.

Definition at line 122 of file GblTrajectory.h.

Referenced by buildLinearEquationSystem(), construct(), GblTrajectory(), getFitToKinkJacobian(), getFitToLocalJacobian(), getJacobian(), and prepare().

unsigned int gbl::GblTrajectory::numMeasurements

private

Total number of measurements.

Definition at line 123 of file GblTrajectory.h.

Referenced by construct(), GblTrajectory(), prepare(), and printTrajectory().

unsigned int gbl::GblTrajectory::numOffsets

private

Number of (points with) offsets on trajectory.

Definition at line 118 of file GblTrajectory.h.

Referenced by construct(), defineOffsets(), GblTrajectory(), prepare(), and printTrajectory().

unsigned int gbl::GblTrajectory::numParameters

private

Number of fit parameters.

Definition at line 121 of file GblTrajectory.h.

Referenced by buildLinearEquationSystem(), construct(), fit(), GblTrajectory(), and printTrajectory().

std::vector<unsigned int> gbl::GblTrajectory::numPoints

private

Number of points on (sub)trajectory.

Definition at line 116 of file GblTrajectory.h.

Referenced by GblTrajectory(), and getJacobian().

unsigned int gbl::GblTrajectory::numTrajectories

private

Number of trajectories (in composed trajectory)

Definition at line 117 of file GblTrajectory.h.

Referenced by calcJacobians(), construct(), defineOffsets(), getJacobian(), getLabels(), prepare(), and printPoints().

std::vector<unsigned int> gbl::GblTrajectory::scatDataIndex

private

mapping points to data blocks from scatterers

Definition at line 133 of file GblTrajectory.h.

Referenced by GblTrajectory(), getScatResults(), and prepare().

unsigned int gbl::GblTrajectory::skippedMeasLabel

private

Label of point with measurements skipped in fit (for unbiased residuals) (or 0)

Definition at line 125 of file GblTrajectory.h.

Referenced by buildLinearEquationSystem(), fit(), GblTrajectory(), getMeasResults(), and getScatResults().

std::vector<GblData> gbl::GblTrajectory::theData

private

List of data blocks.

Definition at line 131 of file GblTrajectory.h.

Referenced by buildLinearEquationSystem(), downWeight(), fit(), GblTrajectory(), getResAndErr(), milleOut(), predict(), prepare(), and printData().

std::vector<unsigned int> gbl::GblTrajectory::theDimension

private

List of active dimensions (0=u1, 1=u2) in fit.

Definition at line 129 of file GblTrajectory.h.

Referenced by construct(), GblTrajectory(), getFitToKinkJacobian(), getFitToLocalJacobian(), getJacobian(), prepare(), and printTrajectory().

BorderedBandMatrix gbl::GblTrajectory::theMatrix

private

(Bordered band) matrix of linear equation system

Definition at line 141 of file GblTrajectory.h.

Referenced by buildLinearEquationSystem(), fit(), getResAndErr(), getResults(), getScatResults(), and printTrajectory().

std::vector<std::vector<GblPoint> > gbl::GblTrajectory::thePoints

private

(list of) List of points on trajectory

Definition at line 130 of file GblTrajectory.h.

Referenced by calcJacobians(), construct(), defineOffsets(), GblTrajectory(), getJacobian(), getLabels(), milleOut(), prepare(), and printPoints().

VVector gbl::GblTrajectory::theVector

private

Vector of linear equation system.

Definition at line 140 of file GblTrajectory.h.

Referenced by buildLinearEquationSystem(), fit(), getResults(), getScatResults(), predict(), and printTrajectory().