example1.cpp 9.72 KB
Newer Older
1
2
3
4
5
6
7
/*
 * example1.cpp
 *
 *  Created on: Aug 24, 2011
 *      Author: kleinwrt
 */

8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
/** \file
 *  Example application.
 *
 *  \author Claus Kleinwort, DESY, 2011 (Claus.Kleinwort@desy.de)
 *
 *  \copyright
 *  Copyright (c) 2011 - 2016 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.
 */

30
31
32
33
34
#include <time.h>
#include "example1.h"
#include "TRandom3.h"
#include "GblTrajectory.h"

35
36
using namespace gbl;

37
TMatrixD gblSimpleJacobian(double ds, double cosl, double bfac) {
38
39
40
41
42
43
44
	/// 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
	 */
45
46
47
48
49
50
51
52
53
54
	TMatrixD jac(5, 5);
	jac.UnitMatrix();
	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;
}

void example1() {
55
56
	/// Simple example.
	/**
57
58
59
60
61
62
	 * 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.
	 *
	 * Equidistant measurement layers and thin scatterers, propagation
	 * with simple jacobian (quadratic in arc length differences).
63
	 * Curvilinear system (U,V,T) as local coordinate system.
64
65
66
	 *
	 * This example simulates and refits tracks in a system of planar detectors
	 * with 2D measurements in a constant magnet field in Z direction using
67
68
	 * the curvilinear system as local system and (Q/P, slopes, offsets) as
	 * local track parameters. The true track parameters are
69
70
71
72
73
	 * 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.
74
	 */
75
76

//MP	MilleBinary mille; // for producing MillePede-II binary file
77
78
79
80
81
82
83
84
85
86
87
88
89
90
	unsigned int nTry = 10000; //: number of tries
	unsigned int nLayer = 10; //: number of detector layers
	std::cout << " Gbltst $Rev$ " << nTry << ", " << nLayer << std::endl;

	TRandom *r = new TRandom3();

	clock_t startTime = clock();
// track direction
	double sinLambda = 0.3;
	double cosLambda = sqrt(1.0 - sinLambda * sinLambda);
	double sinPhi = 0.;
	double cosPhi = sqrt(1.0 - sinPhi * sinPhi);
// tDir = (cosLambda * cosPhi, cosLambda * sinPhi, sinLambda)
// U = Z x T / |Z x T|, V = T x U
91
92
93
94
95
96
97
	TMatrixD uvDir(2, 3);
	uvDir[0][0] = -sinPhi;
	uvDir[0][1] = cosPhi;
	uvDir[0][2] = 0.;
	uvDir[1][0] = -sinLambda * cosPhi;
	uvDir[1][1] = -sinLambda * sinPhi;
	uvDir[1][2] = cosLambda;
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
// measurement resolution
	TVectorD measErr(2);
	measErr[0] = 0.001;
	measErr[1] = 0.001;
	TVectorD measPrec(2); // (independent) precisions
	measPrec[0] = 1.0 / (measErr[0] * measErr[0]);
	measPrec[1] = 1.0 / (measErr[1] * measErr[1]);
	TMatrixDSym measInvCov(2); // inverse covariance matrix
	measInvCov.Zero();
	measInvCov[0][0] = measPrec[0];
	measInvCov[1][1] = measPrec[1];
// scattering error
	TVectorD scatErr(2);
	scatErr[0] = 0.001;
	scatErr[1] = 0.001;
	TVectorD scatPrec(2);
	scatPrec[0] = 1.0 / (scatErr[0] * scatErr[0]);
	scatPrec[1] = 1.0 / (scatErr[1] * scatErr[1]);
116
// (RMS of) CurviLinear track parameters (Q/P, slopes, offsets)
117
118
119
120
121
122
123
124
125
	TVectorD clPar(5);
	TVectorD clErr(5);
	clErr[0] = 0.001;
	clErr[1] = -0.1;
	clErr[2] = 0.2;
	clErr[3] = -0.15;
	clErr[4] = 0.25;
	TMatrixDSym clCov(5), clSeed(5);
	unsigned int seedLabel = 0;
126
127
128
129
// additional parameters
	TVectorD addPar(2);
	addPar[0] = 0.0025;
	addPar[1] = -0.005;
130
131
132
	std::vector<int> globalLabels;
	globalLabels.push_back(4711);
	globalLabels.push_back(4712);
133
134
135
136
// global labels for MP
	/*MP	std::vector<int> globalLabels(2);
	 globalLabels[0] = 11;
	 globalLabels[1] = 12; */
137
138
139
140
141
142
143

	double bfac = 0.2998; // Bz*c for Bz=1
	double step = 1.5 / cosLambda; // constant steps in RPhi

	double Chi2Sum = 0.;
	int NdfSum = 0;
	double LostSum = 0.;
144
	int numFit = 0;
145

146
	for (unsigned int iTry = 1; iTry <= nTry; ++iTry) {
147
		// curvilinear track parameters
148
		for (unsigned int i = 0; i < 5; ++i) {
149
150
151
			clPar[i] = clErr[i] * r->Gaus();
		}
		clCov.Zero();
152
		for (unsigned int i = 0; i < 5; ++i) {
153
154
155
			clCov[i][i] = 1.0 * (clErr[i] * clErr[i]);
		}
//		std::cout << " Try " << iTry << ":" << clPar << std::endl;
156
		TMatrixD addDer(2, 2);
157
158
		addDer.Zero();
		addDer[0][0] = 1.;
159
160
		addDer[1][1] = 1.;
// arclength
161
162
163
		double s = 0.;
		TMatrixD jacPointToPoint(5, 5);
		jacPointToPoint.UnitMatrix();
164
165
// create list of points
		std::vector<GblPoint> listOfPoints;
166
		listOfPoints.reserve(2 * nLayer);
167

168
		for (unsigned int iLayer = 0; iLayer < nLayer; ++iLayer) {
169
170
//			std::cout << " Layer " << iLayer << ", " << s << std::endl;
//     measurement directions
171
			double sinStereo = (iLayer % 2 == 0) ? 0. : 0.1;
172
			double cosStereo = sqrt(1.0 - sinStereo * sinStereo);
173
174
175
			TMatrixD mDirT(3, 2);
			mDirT.Zero();
			mDirT[1][0] = cosStereo;
176
177
178
			mDirT[2][0] = sinStereo;
			mDirT[1][1] = -sinStereo;
			mDirT[2][1] = cosStereo;
179
180
// projection measurement to local (curvilinear uv) directions (duv/dm)
			TMatrixD proM2l = uvDir * mDirT;
181
// projection local (uv) to measurement directions (dm/duv)
182
183
			TMatrixD proL2m = proM2l;
			proL2m.Invert();
184
185
186
187
// point with (independent) measurements (in measurement system)
			GblPoint point(jacPointToPoint);
// measurement - prediction in measurement system with error
			TVectorD meas = proL2m * clPar.GetSub(3, 4);
188
//MP			meas += addDer * addPar; // additional parameters
189
			for (unsigned int i = 0; i < 2; ++i) {
190
191
192
193
194
195
196
				meas[i] += measErr[i] * r->Gaus();
			}
			point.addMeasurement(proL2m, meas, measPrec);
			/* point with (correlated) measurements (in local system)
			 GblPoint point(jacPointToPoint);
			 // measurement - prediction in local system with error
			 TVectorD meas(2);
197
			 for (unsigned int i = 0; i < 2; ++i) {
198
199
200
201
202
203
204
205
206
			 meas[i] = measErr[i] * r->Gaus() + addDer(i,0) * 0.0075;
			 }
			 meas = proM2l * meas + clPar.GetSub(3, 4);
			 TMatrixDSym localInvCov = measInvCov;
			 localInvCov.SimilarityT(proL2m);
			 point.addMeasurement(meas, localInvCov);

			 // additional local parameters? */
//			point.addLocals(addDer);
207
//MP			point.addGlobals(globalLabels, addDer);
208
209
			addDer *= -1.; // Der flips sign every measurement
// add point to trajectory
210
211
			listOfPoints.push_back(point);
			unsigned int iLabel = listOfPoints.size();
212
213
214
215
216
217
218
219
220
221
222
223
224
225
			if (iLabel == seedLabel) {
				clSeed = clCov.Invert();
			}
// propagate to scatterer
			jacPointToPoint = gblSimpleJacobian(step, cosLambda, bfac);
			clPar = jacPointToPoint * clPar;
			clCov = clCov.Similarity(jacPointToPoint);
			s += step;
			if (iLayer < nLayer - 1) {
				TVectorD scat(2);
				scat.Zero();
				// point with scatterer
				GblPoint point(jacPointToPoint);
				point.addScatterer(scat, scatPrec);
226
227
				listOfPoints.push_back(point);
				iLabel = listOfPoints.size();
228
229
230
231
				if (iLabel == seedLabel) {
					clSeed = clCov.Invert();
				}
				// scatter a little
232
				for (unsigned int i = 0; i < 2; ++i) {
233
234
235
236
237
238
239
240
241
242
					clPar[i + 1] += scatErr[i] * r->Gaus();
					clCov[i + 1] += scatErr[i] * scatErr[i];
				}
				// propagate to next measurement layer
				clPar = jacPointToPoint * clPar;
				clCov = clCov.Similarity(jacPointToPoint);
				s += step;
			}
		}
//
243
244
245
		// create trajectory
		//GblTrajectory traj(listOfPoints);
		GblTrajectory traj(listOfPoints, seedLabel, clSeed); // with external seed
246
		//traj.printPoints();
247
248
249
250
		if (not traj.isValid()) {
			std::cout << " Invalid GblTrajectory -> skip" << std::endl;
			continue;
		}
251
252
253
254
255
256
// fit trajectory
		double Chi2;
		int Ndf;
		double lostWeight;
		traj.fit(Chi2, Ndf, lostWeight);
//		std::cout << " Fit: " << Chi2 << ", " << Ndf << ", " << lostWeight << std::endl;
257
258
259
		/*		 TVectorD aCorrection(5);
		 TMatrixDSym aCovariance(5);
		 traj.getResults(1, aCorrection, aCovariance);
260
261
262
263
		 std::cout << " cor " << std::endl;
		 aCorrection.Print();
		 std::cout << " cov " << std::endl;
		 aCovariance.Print(); */
264
    // look at residuals
265
266
     for (unsigned int label=1; label<=listOfPoints.size(); ++label) {
       unsigned int numData=0;
267
       //std::cout << " measResults, label " << label << std::endl;
268
269
       TVectorD residuals(2), measErr(2), resErr(2), downWeights(2);
       traj.getMeasResults(label, numData, residuals, measErr, resErr, downWeights);
270
271
272
273
274
275
       //std::cout << " measResults, numData " << numData << std::endl;
       /* residuals.Print(); measErr.Print(); resErr.Print();
		for (unsigned int i = 0; i < numData; ++i) {
			std::cout << " measResults " << label << " " << i << " " << residuals[i] << " " << measErr[i] << " " << resErr[i] << std::endl;
		} */
     }
276
277
278
279
// debug printout
		//traj.printTrajectory();
		//traj.printPoints();
		//traj.printData();
280
// write to MP binary file
281
//MP		traj.milleOut(mille);
282
283
284
		Chi2Sum += Chi2;
		NdfSum += Ndf;
		LostSum += lostWeight;
285
		numFit++;
286
287
288
289
290
291
	}
	clock_t endTime = clock();
	double diff = endTime - startTime;
	double cps = CLOCKS_PER_SEC;
	std::cout << " Time elapsed " << diff / cps << " s" << std::endl;
	std::cout << " Chi2/Ndf = " << Chi2Sum / NdfSum << std::endl;
292
	std::cout << " Tracks fitted " << numFit << std::endl;
293
294
}