GeneralBrokenLines V03-01-03
using EIGEN
example1.cpp
Go to the documentation of this file.
1/*
2 * example1.cpp
3 *
4 * Created on: Aug 24, 2011
5 * Author: kleinwrt
6 */
7
30#include <time.h>
31#include "GblUtilities.h"
32#include "GblTrajectory.h"
33
34using namespace gbl;
35using namespace Eigen;
36
37void example1() {
39
59//MP MilleBinary mille; // for producing MillePede-II binary file
60 unsigned int nTry = 1000; //: number of tries
61 unsigned int nLayer = 10; //: number of detector layers
62 std::cout << " Gbltst-eigen $Id$ " << nTry << ", " << nLayer << std::endl;
63
64 srand(4711);
65
66 clock_t startTime = clock();
67// track direction
68 double sinLambda = 0.3;
69 double cosLambda = sqrt(1.0 - sinLambda * sinLambda);
70 double sinPhi = 0.;
71 double cosPhi = sqrt(1.0 - sinPhi * sinPhi);
72// tDir = (cosLambda * cosPhi, cosLambda * sinPhi, sinLambda)
73// U = Z x T / |Z x T|, V = T x U
74 Matrix<double, 2, 3> uvDir;
75 uvDir(0, 0) = -sinPhi;
76 uvDir(0, 1) = cosPhi;
77 uvDir(0, 2) = 0.;
78 uvDir(1, 0) = -sinLambda * cosPhi;
79 uvDir(1, 1) = -sinLambda * sinPhi;
80 uvDir(1, 2) = cosLambda;
81// measurement resolution
82 Vector2d measErr;
83 measErr << 0.001, 0.001;
84 Vector2d measPrec; // (independent) precisions
85 measPrec << 1.0 / (measErr(0) * measErr(0)), 1.0 / (measErr(1) * measErr(1));
86 Matrix2d measInvCov; // inverse covariance matrix
87 measInvCov.setZero();
88 measInvCov(0, 0) = measPrec[0];
89 measInvCov(1, 1) = measPrec[1];
90// scattering error
91 Vector2d scatErr;
92 scatErr << 0.001, 0.001;
93 Vector2d scatPrec;
94 scatPrec << 1.0 / (scatErr(0) * scatErr(0)), 1.0 / (scatErr(1) * scatErr(1));
95// (RMS of) CurviLinear track parameters (Q/P, slopes, offsets)
96 Vector5d clPar;
97 Vector5d clErr;
98 clErr << 0.001, -0.1, 0.2, -0.15, 0.25;
99 Matrix5d clCov, clSeed;
100 unsigned int seedLabel = 99999;
101// additional parameters
102 Vector2d addPar;
103 addPar << 0.0025, -0.005;
104 std::vector<int> globalLabels;
105 globalLabels.push_back(4711);
106 globalLabels.push_back(4712);
107// global labels for MP
108 /*MP std::vector<int> globalLabels(2);
109 globalLabels[0] = 11;
110 globalLabels[1] = 12; */
111
112 double bfac = 0.2998; // Bz*c for Bz=1
113 double step = 1.5 / cosLambda; // constant steps in RPhi
114
115 double Chi2Sum = 0.;
116 int NdfSum = 0;
117 double LostSum = 0.;
118 int numFit = 0;
119
120 for (unsigned int iTry = 1; iTry <= nTry; ++iTry) {
121 // curvilinear track parameters
122 for (unsigned int i = 0; i < 5; ++i) {
123 clPar[i] = clErr[i] * unrm();
124 }
125 clCov.setZero();
126 for (unsigned int i = 0; i < 5; ++i) {
127 clCov(i, i) = 1.0 * (clErr[i] * clErr[i]);
128 }
129// std::cout << " Try " << iTry << ":" << clPar << std::endl;
130 Matrix2d addDer;
131 addDer.setZero();
132 addDer(0, 0) = 1.;
133 addDer(1, 1) = 1.;
134// arclength
135 //double s = 0.;
136 Matrix5d jacPointToPoint;
137 jacPointToPoint.setIdentity();
138// create list of points
139 std::vector<GblPoint> listOfPoints;
140 listOfPoints.reserve(2 * nLayer);
141
142 for (unsigned int iLayer = 0; iLayer < nLayer; ++iLayer) {
143// std::cout << " Layer " << iLayer << ", " << s << std::endl;
144// measurement directions
145 double sinStereo = (iLayer % 2 == 0) ? 0. : 0.1;
146 double cosStereo = sqrt(1.0 - sinStereo * sinStereo);
147 Matrix<double, 3, 2> mDirT;
148 mDirT.setZero();
149 mDirT(1, 0) = cosStereo;
150 mDirT(2, 0) = sinStereo;
151 mDirT(1, 1) = -sinStereo;
152 mDirT(2, 1) = cosStereo;
153// projection measurement to local (curvilinear uv) directions (duv/dm)
154 Matrix2d proM2l = uvDir * mDirT;
155// projection local (uv) to measurement directions (dm/duv)
156 Matrix2d proL2m = proM2l.inverse();
157 // point with (independent) measurements (in measurement system)
158 GblPoint pointMeas(jacPointToPoint);
159 // measurement - prediction in measurement system with error
160 Vector2d meas = proL2m * clPar.tail(2);
161 //MP meas += addDer * addPar; // additional parameters
162 for (unsigned int i = 0; i < 2; ++i) {
163 meas[i] += measErr[i] * unrm();
164 }
165 pointMeas.addMeasurement(proL2m, meas, measPrec);
166 /* point with (correlated) measurements (in local system)
167 GblPoint pointMeas(jacPointToPoint);
168 // measurement - prediction in local system with error
169 Vector2d meas;
170 for (unsigned int i = 0; i < 2; ++i) {
171 meas[i] = measErr[i] * unrm(); // + addDer(i, 0) * 0.0075;
172 }
173 meas = proM2l * meas + clPar.tail<2>();
174 Matrix2d localInvCov = proL2m.transpose() * measInvCov * proL2m;
175 pointMeas.addMeasurement(meas, localInvCov); */
176
177 // additional local parameters?
178// point.addLocals(addDer);
179//MP point.addGlobals(globalLabels, addDer);
180 addDer *= -1.; // Der flips sign every measurement
181// add point to trajectory
182 listOfPoints.push_back(pointMeas);
183 unsigned int iLabel = listOfPoints.size();
184 if (iLabel == seedLabel) {
185 clSeed = clCov.inverse();
186 }
187// propagate to scatterer
188 jacPointToPoint = gblSimpleJacobian(step, cosLambda, bfac);
189 //jac2 = gblSimpleJacobian2(step, cosLambda, bfac);
190 clPar = jacPointToPoint * clPar;
191 clCov = jacPointToPoint * clCov * jacPointToPoint.transpose();
192 //s += step;
193 if (iLayer < nLayer - 1) {
194 Vector2d scat(0., 0.);
195 // point with scatterer
196 GblPoint pointScat(jacPointToPoint);
197 pointScat.addScatterer(scat, scatPrec);
198 listOfPoints.push_back(pointScat);
199 iLabel = listOfPoints.size();
200 if (iLabel == seedLabel) {
201 clSeed = clCov.inverse();
202 }
203 // scatter a little
204 for (unsigned int i = 0; i < 2; ++i) {
205 clPar[i + 1] += scatErr[i] * unrm();
206 clCov(i + 1, i + 1) += scatErr[i] * scatErr[i];
207 }
208 // propagate to next measurement layer
209 clPar = jacPointToPoint * clPar;
210 clCov = jacPointToPoint * clCov * jacPointToPoint.transpose();
211 //s += step;
212 }
213 }
214//
215 // create trajectory
216 GblTrajectory traj(listOfPoints);
217 //GblTrajectory traj(listOfPoints, seedLabel, clSeed); // with external seed
218 //traj.printPoints();
219 /*
220 if (not traj.isValid()) {
221 std::cout << " Invalid GblTrajectory -> skip" << std::endl;
222 continue;
223 }*/
224// fit trajectory
225 double Chi2;
226 int Ndf;
227 double lostWeight;
228 traj.fit(Chi2, Ndf, lostWeight);
229 double bandCond = traj.getBandCondition();
230 std::cout << " Fit: " << Chi2 << ", " << Ndf << ", " << lostWeight
231 << ", " << bandCond << std::endl;
232 /* look at (track parameter) corrections
233 VectorXd aCorrection(5);
234 MatrixXd aCovariance(5, 5);
235 traj.getResults(1, aCorrection, aCovariance);
236 std::cout << " cor " << std::endl << aCorrection << std::endl;
237 std::cout << " cov " << std::endl << aCovariance << std::endl;
238 */
239 /* look at residuals
240 for (unsigned int label = 1; label <= listOfPoints.size(); ++label) {
241 unsigned int numData = 0;
242 std::cout << " measResults, label " << label << std::endl;
243 VectorXd residuals(2), measErr(2), resErr(2), downWeights(2);
244 traj.getMeasResults(label, numData, residuals, measErr, resErr,
245 downWeights);
246 std::cout << " measResults, numData " << numData << std::endl;
247 for (unsigned int i = 0; i < numData; ++i) {
248 std::cout << " measResults " << label << " " << i << " "
249 << residuals[i] << " " << measErr[i] << " " << resErr[i]
250 << std::endl;
251 }
252 } */
253// debug printout
254 //traj.printTrajectory(1);
255 //traj.printPoints(1);
256 //traj.printData();
257// write to MP binary file
258//MP traj.milleOut(mille);
259 Chi2Sum += Chi2;
260 NdfSum += Ndf;
261 LostSum += lostWeight;
262 numFit++;
263 }
264
265 clock_t endTime = clock();
266 double diff = endTime - startTime;
267 double cps = CLOCKS_PER_SEC;
268 std::cout << " Time elapsed " << diff / cps << " s" << std::endl;
269 std::cout << " Chi2/Ndf = " << Chi2Sum / NdfSum << std::endl;
270 std::cout << " Tracks fitted " << numFit << std::endl;
271 if (LostSum > 0.)
272 std::cout << " Weight lost " << LostSum << std::endl;
273}
274
GblTrajectory definition.
Definitions for GBL utilities.
Point on trajectory.
Definition: GblPoint.h:62
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.
Definition: GblPoint.h:238
void addScatterer(const Eigen::Vector2d &aResiduals, const Eigen::Matrix2d &aPrecision)
Add a thin scatterer to a point.
Definition: GblPoint.cpp:159
GBL trajectory.
Definition: GblTrajectory.h:50
double getBandCondition() const
Get condition from band (decomposition)
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight, const std::string &optionList="", unsigned int aLabel=0)
Perform fit of (valid) trajectory.
void example1()
Definition: example1.cpp:37
Namespace for the general broken lines package.
double unrm()
unit normal distribution, Box-Muller method, polar form
Eigen::Matrix< double, 5, 1 > Vector5d
Matrix5d gblSimpleJacobian(double ds, double cosl, double bfac)
Simple jacobian.
Eigen::Matrix< double, 5, 5 > Matrix5d
Definition: GblData.h:46