Commit 16fccb66 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

Documentation, examples improved

git-svn-id: http://svnsrv.desy.de/public/GeneralBrokenLines/trunk@90 281f6f2b-e318-4fd1-8bce-1a4ba7aab212
parent f74f3154
......@@ -13,7 +13,13 @@
using namespace gbl;
TMatrixD gblSimpleJacobian(double ds, double cosl, double bfac) {
/* Simple jacobian: quadratic in arc length difference */
/// Simple jacobian: quadratic in arc length difference
/**
* \param [in] ds (3D) arc-length
* \param [in] cosl cos(lambda)
* \param [in] bfac Bz*c
* \return jacobian
*/
TMatrixD jac(5, 5);
jac.UnitMatrix();
jac[1][0] = -bfac * ds * cosl;
......@@ -24,7 +30,8 @@ TMatrixD gblSimpleJacobian(double ds, double cosl, double bfac) {
}
void example1() {
/*
/// Simple example.
/**
* Create points on initial trajectory, create trajectory from points,
* fit and write trajectory to MP-II binary file,
* get track parameter corrections and covariance matrix at points.
......@@ -32,6 +39,15 @@ void example1() {
* Equidistant measurement layers and thin scatterers, propagation
* with simple jacobian (quadratic in arc length differences).
* Curvilinear system (U,V,T) as local coordinate system.
*
* 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
* 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
* from the reference trajectory are therefore always zero and the residuals
* needed (by addMeasurement) are equal to the measurements.
*/
//MP MilleBinary mille; // for producing MillePede-II binary file
......@@ -120,7 +136,6 @@ void example1() {
addDer[1][1] = 1.;
// arclength
double s = 0.;
std::vector<double> sPoint;
TMatrixD jacPointToPoint(5, 5);
jacPointToPoint.UnitMatrix();
// create list of points
......@@ -130,13 +145,14 @@ void example1() {
for (unsigned int iLayer = 0; iLayer < nLayer; ++iLayer) {
// std::cout << " Layer " << iLayer << ", " << s << std::endl;
// measurement directions
double sinStereo = (iLayer % 2 == 0) ? 0. : 0.5;
double sinStereo = (iLayer % 2 == 0) ? 0. : 0.1;
double cosStereo = sqrt(1.0 - sinStereo * sinStereo);
TMatrixD mDirT(3, 2);
mDirT.Zero();
mDirT[0][0] = sinStereo;
mDirT[1][0] = cosStereo;
mDirT[2][1] = 1.;
mDirT[2][0] = sinStereo;
mDirT[1][1] = -sinStereo;
mDirT[2][1] = cosStereo;
// projection measurement to local (curvilinear uv) directions (duv/dm)
TMatrixD proM2l = uvDir * mDirT;
// projection local (uv) to measurement directions (dm/duv)
......@@ -170,7 +186,6 @@ void example1() {
// add point to trajectory
listOfPoints.push_back(point);
unsigned int iLabel = listOfPoints.size();
sPoint.push_back(s);
if (iLabel == seedLabel) {
clSeed = clCov.Invert();
}
......@@ -187,7 +202,6 @@ void example1() {
point.addScatterer(scat, scatPrec);
listOfPoints.push_back(point);
iLabel = listOfPoints.size();
sPoint.push_back(s);
if (iLabel == seedLabel) {
clSeed = clCov.Invert();
}
......
......@@ -16,7 +16,7 @@
* trajectory points are defined with can describe a measurement
* or a (thin) scatterer or both. Measurements are arbitrary
* functions of the local track parameters at a point (e.g. 2D:
* position, 4D: slope+position). The refit provides corrections
* 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.
* Outliers can be down-weighted by use of M-estimators.
......@@ -77,6 +77,14 @@
* -# Optionally write trajectory to MP binary file:\n
* <tt>traj.milleOut(..)</tt>
*
* \section loc_sec Local system and local parameters
* 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.
*
* \section impl_sec Implementation
*
* Matrices are implemented with ROOT (root.cern.ch). User input or output is in the
......
......@@ -20,6 +20,15 @@ def example1():
Equidistant measurement layers and thin scatterers, propagation
with simple jacobian (quadratic in arc length differences).
Curvilinear system (U,V,T) as local coordinate system.
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
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
from the reference trajectory are therefore always zero and the residuals
needed (by addMeasurement) are equal to the measurements.
'''
def gblSimpleJacobian(ds, cosl, bfac):
'''
......@@ -86,7 +95,6 @@ def example1():
clCov[i, i] = clErr[i] ** 2
# arclength
s = 0.
sPoint = []
# point-to-point jacobian (from previous point)
jacPointToPoint = np.eye(5)
# additional (local or global) derivatives
......@@ -97,9 +105,9 @@ def example1():
for iLayer in range(nLayer):
# measurement directions
sinStereo = (0. if iLayer % 2 == 0 else 0.5)
sinStereo = (0. if iLayer % 2 == 0 else 0.1)
cosStereo = math.sqrt(1.0 - sinStereo ** 2)
mDir = np.array([[sinStereo, cosStereo, 0.0], [0., 0, 1.]])
mDir = np.array([[0., cosStereo, sinStereo], [0., -sinStereo, cosStereo]])
# projection measurement to local (curvilinear uv) directions (duv/dm)
proM2l = np.dot(uvDir, mDir.T)
# projection local (uv) to measurement directions (dm/duv)
......@@ -117,7 +125,6 @@ def example1():
addDer = -addDer # locDer flips sign every measurement
# add point to trajectory
iLabel = traj.addPoint(point)
sPoint.append(s)
if iLabel == abs(seedLabel):
clSeed = np.linalg.inv(clCov)
# propagate to scatterer
......@@ -131,7 +138,6 @@ def example1():
point = GblPoint(jacPointToPoint)
point.addScatterer([scat, scatPrec])
iLabel = traj.addPoint(point)
sPoint.append(s)
if iLabel == abs(seedLabel):
clSeed = np.linalg.inv(clCov)
......
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