Multiple measurements per point implemented

git-svn-id: http://svnsrv.desy.de/public/GeneralBrokenLines/trunk@153 281f6f2b-e318-4fd1-8bce-1a4ba7aab212

cpp/CMakeLists.txt

@@ -6,8 +6,8 @@ PROJECT(GBL)
 
 # project version
 SET( ${PROJECT_NAME}_VERSION_MAJOR 2 )
-SET( ${PROJECT_NAME}_VERSION_MINOR 2 )
-SET( ${PROJECT_NAME}_VERSION_PATCH 1 )
+SET( ${PROJECT_NAME}_VERSION_MINOR 3 )
+SET( ${PROJECT_NAME}_VERSION_PATCH 0 )
 
 # make life easier and simply use the ilcsoft default settings
 # load default ilcsoft settings (install prefix, build type, rpath, etc.)

cpp/include/GblData.h

@@ -13,7 +13,7 @@
  *
  *
  *  \copyright
- *  Copyright (c) 2011 - 2016 Deutsches Elektronen-Synchroton,
+ *  Copyright (c) 2011 - 2021 Deutsches Elektronen-Synchroton,
  *  Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
  *  This library is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU Library General Public License as
@@ -57,8 +57,9 @@ enum dataBlockType {
  */
 class GblData {
 public:
-	GblData(unsigned int aLabel, dataBlockType aType, double aMeas,
-			double aPrec, unsigned int aTraj = 0, unsigned int aPoint = 0);
+	GblData(unsigned int aLabel, dataBlockType aType, double aValue,
+			double aPrec, unsigned int aTraj = 0, unsigned int aPoint = 0,
+			unsigned int aMeas = 0);
 	GblData(const GblData&) = default;
 	GblData& operator=(const GblData&) = default;
 	GblData(GblData&&) = default;
@@ -80,11 +81,11 @@ public:
 	 */
 	template<typename LocalDerivative, typename ExtDerivative>
 	void addDerivatives(unsigned int iRow,
-			const std::array<unsigned int, 5>& labDer, const Matrix5d &matDer,
+			const std::array<unsigned int, 5> &labDer, const Matrix5d &matDer,
 			unsigned int iOff,
-			const Eigen::MatrixBase<LocalDerivative>& derLocal,
+			const Eigen::MatrixBase<LocalDerivative> &derLocal,
 			unsigned int extOff,
-			const Eigen::MatrixBase<ExtDerivative>& extDer);
+			const Eigen::MatrixBase<ExtDerivative> &extDer);
 
 	/// Add derivatives from kink.
 	/**
@@ -98,9 +99,9 @@ public:
 	 */
 	template<typename ExtDerivative>
 	void addDerivatives(unsigned int iRow,
-			const std::array<unsigned int, 7>& labDer, const Matrix27d &matDer,
+			const std::array<unsigned int, 7> &labDer, const Matrix27d &matDer,
 			unsigned int extOff,
-			const Eigen::MatrixBase<ExtDerivative>& extDer);
+			const Eigen::MatrixBase<ExtDerivative> &extDer);
 
 	void addDerivatives(const std::vector<unsigned int> &index,
 			const std::vector<double> &derivatives);
@@ -111,14 +112,14 @@ public:
 	void printData() const;
 	unsigned int getLabel() const;
 	dataBlockType getType() const;
-	unsigned int getNumSimple() const;
 	void getLocalData(double &aValue, double &aWeight, unsigned int &numLocal,
-			unsigned int* &indLocal, double* &derLocal);
+			unsigned int *&indLocal, double *&derLocal);
 	void getAllData(double &aValue, double &aErr, unsigned int &numLocal,
-			unsigned int* &indLocal, double* &derLocal, unsigned int &aTraj,
-			unsigned int &aPoint, unsigned int &aRow);
+			unsigned int *&indLocal, double *&derLocal, unsigned int &aTraj,
+			unsigned int &aPoint, unsigned int &aMeas, unsigned int &aRow);
 	void getResidual(double &aResidual, double &aVariance, double &aDownWeight,
-			unsigned int &numLocal, unsigned int* &indLocal, double* &derLocal);
+			unsigned int &numLocal, unsigned int *&indLocal, double *&derLocal);
+	void getResidual(double &aResidual, double &aVariance);
 
 private:
 	unsigned int theLabel; ///< Label (of corresponding point)
@@ -128,6 +129,7 @@ private:
 	double thePrecision; ///< Precision (1/sigma**2)
 	unsigned int theTrajectory; ///< Trajectory number
 	unsigned int thePoint; ///< Point number (on trajectory)
+	unsigned int theMeas; ///< Measurement number (at point)
 	unsigned int theDWMethod; ///< Down-weighting method (0: None, 1: Tukey, 2: Huber, 3: Cauchy)
 	double theDownWeight; ///< Down-weighting factor (0-1)
 	double thePrediction; ///< Prediction from fit
@@ -142,9 +144,9 @@ private:
 
 template<typename LocalDerivative, typename ExtDerivative>
 void GblData::addDerivatives(unsigned int iRow,
-		const std::array<unsigned int, 5>& labDer, const Matrix5d &matDer,
-		unsigned int iOff, const Eigen::MatrixBase<LocalDerivative>& derLocal,
-		unsigned int extOff, const Eigen::MatrixBase<ExtDerivative>& extDer) {
+		const std::array<unsigned int, 5> &labDer, const Matrix5d &matDer,
+		unsigned int iOff, const Eigen::MatrixBase<LocalDerivative> &derLocal,
+		unsigned int extOff, const Eigen::MatrixBase<ExtDerivative> &extDer) {
 
 	unsigned int nParMax = 5 + derLocal.cols() + extDer.cols();
 	theRow = iRow - iOff;
@@ -208,8 +210,8 @@ void GblData::addDerivatives(unsigned int iRow,
 
 template<typename ExtDerivative>
 void GblData::addDerivatives(unsigned int iRow,
-		const std::array<unsigned int, 7>& labDer, const Matrix27d &matDer,
-		unsigned int extOff, const Eigen::MatrixBase<ExtDerivative>& extDer) {
+		const std::array<unsigned int, 7> &labDer, const Matrix27d &matDer,
+		unsigned int extOff, const Eigen::MatrixBase<ExtDerivative> &extDer) {
 
 	unsigned int nParMax = 7 + extDer.cols();
 	theRow = iRow;

cpp/include/GblMeasurement.h

+/*
+ * GblMeasurement.h
+ *
+ *  Created on: 31 Mar 2021
+ *      Author: kleinwrt
+ */
+
+/** \file
+ *  GblMeasurement (with optional derivatives) definition.
+ *
+ *  \author Claus Kleinwort, DESY, 2021 (Claus.Kleinwort@desy.de)
+ *
+ *
+ *  \copyright
+ *  Copyright (c) 2021 Deutsches Elektronen-Synchroton,
+ *  Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
+ *  This library is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU Library General Public License as
+ *  published by the Free Software Foundation; either version 2 of the
+ *  License, or (at your option) any later version. \n\n
+ *  This library is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU Library General Public License for more details. \n\n
+ *  You should have received a copy of the GNU Library General Public
+ *  License along with this program (see the file COPYING.LIB for more
+ *  details); if not, write to the Free Software Foundation, Inc.,
+ *  675 Mass Ave, Cambridge, MA 02139, USA.
+ */
+
+#ifndef GBLMEASUREMENT_H_
+#define GBLMEASUREMENT_H_
+
+#include<iostream>
+#include<vector>
+#include <array>
+#include<math.h>
+#include "VMatrix.h"
+
+#ifdef GBL_EIGEN_SUPPORT_ROOT
+#include "TVectorD.h"
+#include "TMatrixD.h"
+#include "TMatrixDSym.h"
+#include "TMatrixDSymEigen.h"
+#endif
+
+#include "Eigen/Dense"
+
+namespace gbl {
+
+typedef Eigen::Matrix<double, 5, 1> Vector5d;
+typedef Eigen::Matrix<double, 2, 3> Matrix23d;
+typedef Eigen::Matrix<double, 2, 5> Matrix25d;
+typedef Eigen::Matrix<double, 2, 7> Matrix27d;
+typedef Eigen::Matrix<double, 3, 2> Matrix32d;
+typedef Eigen::Matrix<double, 5, 5> Matrix5d;
+
+/// Measurement at point.
+/**
+ * User supplied measurement at point on (initial) trajectory.
+ *
+ * Must have measurement (1D - 5D). May have:
+ *
+ *   -# Additional local parameters (with derivatives). Fitted together with track parameters.
+ *   -# Additional global parameters (with labels and derivatives). Not fitted, only passed
+ *      on to (binary) file for fitting with Millepede-II.
+ */
+class GblMeasurement {
+public:
+	GblMeasurement(const Eigen::MatrixXd &aProjection,
+			const Eigen::VectorXd &aResiduals,
+			const Eigen::MatrixXd &aPrecision, double minPrecision = 0.);
+	GblMeasurement(const Eigen::VectorXd &aResiduals,
+			const Eigen::MatrixXd &aPrecision, double minPrecision = 0.);
+	GblMeasurement(const GblMeasurement&) = default;
+	GblMeasurement& operator=(const GblMeasurement&) = default;
+	GblMeasurement(GblMeasurement&&) = default;
+	GblMeasurement& operator=(GblMeasurement&&) = default;
+	virtual ~GblMeasurement();
+#ifdef GBL_EIGEN_SUPPORT_ROOT
+	// input via ROOT
+	GblMeasurement(const TMatrixD &aProjection, const TVectorD &aResiduals,
+			const TVectorD &aPrecision, double minPrecision = 0.);
+	GblMeasurement(const TMatrixD &aProjection, const TVectorD &aResiduals,
+			const TMatrixDSym &aPrecision, double minPrecision = 0.);
+	GblMeasurement(const TVectorD &aResiduals, const TVectorD &aPrecision,
+			double minPrecision = 0.);
+	GblMeasurement(const TVectorD &aResiduals, const TMatrixDSym &aPrecision,
+			double minPrecision = 0.);
+	void addLocals(const TMatrixD &aDerivatives);
+	void addGlobals(const std::vector<int> &aLabels,
+			const TMatrixD &aDerivatives);
+#endif
+
+	// input via Eigen
+	void addLocals(const Eigen::MatrixXd &aDerivatives);
+	void addGlobals(const std::vector<int> &aLabels,
+			const Eigen::MatrixXd &aDerivatives);
+	void setEnabled(bool flag);
+	bool isEnabled() const;
+	unsigned int getMeasDim() const;
+	double getMeasPrecMin() const;
+	void getMeasurement(Matrix5d &aProjection, Vector5d &aResiduals,
+			Vector5d &aPrecision) const;
+	void getMeasTransformation(Eigen::MatrixXd &aTransformation) const;
+	unsigned int getNumLocals() const;
+	const Eigen::MatrixXd& getLocalDerivatives() const;
+	unsigned int getNumGlobals() const;
+	void printMeasurement(unsigned int level = 0) const;
+	void getGlobalLabels(std::vector<int> &aLabels) const;
+	void getGlobalDerivatives(Eigen::MatrixXd &aDerivatives) const;
+	void getGlobalLabelsAndDerivatives(unsigned int aRow,
+			std::vector<int> &aLabels, std::vector<double> &aDerivatives) const;
+
+private:
+	bool enabled; ///< Enabled flag (to be used)
+	unsigned int measDim; ///< Dimension of measurement (1-5), 0 indicates absence of measurement
+	double measPrecMin; ///< Minimal measurement precision (for usage)
+	Matrix5d measProjection; ///< Projection from measurement to local system
+
+	Vector5d measResiduals; ///< Measurement residuals
+	Vector5d measPrecision; ///< Measurement precision (diagonal of inverse covariance matrix)
+	bool transFlag; ///< Transformation exists?
+	Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic,
+			Eigen::ColMajor /* default */, 5, 5> measTransformation; ///< Transformation of diagonalization (of meas. precision matrix)
+	Eigen::MatrixXd localDerivatives; ///< Derivatives of measurement vs additional local (fit) parameters
+	std::vector<int> globalLabels; ///< Labels of global (MP-II) derivatives
+	Eigen::MatrixXd globalDerivatives; ///< Derivatives of measurement vs additional global (MP-II) parameters
+};
+}
+
+#endif /* GBLMEASUREMENT_H_ */

cpp/include/GblPoint.h

@@ -13,7 +13,7 @@
  *
  *
  *  \copyright
- *  Copyright (c) 2011 - 2017 Deutsches Elektronen-Synchroton,
+ *  Copyright (c) 2011 - 2021 Deutsches Elektronen-Synchroton,
  *  Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
  *  This library is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU Library General Public License as
@@ -32,16 +32,12 @@
 #ifndef GBLPOINT_H_
 #define GBLPOINT_H_
 
+#include "GblMeasurement.h"
+
 #include<iostream>
 #include<vector>
 #include<math.h>
 #include <stdexcept>
-#ifdef GBL_EIGEN_SUPPORT_ROOT
-#include "TVectorD.h"
-#include "TMatrixD.h"
-#include "TMatrixDSym.h"
-#include "TMatrixDSymEigen.h"
-#endif
 
 #include "Eigen/Dense"
 
@@ -60,11 +56,8 @@ typedef Eigen::Matrix<double, 5, 5> Matrix5d;
  *
  * Must have jacobian for propagation from previous point. May have:
  *
- *   -# Measurement (1D - 5D)
+ *   -# Measurement(s) (1D - 5D)
  *   -# Scatterer (thin, 2D kinks)
- *   -# Additional local parameters (with derivatives). Fitted together with track parameters.
- *   -# Additional global parameters (with labels and derivatives). Not fitted, only passed
- *      on to (binary) file for fitting with Millepede-II.
  */
 class GblPoint {
 public:
@@ -110,9 +103,9 @@ public:
 	template<typename Projection, typename Residuals, typename Precision,
 			typename std::enable_if<(Precision::ColsAtCompileTime != 1)>::type* =
 					nullptr>
-	void addMeasurement(const Eigen::MatrixBase<Projection>& aProjection,
-			const Eigen::MatrixBase<Residuals>& aResiduals,
-			const Eigen::MatrixBase<Precision>& aPrecision,
+	void addMeasurement(const Eigen::MatrixBase<Projection> &aProjection,
+			const Eigen::MatrixBase<Residuals> &aResiduals,
+			const Eigen::MatrixBase<Precision> &aPrecision,
 			double minPrecision = 0.);
 
 	/// Add a measurement to a point.
@@ -130,9 +123,9 @@ public:
 	template<typename Projection, typename Residuals, typename Precision,
 			typename std::enable_if<(Precision::ColsAtCompileTime == 1)>::type* =
 					nullptr>
-	void addMeasurement(const Eigen::MatrixBase<Projection>& aProjection,
-			const Eigen::MatrixBase<Residuals>& aResiduals,
-			const Eigen::MatrixBase<Precision>& aPrecision,
+	void addMeasurement(const Eigen::MatrixBase<Projection> &aProjection,
+			const Eigen::MatrixBase<Residuals> &aResiduals,
+			const Eigen::MatrixBase<Precision> &aPrecision,
 			double minPrecision = 0.);
 
 	/// Add a measurement to a point.
@@ -148,8 +141,8 @@ public:
 	 */
 	template<typename Residuals, typename Precision, typename std::enable_if<
 			(Precision::ColsAtCompileTime != 1)>::type* = nullptr>
-	void addMeasurement(const Eigen::MatrixBase<Residuals>& aResiduals,
-			const Eigen::MatrixBase<Precision>& aPrecision,
+	void addMeasurement(const Eigen::MatrixBase<Residuals> &aResiduals,
+			const Eigen::MatrixBase<Precision> &aPrecision,
 			double minPrecision = 0.);
 
 	/// Add a measurement to a point.
@@ -164,8 +157,8 @@ public:
 	 */
 	template<typename Residuals, typename Precision, typename std::enable_if<
 			(Precision::ColsAtCompileTime == 1)>::type* = nullptr>
-	void addMeasurement(const Eigen::MatrixBase<Residuals>& aResiduals,
-			const Eigen::MatrixBase<Precision>& aPrecision,
+	void addMeasurement(const Eigen::MatrixBase<Residuals> &aResiduals,
+			const Eigen::MatrixBase<Precision> &aPrecision,
 			double minPrecision = 0.);
 
 	/// Add a (thin) scatterer to a point.
@@ -188,7 +181,7 @@ public:
 	template<typename Precision, typename std::enable_if<
 			(Precision::ColsAtCompileTime == 2)>::type* = nullptr>
 	void addScatterer(const Eigen::Vector2d &aResiduals,
-			const Eigen::MatrixBase<Precision>& aPrecision);
+			const Eigen::MatrixBase<Precision> &aPrecision);
 
 	/// Add a (thin) scatterer to a point.
 	/**
@@ -210,7 +203,7 @@ public:
 	template<typename Precision, typename std::enable_if<
 			(Precision::ColsAtCompileTime == 1)>::type* = nullptr>
 	void addScatterer(const Eigen::Vector2d &aResiduals,
-			const Eigen::MatrixBase<Precision>& aPrecision);
+			const Eigen::MatrixBase<Precision> &aPrecision);
 
 	/// Add local derivatives to a point.
 	/**
@@ -219,7 +212,7 @@ public:
 	 * \param [in] aDerivatives Local derivatives (matrix)
 	 */
 	template<typename Derivative>
-	void addLocals(const Eigen::MatrixBase<Derivative>& aDerivatives);
+	void addLocals(const Eigen::MatrixBase<Derivative> &aDerivatives);
 	template<typename Derivative>
 
 	/// Add global derivatives to a point.
@@ -230,30 +223,23 @@ public:
 	 * \param [in] aDerivatives Global derivatives (matrix)
 	 */
 	void addGlobals(const std::vector<int> &aLabels,
-			const Eigen::MatrixBase<Derivative>& aDerivatives);
+			const Eigen::MatrixBase<Derivative> &aDerivatives);
 	//
-	unsigned int hasMeasurement() const;
-	double getMeasPrecMin() const;
-	void getMeasurement(Matrix5d &aProjection, Vector5d &aResiduals,
-			Vector5d &aPrecision) const;
-	void getMeasTransformation(Eigen::MatrixXd &aTransformation) const;
+	unsigned int numMeasurements() const;
 	bool hasScatterer() const;
 	void getScatterer(Eigen::Matrix2d &aTransformation,
 			Eigen::Vector2d &aResiduals, Eigen::Vector2d &aPrecision) const;
 	void getScatTransformation(Eigen::Matrix2d &aTransformation) const;
-	unsigned int getNumLocals() const;
-	const Eigen::MatrixXd& getLocalDerivatives() const;
-	unsigned int getNumGlobals() const;
-	void getGlobalLabels(std::vector<int> &aLabels) const;
-	void getGlobalDerivatives(Eigen::MatrixXd &aDerivatives) const;
-	void getGlobalLabelsAndDerivatives(unsigned int aRow,
-			std::vector<int> &aLabels, std::vector<double> &aDerivatives) const;
 	unsigned int getLabel() const;
 	int getOffset() const;
 	const Matrix5d& getP2pJacobian() const;
 	void getDerivatives(int aDirection, Eigen::Matrix2d &matW,
 			Eigen::Matrix2d &matWJ, Eigen::Vector2d &vecWd) const;
 	void printPoint(unsigned int level = 0) const;
+	std::vector<GblMeasurement>::iterator getMeasBegin();
+	std::vector<GblMeasurement>::iterator getMeasEnd();
+	void getGlobalLabelsAndDerivatives(unsigned int aMeas, unsigned int aRow,
+			std::vector<int> &aLabels, std::vector<double> &aDerivatives) const;
 
 private:
 	friend class GblTrajectory; // to have the following setters private
@@ -267,104 +253,68 @@ private:
 	Matrix5d p2pJacobian; ///< Point-to-point jacobian from previous point
 	Matrix5d prevJacobian; ///< Jacobian to previous scatterer (or first measurement)
 	Matrix5d nextJacobian; ///< Jacobian to next scatterer (or last measurement)
-	unsigned int measDim; ///< Dimension of measurement (1-5), 0 indicates absence of measurement
-	double measPrecMin; ///< Minimal measurement precision (for usage)
-	Matrix5d measProjection; ///< Projection from measurement to local system
-
-	Vector5d measResiduals; ///< Measurement residuals
-	Vector5d measPrecision; ///< Measurement precision (diagonal of inverse covariance matrix)
-	bool transFlag; ///< Transformation exists?
-	Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic,
-			Eigen::ColMajor /* default */, 5, 5> measTransformation; ///< Transformation of diagonalization (of meas. precision matrix)
 	bool scatFlag; ///< Scatterer present?
 	Eigen::Matrix2d scatTransformation; ///< Transformation of diagonalization (of scat. precision matrix)
 	Eigen::Vector2d scatResiduals; ///< Scattering residuals (initial kinks if iterating)
 	Eigen::Vector2d scatPrecision; ///< Scattering precision (diagonal of inverse covariance matrix)
-	Eigen::MatrixXd localDerivatives; ///< Derivatives of measurement vs additional local (fit) parameters
-	std::vector<int> globalLabels; ///< Labels of global (MP-II) derivatives
-	Eigen::MatrixXd globalDerivatives; ///< Derivatives of measurement vs additional global (MP-II) parameters
+	std::vector<GblMeasurement> theMeasurements; ///< List of measurements at point
+
 };
 
 template<typename Projection, typename Residuals, typename Precision,
 		typename std::enable_if<(Precision::ColsAtCompileTime != 1)>::type*>
-void GblPoint::addMeasurement(const Eigen::MatrixBase<Projection>& aProjection,
-		const Eigen::MatrixBase<Residuals>& aResiduals,
-		const Eigen::MatrixBase<Precision>& aPrecision, double minPrecision) {
+void GblPoint::addMeasurement(const Eigen::MatrixBase<Projection> &aProjection,
+		const Eigen::MatrixBase<Residuals> &aResiduals,
+		const Eigen::MatrixBase<Precision> &aPrecision, double minPrecision) {
 	static_assert(static_cast<int>(Residuals::ColsAtCompileTime) == 1, "addMeasurement: cols(Residuals) must be 1 (vector)");
 	static_assert(static_cast<int>(Residuals::RowsAtCompileTime) <= 5 or static_cast<int>(Residuals::RowsAtCompileTime) == Eigen::Dynamic, "addMeasurement: rows(Residuals) must be 1-5 or dynamic");
 	static_assert(static_cast<int>(Residuals::RowsAtCompileTime) == static_cast<int>(Precision::RowsAtCompileTime), "addMeasurement: rows(Residuals) and rows(Precision) must be equal");
 	static_assert(static_cast<int>(Residuals::RowsAtCompileTime) == static_cast<int>(Projection::RowsAtCompileTime), "addMeasurement: rows(Residuals) and rows(Projection) must be equal");
 	static_assert(static_cast<int>(Precision::RowsAtCompileTime) == static_cast<int>(Precision::ColsAtCompileTime), "addMeasurement: rows(Precision) and cols(Precision) must be equal");
 	static_assert(static_cast<int>(Projection::RowsAtCompileTime) == static_cast<int>(Projection::ColsAtCompileTime), "addMeasurement: rows(Projection) and cols(Projection) must be equal");
-	measDim = aResiduals.rows();
-	measPrecMin = minPrecision;
-	// arbitrary precision matrix
-	Eigen::SelfAdjointEigenSolver<typename Precision::PlainObject> measEigen {
-			aPrecision };
-	measTransformation = measEigen.eigenvectors().transpose();
-	transFlag = true;
-	measResiduals.tail(measDim) = measTransformation * aResiduals;
-	measPrecision.tail(measDim) = measEigen.eigenvalues();
-	measProjection.bottomRightCorner(measDim, measDim) = measTransformation
-			* aProjection;
+	theMeasurements.emplace_back(aProjection, aResiduals, aPrecision,
+			minPrecision);
 }
 
 template<typename Projection, typename Residuals, typename Precision,
 		typename std::enable_if<(Precision::ColsAtCompileTime == 1)>::type*>
-void GblPoint::addMeasurement(const Eigen::MatrixBase<Projection>& aProjection,
-		const Eigen::MatrixBase<Residuals>& aResiduals,
-		const Eigen::MatrixBase<Precision>& aPrecision, double minPrecision) {
+void GblPoint::addMeasurement(const Eigen::MatrixBase<Projection> &aProjection,
+		const Eigen::MatrixBase<Residuals> &aResiduals,
+		const Eigen::MatrixBase<Precision> &aPrecision, double minPrecision) {
 	static_assert(static_cast<int>(Residuals::ColsAtCompileTime) == 1, "addMeasurement: cols(Residuals) must be 1 (vector)");
 	static_assert(static_cast<int>(Residuals::RowsAtCompileTime) <= 5 or static_cast<int>(Residuals::RowsAtCompileTime) == Eigen::Dynamic, "addMeasurement: rows(Residuals) must be 1-5 or dynamic");
 	static_assert(static_cast<int>(Residuals::RowsAtCompileTime) == static_cast<int>(Precision::RowsAtCompileTime), "addMeasurement: rows(Residuals) and rows(Precision) must be equal");
 	static_assert(static_cast<int>(Residuals::RowsAtCompileTime) == static_cast<int>(Projection::RowsAtCompileTime), "addMeasurement: rows(Residuals) and rows(Projection) must be equal");
 	static_assert(static_cast<int>(Projection::RowsAtCompileTime) == static_cast<int>(Projection::ColsAtCompileTime), "addMeasurement: rows(Projection) and cols(Projection) must be equal");
-	measDim = aResiduals.rows();
-	measPrecMin = minPrecision;
-	// diagonal precision matrix
-	measResiduals.tail(measDim) = aResiduals;
-	measPrecision.tail(measDim) = aPrecision;
-	measProjection.bottomRightCorner(measDim, measDim) = aProjection;
+	//theMeasurements.push_back(GblMeasurement(aProjection, aResiduals, aPrecision, minPrecision));
+	theMeasurements.emplace_back(aProjection, aResiduals, aPrecision,
+			minPrecision);
 }
 
 template<typename Residuals, typename Precision, typename std::enable_if<
 		(Precision::ColsAtCompileTime != 1)>::type*>
-void GblPoint::addMeasurement(const Eigen::MatrixBase<Residuals>& aResiduals,
-		const Eigen::MatrixBase<Precision>& aPrecision, double minPrecision) {
+void GblPoint::addMeasurement(const Eigen::MatrixBase<Residuals> &aResiduals,
+		const Eigen::MatrixBase<Precision> &aPrecision, double minPrecision) {
 	static_assert(static_cast<int>(Residuals::ColsAtCompileTime) == 1, "addMeasurement: cols(Residuals) must be 1 (vector)");
 	static_assert(static_cast<int>(Residuals::RowsAtCompileTime) <= 5 or static_cast<int>(Residuals::RowsAtCompileTime) == Eigen::Dynamic, "addMeasurement: rows(Residuals) must be 1-5 or dynamic");
 	static_assert(static_cast<int>(Residuals::RowsAtCompileTime) == static_cast<int>(Precision::RowsAtCompileTime), "addMeasurement: rows(Residuals) and rows(Precision) must be equal");
-	measDim = aResiduals.rows();
-	measPrecMin = minPrecision;
-	// arbitrary precision matrix
-	Eigen::SelfAdjointEigenSolver<typename Precision::PlainObject> measEigen {
-			aPrecision };
-	measTransformation = measEigen.eigenvectors().transpose();
-	transFlag = true;
-	measResiduals.tail(measDim) = measTransformation * aResiduals;
-	measPrecision.tail(measDim) = measEigen.eigenvalues();
-	measProjection.bottomRightCorner(measDim, measDim) = measTransformation;
+	theMeasurements.emplace_back(aResiduals, aPrecision, minPrecision);
 }
 
 template<typename Residuals, typename Precision, typename std::enable_if<
 		(Precision::ColsAtCompileTime == 1)>::type*>
-void GblPoint::addMeasurement(const Eigen::MatrixBase<Residuals>& aResiduals,
-		const Eigen::MatrixBase<Precision>& aPrecision, double minPrecision) {
+void GblPoint::addMeasurement(const Eigen::MatrixBase<Residuals> &aResiduals,
+		const Eigen::MatrixBase<Precision> &aPrecision, double minPrecision) {
 	static_assert(static_cast<int>(Residuals::ColsAtCompileTime) == 1, "addMeasurement: cols(Residuals) must be 1 (vector)");
 	static_assert(static_cast<int>(Residuals::RowsAtCompileTime) <= 5 or static_cast<int>(Residuals::RowsAtCompileTime) == Eigen::Dynamic, "addMeasurement: rows(Residuals) must be 1-5 or dynamic");
 	static_assert(static_cast<int>(Residuals::RowsAtCompileTime) == static_cast<int>(Precision::RowsAtCompileTime), "addMeasurement: rows(Residuals) and rows(Precision) must be equal");
-	measDim = aResiduals.rows();
-	measPrecMin = minPrecision;
-	// diagonal precision matrix
-	measResiduals.tail(measDim) = aResiduals;
-	measPrecision.tail(measDim) = aPrecision;
-	measProjection.setIdentity();
+	theMeasurements.emplace_back(aResiduals, aPrecision, minPrecision);
 }
 
 template<typename Precision, typename std::enable_if<
 		(Precision::ColsAtCompileTime == 2)>::type*>
 void GblPoint::addScatterer(const Eigen::Vector2d &aResiduals,
-		const Eigen::MatrixBase<Precision>& aPrecision) {
+		const Eigen::MatrixBase<Precision> &aPrecision) {
 	static_assert(static_cast<int>(Precision::RowsAtCompileTime) == 2 or static_cast<int>(Precision::RowsAtCompileTime) == Eigen::Dynamic, "addScatterer: rows(Precision) must be 2 or dynamic");
 	scatFlag = true;
 	// arbitrary precision matrix
@@ -379,7 +329,7 @@ void GblPoint::addScatterer(const Eigen::Vector2d &aResiduals,
 template<typename Precision, typename std::enable_if<
 		(Precision::ColsAtCompileTime == 1)>::type*>
 void GblPoint::addScatterer(const Eigen::Vector2d &aResiduals,
-		const Eigen::MatrixBase<Precision>& aPrecision) {
+		const Eigen::MatrixBase<Precision> &aPrecision) {
 	static_assert(static_cast<int>(Precision::RowsAtCompileTime) == 2 or static_cast<int>(Precision::RowsAtCompileTime) == Eigen::Dynamic, "addScatterer: rows(Precision) must be 2 or dynamic");
 	scatFlag = true;
 	scatResiduals = aResiduals;
@@ -388,30 +338,16 @@ void GblPoint::addScatterer(const Eigen::Vector2d &aResiduals,
 }
 
 template<typename Derivative>
-void GblPoint::addLocals(const Eigen::MatrixBase<Derivative>& aDerivatives) {
-	if (measDim) {
-		localDerivatives.resize(aDerivatives.rows(), aDerivatives.cols());
-		if (transFlag) {
-			localDerivatives = measTransformation * aDerivatives;
-		} else {
-			localDerivatives = aDerivatives;
-		}
-	}
+void GblPoint::addLocals(const Eigen::MatrixBase<Derivative> &aDerivatives) {
+	if (theMeasurements.size())
+		theMeasurements.back().addLocals(aDerivatives);
 }
 
 template<typename Derivative>
 void GblPoint::addGlobals(const std::vector<int> &aLabels,
-		const Eigen::MatrixBase<Derivative>& aDerivatives) {
-	if (measDim) {
-		globalLabels = aLabels;
-		globalDerivatives.resize(aDerivatives.rows(), aDerivatives.cols());
-		if (transFlag) {
-			globalDerivatives = measTransformation * aDerivatives;
-		} else {
-			globalDerivatives = aDerivatives;
-		}
-
-	}
+		const Eigen::MatrixBase<Derivative> &aDerivatives) {
+	if (theMeasurements.size())
+		theMeasurements.back().addGlobals(aLabels, aDerivatives);
 }
 
 }

cpp/include/GblTrajectory.h

@@ -13,7 +13,7 @@
  *
  *
  *  \copyright
- *  Copyright (c) 2011 - 2018 Deutsches Elektronen-Synchroton,
+ *  Copyright (c) 2011 - 2021 Deutsches Elektronen-Synchroton,
  *  Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
  *  This library is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU Library General Public License as
@@ -52,7 +52,7 @@ public:
 	GblTrajectory(const std::vector<GblPoint> &aPointList, bool flagCurv = true,
 			bool flagU1dir = true, bool flagU2dir = true);
 	GblTrajectory(
-			const std::vector<std::pair<std::vector<GblPoint>, Eigen::MatrixXd> >& aPointsAndTransList);
+			const std::vector<std::pair<std::vector<GblPoint>, Eigen::MatrixXd> > &aPointsAndTransList);
 
 	/// Create new (simple) trajectory from list of points with external seed.
 	/**
@@ -69,7 +69,7 @@ public:
 	 */
 	template<typename Seed>
 	GblTrajectory(const std::vector<GblPoint> &aPointList, unsigned int aLabel,
-			const Eigen::MatrixBase<Seed>& aSeed, bool flagCurv = true,
+			const Eigen::MatrixBase<Seed> &aSeed, bool flagCurv = true,
 			bool flagU1dir = true, bool flagU2dir = true);
 
 	/// Create new composed trajectory from list of points and transformations with arbitrary external measurements.
@@ -88,10 +88,10 @@ public:
 			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);