Added Si tracker example, fixed down-weighted residuals

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

cpp/CMakeLists.txt

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

cpp/examples/example1.cpp

@@ -11,7 +11,7 @@
  *  \author Claus Kleinwort, DESY, 2011 (Claus.Kleinwort@desy.de)
  *
  *  \copyright
- *  Copyright (c) 2011 - 2017 Deutsches Elektronen-Synchroton,
+ *  Copyright (c) 2011 - 2018 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
@@ -75,7 +75,7 @@ double unrm() {
 }
 
 void example1() {
-	/// Simple example.
+	/// Simple technical example (curvilinear as local system).
 	/**
 	 * Create points on initial trajectory, create trajectory from points,
 	 * fit and write trajectory to MP-II binary file,

cpp/examples/example2.cpp

@@ -11,7 +11,7 @@
  *  \author Claus Kleinwort, DESY, 2011 (Claus.Kleinwort@desy.de)
  *
  *  \copyright
- *  Copyright (c) 2011 - 2017 Deutsches Elektronen-Synchroton,
+ *  Copyright (c) 2011 - 2018 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
@@ -75,7 +75,7 @@ double unrm2() {
 }
 
 void example2() {
-	/// Simple example.
+	/// Simple technical example (measurement as local system).
 	/**
 	 * Create points on initial trajectory, create trajectory from points,
 	 * fit and write trajectory to MP-II binary file,

cpp/examples/example3.cpp

@@ -11,7 +11,7 @@
  *  \author Claus Kleinwort, DESY, 2011 (Claus.Kleinwort@desy.de)
  *
  *  \copyright
- *  Copyright (c) 2011 - 2017 Deutsches Elektronen-Synchroton,
+ *  Copyright (c) 2011 - 2018 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
@@ -75,7 +75,7 @@ double unrm3() {
 }
 
 void example3() {
-	/// Simple example.
+	/// Simple technical example (measure scattering).
 	/**
 	 * Create points on initial trajectory, create trajectory from points,
 	 * fit trajectory,

cpp/examples/exampleSit.h

+/*
+ * exampleSit.h
+ *
+ *  Created on: 11 Oct 2018
+ *      Author: kleinwrt
+ */
+
+/** \file
+ *  Definitions for exampleSit.
+ *
+ *  \author Claus Kleinwort, DESY, 2011 (Claus.Kleinwort@desy.de)
+ *  \author Gregor Mittag, DESY, 2017 (templates and other optimizations)
+ *
+ *
+ *  \copyright
+ *  Copyright (c) 2011 - 2017 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 EXAMPLESIT_H_
+#define EXAMPLESIT_H_
+
+#include "Eigen/Dense"
+
+//! Namespace for the general broken lines package
+namespace gbl {
+
+double unrm();
+double unif();
+Eigen::Matrix<double, 5, 5> gblSimpleJacobian(double ds, double cosl,
+		double bfac);
+double gblMultipleScatteringError(double qbyp, double xbyx0);
+
+/// Prediction on helix
+/**
+ * Prediction at intersection of helix and measurement plane
+ */
+class GblHelixPrediction {
+public:
+	GblHelixPrediction(double sArc, const Eigen::Vector2d& aPred,
+			const Eigen::Vector3d& tDir, const Eigen::Vector3d& uDir,
+			const Eigen::Vector3d& vDir, const Eigen::Vector3d& nDir,
+			const Eigen::Vector3d& aPos);
+	virtual ~GblHelixPrediction();
+	double getArcLength() const;
+	const Eigen::Vector2d& getMeasPred() const;
+	const Eigen::Vector3d& getMeasPos() const;
+	const Eigen::Matrix3d& getMeasSystemDirs() const;
+	double getCosIncidence() const;
+	const Eigen::Matrix<double, 3, 6> getRigidBodyDerGlobal(
+			Eigen::Vector3d& rotCenter) const;
+
+private:
+	const double sarc; ///< arc-length at prediction
+	const Eigen::Vector2d pred; ///< prediction for measurement (u,v)
+	const Eigen::Vector3d tdir; ///< track direction at prediction
+	const Eigen::Vector3d udir; ///< measurement direction for u
+	const Eigen::Vector3d vdir; ///< measurement direction for v
+	const Eigen::Vector3d ndir; ///< normal to measurement plane
+	const Eigen::Vector3d pos; ///< position at prediction
+	Eigen::Matrix3d global2meas; ///< transformation into measurement system
+};
+
+///Simple helix
+/**
+ * Circle in XY plane, straight line in ZS.
+ */
+class GblSimpleHelix {
+public:
+	GblSimpleHelix(double aRinv, double aPhi0, double aDca, double aDzds,
+			double aZ0);
+	virtual ~GblSimpleHelix();
+	double getArcLengthXY(double xPos, double yPos) const;
+	void moveToXY(double xPos, double yPos, double& newPhi0, double& newDca,
+			double& newZ0) const;
+	GblHelixPrediction getPrediction(const Eigen::Vector3d& refPos,
+			const Eigen::Vector3d& uDir, const Eigen::Vector3d& vDir) const;
+
+private:
+	const double rinv; ///< curvature (1/Radius)
+	const double phi0; ///< azimuth at PCA (point of closest approach to origin in XY plane, defines arc-length S=0)
+	const double dca; ///< distance to origin in XY plane at PCA
+	const double dzds; ///< slope in ZS plane (dZ/dS)
+	const double z0; ///< offset in ZS plane
+	const double cosPhi0; ///< cos(phi0)
+	const double sinPhi0; ///< sin(phi0)
+	const double xRelCenter; ///< X position of circle center / R
+	const double yRelCenter; ///< Y position of circle center / R
+};
+
+/// Silicon tracker layer
+/**
+ * Silicon tracker layer with 1D or 2D measurement at fixed X-position.
+ * The measurement directions in the YZ plane can be orthogonal or non-orthogonal
+ * (but must be different).
+ */
+class GblSiliconLayer {
+public:
+	GblSiliconLayer(const std::string aName, double xPos, double yPos,
+			double zPos, double thickness, double uAngle, double uRes); // 1D measurement
+	GblSiliconLayer(const std::string aName, double xPos, double yPos,
+			double zPos, double thickness, double uAngle, double uRes,
+			double vAngle, double vRes); // 2D measurement
+	virtual ~GblSiliconLayer();
+	void print() const;
+	double getRadiationLength() const;
+	Eigen::Vector2d getResolution() const;
+	Eigen::Vector2d getPrecision() const;
+	Eigen::Vector3d getCenter() const;
+	GblHelixPrediction intersectWithHelix(GblSimpleHelix hlx) const;
+
+private:
+	std::string name; ///< name
+	int measDim; ///< measurement dimension (1 or 2)
+	double xbyx0; ///< normalized material thickness
+	Eigen::Vector3d center; ///< center
+	Eigen::Vector2d resolution; ///< measurements resolution
+	Eigen::Vector2d precision; ///< measurements precision
+	Eigen::Vector3d udir; ///< 1. measurement direction
+	Eigen::Vector3d vdir; ///< 2. measurement direction
+};
+
+}
+
+#endif /* EXAMPLESIT_H_ */

cpp/include/GblTrajectory.h

@@ -13,7 +13,7 @@
  *
  *
  *  \copyright
- *  Copyright (c) 2011 - 2017 Deutsches Elektronen-Synchroton,
+ *  Copyright (c) 2011 - 2018 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
@@ -186,6 +186,7 @@ private:
 	Eigen::MatrixXd externalDerivatives; // Derivatives for external measurements of composed trajectory
 	Eigen::VectorXd externalMeasurements; // Residuals for external measurements of composed trajectory
 	Eigen::VectorXd externalPrecisions; // Precisions for external measurements of composed trajectory
+	// linear equation system
 	VVector theVector; ///< Vector of linear equation system
 	BorderedBandMatrix theMatrix; ///< (Bordered band) matrix of linear equation system
 

cpp/include/MilleBinary.h

@@ -11,7 +11,7 @@
  *  \author Claus Kleinwort, DESY, 2011 (Claus.Kleinwort@desy.de)
  *
  *  \copyright
- *  Copyright (c) 2011 - 2017 Deutsches Elektronen-Synchroton,
+ *  Copyright (c) 2011 - 2018 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
@@ -50,10 +50,11 @@ namespace gbl {
  *     1   RMEAS, measured value   0                            -+
  *     2   local derivative        index of local derivative     |
  *     3   local derivative        index of local derivative     |
- *     4    ...                                                  | block
+ *     4   ...                                                   | block
  *         SIGMA, error (>0)       0                             |
  *         global derivative       label of global derivative    |
- *         global derivative       label of global derivative   -+
+ *         global derivative       label of global derivative    |
+ *         ...                                                  -+
  *         RMEAS, measured value   0
  *         local derivative        index of local derivative
  *         local derivative        index of local derivative
@@ -64,6 +65,18 @@ namespace gbl {
  *         ...
  *         global derivative       label of global derivative
  *\endverbatim
+ *
+ *  Special data block (other/debug information).
+ *  Contains no local derivatives and (error) SIGMA is negative (-Number of SPecial data words).
+ *
+ *\verbatim
+ *         real array              integer array
+ *         0.0                     0                   -+
+ *         -float(NSP)             0                    |
+ *         special data            special data         | special block (2+NSP words)
+ *         special data            special data         |
+ *         ...                                         -+
+ *\endverbatim
  */
 class MilleBinary {
 public:

cpp/src/GblTrajectory.cpp

@@ -11,7 +11,7 @@
  *  \author Claus Kleinwort, DESY, 2011 (Claus.Kleinwort@desy.de)
  *
  *  \copyright
- *  Copyright (c) 2011 - 2017 Deutsches Elektronen-Synchroton,
+ *  Copyright (c) 2011 - 2018 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
@@ -130,6 +130,10 @@
  *  based on std::vector<> for variable sized matrices.
  *  Several GBL methods are implemented as templates to allow more EIGEN compile time optimization.
  *
+ *  \section example_sec Examples
+ *  Technical examples are given in example1.cpp, example2.cpp and example3.cpp.
+ *  An example silicon tracker is described in exampleSit.cpp.
+ *
  *  \section ref_sec References
  *    - V. Blobel, C. Kleinwort, F. Meier,
  *      Fast alignment of a complex tracking detector using advanced track models,
@@ -591,8 +595,7 @@ void GblTrajectory::getFitToLocalJacobian(std::array<unsigned int, 5>& anIndex,
 
 	int nOffset = aPoint.getOffset();
 
-	for (unsigned int i = 0; i < 5; ++i)
-		anIndex[i] = 0;
+	anIndex = {}; // reset to 0
 	aJacobian.setZero();
 	if (nOffset < 0) // need interpolation
 			{
@@ -687,8 +690,7 @@ void GblTrajectory::getFitToKinkJacobian(std::array<unsigned int, 7>& anIndex,
 
 	int nOffset = aPoint.getOffset();
 
-	for (unsigned int i = 0; i < 7; ++i)
-		anIndex[i] = 0;
+	anIndex = {}; // reset to 0
 	aJacobian.setZero();
 
 	Matrix2d prevW, prevWJ, nextW, nextWJ;
@@ -978,6 +980,7 @@ void GblTrajectory::getResAndErr(unsigned int aData, bool used,
 	}
 	MatrixXd aMat = theMatrix.getBlockMatrix(numLocal, indLocal); // compressed (covariance) matrix
 	double aFitVar = aVec.transpose() * aMat * aVec; // variance from track fit
+	aFitVar *= aDownWeight; // account for down-weighting (of measurement in fit)
 	aMeasError = sqrt(aMeasVar); // error of measurement
 	if (used)
 		aResError = (aFitVar < aMeasVar ? sqrt(aMeasVar - aFitVar) : 0.); // error of (biased) residual

python/gblfit.py

@@ -12,7 +12,7 @@ Created on Jul 27, 2011
 # \author Claus Kleinwort, DESY, 2011 (Claus.Kleinwort@desy.de)
 #
 #  \copyright
-#  Copyright (c) 2011 - 2016 Deutsches Elektronen-Synchroton,
+#  Copyright (c) 2011 - 2018 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,6 +32,7 @@ import math
 from gblnum import BorderedBandMatrix
 from mille import MilleRecord
 
+
 ## User supplied point on (initial) trajectory.
 #
 #  Must have jacobians for propagation to previous or next point with offsets (first, last
@@ -62,6 +63,8 @@ class GblPoint(object):
     self.__measurement = None
     ## dimension of measurement (2D, 4D or 5D); int
     self.__measDim = 0
+    ## minimal precision to accept measurement                                                                                <
+    self.__measMinPrec = 0.
     ## transformation (to eigen-vectors of precision matrix); matrix(float)
     self.__measTransformation = None    
     ## scatterer at point: (transformation or None,) initial kinks, precision (inverse covariance matrix); list(matrix(float))
@@ -85,10 +88,12 @@ class GblPoint(object):
   # 
   #  @param aMeasurement measurement (projection (or None), residuals, precision
   #                       (diagonal of or full matrix)); list(matrix(float))
+  #  @param minPrecision Minimal precision to accept measurement
   #
-  def addMeasurement(self, aMeasurement):
+  def addMeasurement(self, aMeasurement, minPrecision=0.):
     self.__measurement = aMeasurement
     self.__measDim = aMeasurement[1].shape[0]
+    self.__measMinPrec = minPrecision
     if (aMeasurement[2].ndim == 2):  # full precision matrix, need to diagonalize
       eigenVal, eigenVec = np.linalg.eigh(aMeasurement[2])
       self.__measTransformation = eigenVec.T
@@ -121,6 +126,13 @@ class GblPoint(object):
   #
   def getMeasDim(self):
     return self.__measDim
+  
+  ## Retrieve minimal precision to accept measurement.
+  #
+  #  @return minimal precision to accept measurement; float
+  #
+  def getMeasMinPrec(self):
+    return self.__measMinPrec
 
   ## Add a (thin) scatterer to a point.
   #  
@@ -293,6 +305,7 @@ class GblPoint(object):
 
 #------------------------------------------------------------------------------ 
 
+
 ## Data (block) containing value, precision and derivatives for measurements, kinks and external seed.
 #  
 #  Created from attributes of GblPoints, used to construct linear equation system for track fit.
@@ -578,6 +591,9 @@ class GblData(object):
 #  points (corrections, covariance matrix) are not available and omission of
 #  measurements from a point is not possible.
 #  
+#  \section example_sec Examples
+# Technical examples are given in gbltst.py, an example silicon tracker in gblsit.py.
+#
 #  \section ref_sec References:  
 #    - V. Blobel, C. Kleinwort, F. Meier,
 #      Fast alignment of a complex tracking detector using advanced track models,
@@ -586,6 +602,7 @@ class GblData(object):
 #      NIM A, 673 (2012), 107-110, doi:10.1016/j.nima.2012.01.024
 # 
 
+
 ## General Broken Lines Trajectory.   
 #      
 class GblTrajectory(object):
@@ -664,7 +681,10 @@ class GblTrajectory(object):
     for d in self.__data:
       d.printData() 
 
-  ## Get data of trajectory.          
+  ## Get data of trajectory. 
+  #
+  # @return data blocks; list
+  #         
   def getData(self):
     return self.__data
  
@@ -887,6 +907,7 @@ class GblTrajectory(object):
     aVec = np.array(derLocal)  # compressed vector
     aMat = self.__matrix.getBlockMatrix(indLocal)  # compressed matrix     
     aFitVar = np.dot(aVec, np.dot(aMat, aVec.T))  # variance from track fit
+    aFitVar *= aDownWeight  # account for down-weighting (of measurement in fit)
     aMeasError = math.sqrt(aMeasVar)  # error of measurement
     if used:
       aResError = math.sqrt(aMeasVar - aFitVar) if aFitVar < aMeasVar else 0.  # error of biased residual
@@ -988,7 +1009,7 @@ class GblTrajectory(object):
       self.__points[-1].setOffset(nOffsets)
       self.__numOffsets = nOffsets + 1
       self.__numParameters = self.__numOffsets * len(self.__dimensions) \
-                           + self.__numCurvature + self.__numLocals       
+                           +self.__numCurvature + self.__numLocals       
 
     ## Calculate Jacobians to previous/next scatterer from point to point ones.
     def calcJacobians():
@@ -1024,6 +1045,7 @@ class GblTrajectory(object):
         if (aPoint.hasMeasurement()):
           nLabel = aPoint.getLabel()
           measDim = aPoint.getMeasDim()
+          measPrecision = aPoint.getMeasMinPrec()
           localDer = aPoint.getLocalDerivatives()
           globalLab = aPoint.getGlobalLabels()
           globalDer = aPoint.getGlobalDerivatives()
@@ -1032,7 +1054,7 @@ class GblTrajectory(object):
           labDer, matDer = self.__getFitToLocalJacobian(aPoint, measDim, nJacobian)
           matPDer = matDer if matP is None else np.dot(matP, matDer)
           for i in range(measDim):
-            if (aPrec[i] > 0.):
+            if (aPrec[i] > measPrecision):
               aData = GblData(nLabel, 1, aMeas[i], aPrec[i])
               aData.addDerivatives(i, labDer, matPDer, localDer, \
                                    globalLab, globalDer)

python/gblsit.py

+'''
+Created on 28 Sep 2018
+
+@author: kleinwrt
+'''
+
+## \file
+# silicon tracker example
+#
+# \author Claus Kleinwort, DESY, 2018 (Claus.Kleinwort@desy.de)
+#
+#  \copyright
+#  Copyright (c) 2018 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.
+
+import numpy as np
+import math
+import random
+import time
+from gblfit import GblPoint, GblTrajectory
+
+
+## Simulate and reconstruct helical tracks in silicon pixel and (1D or 2D) strip detectors. 
+#
+#  Create points on initial trajectory, create trajectory from points,
+#  fit and write trajectory to MP-II binary file (for rigid body alignment).
+#
+#  Setup:
+#   - Beam (mainly) in X direction
+#   - Constant magnetic field in Z direction
+#   - Silicon sensors measuring in YZ plane, orthogonal (pixel) or non-orthogonal (stereo strips) measurement systems
+#   - Multiple scattering in sensors (air inbetween ignored)
+#   - Curvilinear system (T,U,V) as local coordinate system and (q/p, slopes, offsets) as local track parameters
+#
+# \remark To exercise (mis)alignment different sets of layers (with different geometry) for simulation and reconstruction can be used.
+def exampleSit():
+
+  #
+  random.seed(47117)
+  
+  # magnetic field
+  bfac = 0.003  # B*c for 1 T
+  #bfac = 0.  # B*c for 0 T
+  # detector layers: name, position (x,y,z), thickness (X/X_0), (1 or 2) measurements (direction in YZ, resolution)  
+  det = gblSiliconDet([ ['PIX1', (2.0, 0., 0.), 0.0033, [(0., 0.0010), (90., 0.0020)] ],  # pixel
+                        ['PIX2', (3.0, 0., 0.), 0.0033, [(0., 0.0010), (90., 0.0020)] ],  # pixel
+                        ['PIX3', (4.0, 0., 0.), 0.0033, [(0., 0.0010), (90., 0.0020)] ],  # pixel
+                        ['S2D4', (6.0, 0., 0.), 0.0033, [(0., 0.0025), (5., 0.0025)] ],  # strip 2D, 5 deg stereo angle
+                        ['S2D5', (8.0, 0., 0.), 0.0033, [(0., 0.0025), (-5., 0.0025)] ],  # strip 2D, -5 deg stereo angle
+                        ['S2D6', (10.0, 0., 0.), 0.0033, [(0., 0.0025), (5., 0.0025)] ],  # strip 2D, 5 deg stereo angle
+                        ['S2D7', (12.0, 0., 0.), 0.0033, [(0., 0.0025), (-5., 0.0025)] ],  # strip 2D, -5 deg stereo angle
+                        ['S1D8', (15.0, 0., 0.), 0.0033, [(0., 0.0040)] ]  # strip 1D
+                      ], bfac)
+               
+  nTry = 1000  #: number of tries
+  qbyp = 0.2  # 5 GeV
+  binaryFile = open("milleBinary.dat", "wb")
+  #binaryFile = None
+  #
+  print " Gblsit $Rev$ ", nTry
+  #
+  start = time.clock()
+  Chi2Sum = 0.
+  NdfSum = 0
+  LostSum = 0.
+  #
+  for iTry in range(nTry):
+    # helix parameter for track generation
+    dca = random.gauss(0., 0.1)
+    z0 = random.gauss(0., 0.1)
+    phi0 = 0.5 * random.uniform(-1., 1.)
+    dzds = 0.3 * random.uniform(-1., 1.)
+    curv = bfac * qbyp * math.sqrt(1. + dzds * dzds)
+    genPar = [curv, phi0, dca, dzds, z0]
+    #print " genPar ", iTry, genPar
+    # generate hits
+    genHits = det.generateHits(qbyp, genPar)
+    # seed (with true parameters)
+    seedPar = genPar
+    seed = gblSimpleHelix(seedPar)
+    sOld = 0.
+    cosLambda = 1. / math.sqrt(1. + seedPar[3] ** 2)
+    sinLambda = seedPar[3] * cosLambda
+    # construct GBL trajectory
+    traj = GblTrajectory(bfac != 0.)
+    # add GBL points
+    for l, layer in enumerate(det.getLayers()):
+      # prediction from seeding helix
+      pred = layer.intersectWithHelix(seed)
+      measPred = pred.getMeasPred() 
+      sArc = pred.getArcLength()      
+      # residuals
+      res = np.array([genHits[l][0] - measPred[0], genHits[l][1] - measPred[1]])
+      # measurement precision
+      measPrec = np.array(layer.getPrecision())
+      # Curvilinear system: track direction T, U = Z x T / |Z x T|, V = T x U
+      # as local system
+      phi = seedPar[1] + seedPar[0] * sArc
+      cosPhi = math.cos(phi); sinPhi = math.sin(phi)
+      tuvDir = np.array([[cosLambda * cosPhi, cosLambda * sinPhi, sinLambda], \
+                     [-sinPhi, cosPhi, 0.], \
+                     [-sinLambda * cosPhi, -sinLambda * sinPhi, cosLambda]])
+      # projection matrix (local to measurement)
+      proL2m = np.linalg.inv(np.dot(tuvDir[1:3, :], np.linalg.inv(pred.getMeasSystemDirs())[:, :2]))
+      # propagation
+      jacPointToPoint = gblSimpleJacobian(sArc - sOld, cosLambda, bfac)
+      sOld = sArc
+      # point with (independent) measurements (in measurement system)
+      point = GblPoint(jacPointToPoint)
+      point.addMeasurement([proL2m, res, measPrec])
+      # global parameters for rigid body alignment?
+      if binaryFile is not None:
+        labels = [l * 10 + 1, l * 10 + 2, l * 10 + 3, l * 10 + 4, l * 10 + 5, l * 10 + 6]
+        labGlobal = np.array([labels, labels])
+        derGlobal = pred.getRigidBodyDerGlobal(layer.getCenter())[:2]
+        point.addGlobals(labGlobal, derGlobal)
+      # add scatterer to point
+      radlen = layer.getRadiationLength() / abs(pred.getCosIncidence())
+      scatErr = gblMultipleScatteringError(qbyp, radlen)  # simple model
+      if scatErr > 0.:
+        scat = np.array([0., 0.])
+        scatP = np.array([1. / scatErr ** 2, 1. / scatErr ** 2])
+        point.addScatterer([scat, scatP])
+      # add point to trajectory      
+      traj.addPoint(point)
+
+    # fit trajectory
+    Chi2, Ndf, Lost = traj.fit()
+    print " Record, Chi2, Ndf, Lost", iTry, Chi2, Ndf, Lost
+    # sum up    
+    Chi2Sum += Chi2
+    NdfSum += Ndf
+    LostSum += Lost
+    # write to binary file
+    if binaryFile is not None:
+      traj.milleOut(binaryFile)
+    
+  end = time.clock()
+  print " Time [s] ", end - start
+  print " Chi2Sum/NdfSum ", Chi2Sum / NdfSum
+  print " LostSum/nTry ", LostSum / nTry
+
+
+## Simple jacobian.
+#  
+#  Simple jacobian for (q/p, slopes, offsets) in curvilinear system
+#  quadratic in arc length difference.
+#
+#  @param ds arc length difference; float
+#  @param cosl cos(lambda); float
+#  @param bfac Bz*c; float
+#  @return jacobian to move by 'ds' on trajectory matrix(float)
+#
+def gblSimpleJacobian(ds, cosl, bfac):
+  jac = np.eye(5)
+  jac[1, 0] = -bfac * ds * cosl
+  jac[3, 0] = -0.5 * bfac * ds * ds * cosl
+  jac[3, 1] = ds
+  jac[4, 2] = ds  
+  return jac     
+
+
+## Multiple scattering error
+#
+# Simple model (Rossi, Greisen)
+#
+# @param[in] qbyp    q/p [1/GeV]; float
+# @param[in] xbyx0   thickness / radiation length; float
+#
+def gblMultipleScatteringError(qbyp, xbyx0):
+  return 0.015 * abs(qbyp) * math.sqrt(xbyx0)
+  
+
+## Silicon layer
+class gblSiliconLayer(object):
+
+  ## Constructor
+  #
+  # @param[in]  layer   layer description; list