Commit f2c3b8e3 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

cpp: diagonalization of scatterers implemented

git-svn-id: http://svnsrv.desy.de/public/GeneralBrokenLines/trunk@93 281f6f2b-e318-4fd1-8bce-1a4ba7aab212
parent 4ae9515c
......@@ -42,7 +42,8 @@ void example1() {
*
* This example simulates and refits tracks in a system of planar detectors
* with 2D measurements in a constant magnet field in Z direction using
* the curvilinear system as local system. The true track parameters are
* the curvilinear system as local system and (Q/P, slopes, offsets) as
* local track parameters. The true track parameters are
* randomly smeared with respect to a (constant and straight) reference
* trajectory with direction (lambda, phi) and are used (only) for the
* on-the-fly simulation of the measurements and scatterers. The predictions
......@@ -220,6 +221,7 @@ void example1() {
// create trajectory
//GblTrajectory traj(listOfPoints);
GblTrajectory traj(listOfPoints, seedLabel, clSeed); // with external seed
//traj.printPoints();
if (not traj.isValid()) {
std::cout << " Invalid GblTrajectory -> skip" << std::endl;
continue;
......@@ -237,6 +239,15 @@ void example1() {
aCorrection.Print();
std::cout << " cov " << std::endl;
aCovariance.Print(); */
/* look at residuals
for (unsigned int label=1; label<=listOfPoints.size(); ++label) {
unsigned int numData=0;
std::cout << " measResults, label " << label << std::endl;
TVectorD residuals(2), measErr(2), resErr(2), downWeights(2);
traj.getMeasResults(label, numData, residuals, measErr, resErr, downWeights);
std::cout << " measResults, numData " << numData << std::endl;
// residuals.Print(); measErr.Print(); resErr.Print();
} */
// debug printout
//traj.printTrajectory();
//traj.printPoints();
......
......@@ -52,14 +52,21 @@ public:
const TVectorD &aPrecision, double minPrecision = 0.);
void addMeasurement(const TMatrixD &aProjection, const TVectorD &aResiduals,
const TMatrixDSym &aPrecision, double minPrecision = 0.);
void addMeasurement(const TVectorD &aResiduals, const TVectorD &aPrecision,
double minPrecision = 0.);
void addMeasurement(const TVectorD &aResiduals,
const TMatrixDSym &aPrecision, double minPrecision = 0.);
unsigned int hasMeasurement() const;
void getMeasurement(SMatrix55 &aProjection, SVector5 &aResiduals,
SVector5 &aPrecision) const;
void getMeasTransformation(TMatrixD &aTransformation) const;
void addScatterer(const TVectorD &aResiduals, const TVectorD &aPrecision);
void addScatterer(const TVectorD &aResiduals,
const TMatrixDSym &aPrecision);
bool hasScatterer() const;
void getScatterer(SVector2 &aResiduals, SVector2 &aPrecision) const;
void getScatterer(SMatrix22 &aTransformation, SVector2 &aResiduals,
SVector2 &aPrecision) const;
void getScatTransformation(TMatrixD &aTransformation) const;
void addLocals(const TMatrixD &aDerivatives);
unsigned int getNumLocals() const;
const TMatrixD& getLocalDerivatives() const;
......@@ -90,8 +97,9 @@ private:
SVector5 measResiduals; ///< Measurement residuals
SVector5 measPrecision; ///< Measurement precision (diagonal of inverse covariance matrix)
bool transFlag; ///< Transformation exists?
TMatrixD measTransformation; ///< Transformation of diagonalization (of precision matrix)
TMatrixD measTransformation; ///< Transformation of diagonalization (of meas. precision matrix)
bool scatFlag; ///< Scatterer present?
SMatrix22 scatTransformation; ///< Transformation of diagonalization (of scat. precision matrix)
SVector2 scatResiduals; ///< Scattering residuals (initial kinks if iterating)
SVector2 scatPrecision; ///< Scattering precision (diagonal of inverse covariance matrix)
TMatrixD localDerivatives; ///< Derivatives of measurement vs additional local (fit) parameters
......
......@@ -128,7 +128,7 @@ void GblData::addDerivatives(const std::vector<unsigned int> &index,
}
}
/// Calculate prediction for data from fit.
/// Calculate prediction for data from fit (by GblTrajectory::fit).
void GblData::setPrediction(const VVector &aVector) {
thePrediction = 0.;
......@@ -137,7 +137,7 @@ void GblData::setPrediction(const VVector &aVector) {
}
}
/// Outlier down weighting with M-estimators.
/// Outlier down weighting with M-estimators (by GblTrajectory::fit).
/**
* \param [in] aMethod M-estimator (1: Tukey, 2:Huber, 3:Cauchy)
*/
......
......@@ -35,11 +35,11 @@ GblPoint::GblPoint(const SMatrix55 &aJacobian) :
GblPoint::~GblPoint() {
}
/// Add a mesurement to a point.
/// Add a measurement to a point.
/**
* Add measurement with diagonal precision (inverse covariance) matrix.
* Add measurement (in meas. system) with diagonal precision (inverse covariance) matrix.
* ((up to) 2D: position, 4D: slope+position, 5D: curvature+slope+position)
* \param [in] aProjection Projection from measurement to local system
* \param [in] aProjection Projection from local to measurement system
* \param [in] aResiduals Measurement residuals
* \param [in] aPrecision Measurement precision (diagonal)
* \param [in] minPrecision Minimal precision to accept measurement
......@@ -59,12 +59,12 @@ void GblPoint::addMeasurement(const TMatrixD &aProjection,
}
}
/// Add a mesurement to a point.
/// Add a measurement to a point.
/**
* Add measurement on local system with arbitrary precision (inverse covariance) matrix.
* Add measurement (in meas. system) with arbitrary precision (inverse covariance) matrix.
* Will be diagonalized.
* ((up to) 2D: position, 4D: slope+position, 5D: curvature+slope+position)
* \param [in] aProjection Projection from measurement to local system
* \param [in] aProjection Projection from local to measurement system
* \param [in] aResiduals Measurement residuals
* \param [in] aPrecision Measurement precision (matrix)
* \param [in] minPrecision Minimal precision to accept measurement
......@@ -92,9 +92,29 @@ void GblPoint::addMeasurement(const TMatrixD &aProjection,
}
}
/// Add a mesurement to a point.
/// Add a measurement to a point.
/**
* Add measurement on local system with arbitrary precision (inverse covariance) matrix.
* Add measurement in local system with diagonal precision (inverse covariance) matrix.
* ((up to) 2D: position, 4D: slope+position, 5D: curvature+slope+position)
* \param [in] aResiduals Measurement residuals
* \param [in] aPrecision Measurement precision (diagonal)
* \param [in] minPrecision Minimal precision to accept measurement
*/
void GblPoint::addMeasurement(const TVectorD &aResiduals,
const TVectorD &aPrecision, double minPrecision) {
measDim = aResiduals.GetNrows();
unsigned int iOff = 5 - measDim;
for (unsigned int i = 0; i < measDim; ++i) {
measResiduals(iOff + i) = aResiduals[i];
measPrecision(iOff + i) = (
aPrecision[i] >= minPrecision ? aPrecision[i] : 0.);
}
measProjection = ROOT::Math::SMatrixIdentity();
}
/// Add a measurement to a point.
/**
* Add measurement in local system with arbitrary precision (inverse covariance) matrix.
* Will be diagonalized.
* ((up to) 2D: position, 4D: slope+position, 5D: curvature+slope+position)
* \param [in] aResiduals Measurement residuals
......@@ -144,8 +164,22 @@ void GblPoint::getMeasurement(SMatrix55 &aProjection, SVector5 &aResiduals,
aPrecision = measPrecision;
}
/// Get measurement transformation (from diagonalization).
/**
* \param [out] aTransformation Transformation matrix
*/
void GblPoint::getMeasTransformation(TMatrixD &aTransformation) const {
aTransformation.ResizeTo(measDim, measDim);
if (transFlag) {
aTransformation = measTransformation;
} else {
aTransformation.UnitMatrix();
}
}
/// Add a (thin) scatterer to a point.
/**
* Add scatterer with diagonal precision (inverse covariance) matrix.
* Changes local track direction.
*
* \param [in] aResiduals Scatterer residuals
......@@ -158,6 +192,40 @@ void GblPoint::addScatterer(const TVectorD &aResiduals,
scatResiduals(1) = aResiduals[1];
scatPrecision(0) = aPrecision[0];
scatPrecision(1) = aPrecision[1];
scatTransformation = ROOT::Math::SMatrixIdentity();
}
/// Add a (thin) scatterer to a point.
/**
* Add scatterer with arbitrary precision (inverse covariance) matrix.
* Will be diagonalized. Changes local track direction.
*
* The precision matrix for the local slopes is defined by the
* angular scattering error theta_0 and the scalar products c_1, c_2 of the
* offset directions in the local frame with the track direction:
*
* (1 - c_1*c_1 - c_2*c_2) | 1 - c_1*c_1 - c_1*c_2 |
* P = ~~~~~~~~~~~~~~~~~~~~~~~ * | |
* theta_0*theta_0 | - c_1*c_2 1 - c_2*c_2 |
*
* \param [in] aResiduals Scatterer residuals
* \param [in] aPrecision Scatterer precision (matrix)
*/
void GblPoint::addScatterer(const TVectorD &aResiduals,
const TMatrixDSym &aPrecision) {
scatFlag = true;
TMatrixDSymEigen scatEigen(aPrecision);
TMatrixD aTransformation = scatEigen.GetEigenVectors();
aTransformation.T();
TVectorD transResiduals = aTransformation * aResiduals;
TVectorD transPrecision = scatEigen.GetEigenValues();
for (unsigned int i = 0; i < 2; ++i) {
scatResiduals(i) = transResiduals[i];
scatPrecision(i) = transPrecision[i];
for (unsigned int j = 0; j < 2; ++j) {
scatTransformation(i, j) = aTransformation(i, j);
}
}
}
/// Check for scatterer at a point.
......@@ -167,14 +235,34 @@ bool GblPoint::hasScatterer() const {
/// Retrieve scatterer of a point.
/**
* \param [out] aTransformation Scatterer transformation from diagonalization
* \param [out] aResiduals Scatterer residuals
* \param [out] aPrecision Scatterer precision (diagonal)
*/
void GblPoint::getScatterer(SVector2 &aResiduals, SVector2 &aPrecision) const {
void GblPoint::getScatterer(SMatrix22 &aTransformation, SVector2 &aResiduals,
SVector2 &aPrecision) const {
aTransformation = scatTransformation;
aResiduals = scatResiduals;
aPrecision = scatPrecision;
}
/// Get scatterer transformation (from diagonalization).
/**
* \param [out] aTransformation Transformation matrix
*/
void GblPoint::getScatTransformation(TMatrixD &aTransformation) const {
aTransformation.ResizeTo(2, 2);
if (scatFlag) {
for (unsigned int i = 0; i < 2; ++i) {
for (unsigned int j = 0; j < 2; ++j) {
aTransformation(i, j) = scatTransformation(i, j);
}
}
} else {
aTransformation.UnitMatrix();
}
}
/// Add local derivatives to a point.
/**
* Point needs to have a measurement.
......@@ -236,7 +324,7 @@ const TMatrixD& GblPoint::getGlobalDerivatives() const {
return globalDerivatives;
}
/// Define label of point
/// Define label of point (by GBLTrajectory constructor)
/**
* \param [in] aLabel Label identifying point
*/
......@@ -249,7 +337,7 @@ unsigned int GblPoint::getLabel() const {
return theLabel;
}
/// Define offset for point
/// Define offset for point (by GBLTrajectory constructor)
/**
* \param [in] anOffset Offset number
*/
......@@ -267,7 +355,7 @@ const SMatrix55& GblPoint::getP2pJacobian() const {
return p2pJacobian;
}
/// Define jacobian to previous scatterer
/// Define jacobian to previous scatterer (by GBLTrajectory constructor)
/**
* \param [in] aJac Jacobian
*/
......@@ -284,7 +372,7 @@ void GblPoint::addPrevJacobian(const SMatrix55 &aJac) {
prevJacobian.Place_at(-DCAB * CA, 3, 0);
}
/// Define jacobian to next scatterer
/// Define jacobian to next scatterer (by GBLTrajectory constructor)
/**
* \param [in] aJac Jacobian
*/
......@@ -320,7 +408,9 @@ void GblPoint::getDerivatives(int aDirection, SMatrix22 &matW, SMatrix22 &matWJ,
if (!matW.InvertFast()) {
std::cout << " GblPoint::getDerivatives failed to invert matrix: "
<< matW << std::endl;
std::cout << " Possible reason for singular matrix: multiple GblPoints at same arc-length" << std::endl;
std::cout
<< " Possible reason for singular matrix: multiple GblPoints at same arc-length"
<< std::endl;
throw std::overflow_error("Singular matrix inversion exception");
}
matWJ = matW * matJ;
......
......@@ -19,6 +19,7 @@
* position, 4D: direction+position). The refit provides corrections
* to the local track parameters (in the local system) and the
* corresponding covariance matrix at any of those points.
* Non-diagonal covariance matrices will be diagonalized internally.
* Outliers can be down-weighted by use of M-estimators.
*
* The broken lines trajectory is defined by (2D) offsets at the
......@@ -81,9 +82,8 @@
* At each point on the trajectory a local coordinate system with local track
* parameters has to be defined. The first of the five parameters describes
* the bending, the next two the direction and the last two the position (offsets).
* The covariance matrix due to multiple scattering for the direction parameters
* has to be diagonal. The curvilinear system (T,U,V) with parameters
* (q/p, lambda, phi, x_t, y_t) is well suited.
* The curvilinear system (T,U,V) with parameters (q/p, lambda, phi, x_t, y_t)
* is well suited.
*
* \section impl_sec Implementation
*
......@@ -567,8 +567,9 @@ unsigned int GblTrajectory::getResults(int aSignedLabel, TVectorD &localPar,
/// Get residuals at point from measurement.
/**
* Get residual, error of measurement and residual and down-weighting
* Get (diagonalized) residual, error of measurement and residual and down-weighting
* factor for measurement at point
*
* \param [in] aLabel Label of point on trajectory
* \param [out] numData Number of data blocks from measurement at point
* \param [out] aResiduals Measurements-Predictions
......@@ -595,8 +596,9 @@ unsigned int GblTrajectory::getMeasResults(unsigned int aLabel,
/// Get (kink) residuals at point from scatterer.
/**
* Get residual, error of measurement and residual and down-weighting
* Get (diagonalized) residual, error of measurement and residual and down-weighting
* factor for scatterering kinks at point
*
* \param [in] aLabel Label of point on trajectory
* \param [out] numData Number of data blocks from scatterer at point
* \param [out] aResiduals (kink)Measurements-(kink)Predictions
......@@ -812,6 +814,7 @@ void GblTrajectory::prepare() {
}
// pseudo measurements from kinks
SMatrix22 matT;
scatDataIndex[0] = nData;
scatDataIndex[1] = nData;
// loop over trajectories
......@@ -821,12 +824,12 @@ void GblTrajectory::prepare() {
SVector2 aMeas, aPrec;
unsigned int nLabel = itPoint->getLabel();
if (itPoint->hasScatterer()) {
itPoint->getScatterer(aMeas, aPrec);
itPoint->getScatterer(matT, aMeas, aPrec);
TMatrixD transDer;
std::vector<unsigned int> labDer(7);
SMatrix27 matDer;
SMatrix27 matDer, matTDer;
getFitToKinkJacobian(labDer, matDer, *itPoint);
matTDer = matT * matDer;
if (numInnerTrans > 0) {
// transform for external parameters
TMatrixD proDer(nDim, 5);
......@@ -845,7 +848,7 @@ void GblTrajectory::prepare() {
// match
labDer[ilabel] = 0; // mark as related to external parameters
for (unsigned int k = 0; k < nDim; ++k) {
proDer(k, ifirst) = matDer(k, ilabel);
proDer(k, ifirst) = matTDer(k, ilabel);
}
}
}
......@@ -858,7 +861,7 @@ void GblTrajectory::prepare() {
unsigned int iDim = theDimension[i];
if (aPrec(iDim) > 0.) {
GblData aData(nLabel, aMeas(iDim), aPrec(iDim));
aData.addDerivatives(iDim, labDer, matDer, numLocals,
aData.addDerivatives(iDim, labDer, matTDer, numLocals,
transDer);
theData.push_back(aData);
nData++;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment