exampleSit.cpp 12.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
/*
 * exampleSit.cpp
 *
 *  Created on: 11 Oct 2018
 *      Author: kleinwrt
 */

/** \file
 *  Example silicon tracker application.
 *
 *  \author Claus Kleinwort, DESY, 2018 (Claus.Kleinwort@desy.de)
 *
 *  \copyright
Claus Kleinwort's avatar
Claus Kleinwort committed
14
 *  Copyright (c) 2018-2021 Deutsches Elektronen-Synchroton,
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
 *  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.
 */

#include <time.h>
#include "exampleSit.h"
#include "GblTrajectory.h"

using namespace gbl;
using namespace Eigen;

/// Silicon tracker example
/**
 * 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
Claus Kleinwort's avatar
Claus Kleinwort committed
48
 *   - Multiple scattering in sensors (air in between ignored)
49
50
51
52
 *   - 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.
Claus Kleinwort's avatar
Claus Kleinwort committed
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
 *
 * Example steering file for Millepede-II (B=0):
 * \code{.unparsed}
 * Cfiles
 * milleBinaryISN.dat
 *
 * method inversion 3 0.1
 * chiscut 30. 6.
 * printcounts
 * ! fix first pixel and last stereo layer as reference
 * parameter
 *   1  0.  -1.
 *   2  0.  -1.
 *   3  0.  -1.
 *   4  0.  -1.
 *   5  0.  -1.
 *   6  0.  -1.
 *  61  0.  -1.
 *  62  0.  -1.
 *  63  0.  -1.
 *  64  0.  -1.
 *  65  0.  -1.
 *  66  0.  -1.
 * end
 * \endcode
78
79
80
81
82
 */
void exampleSit() {

	// detector layers (ordered in X):
	// name, position (x,y,z), thickness (X/X_0), (1 or 2) measurements (direction in YZ, resolution)
Claus Kleinwort's avatar
Claus Kleinwort committed
83
	std::vector<GblDetectorLayer> layers;
84
	layers.push_back(
Claus Kleinwort's avatar
Claus Kleinwort committed
85
			CreateLayerSit("PIX1", 0, 2.0, 0., 0., 0.0033, 0., 0.0010, 90.,
86
87
					0.0020)); // pixel
	layers.push_back(
Claus Kleinwort's avatar
Claus Kleinwort committed
88
			CreateLayerSit("PIX2", 1, 3.0, 0., 0., 0.0033, 0., 0.0010, 90.,
89
90
					0.0020)); // pixel
	layers.push_back(
Claus Kleinwort's avatar
Claus Kleinwort committed
91
			CreateLayerSit("PIX3", 2, 4.0, 0., 0., 0.0033, 0., 0.0010, 90.,
92
93
					0.0020)); // pixel
	layers.push_back(
Claus Kleinwort's avatar
Claus Kleinwort committed
94
			CreateLayerSit("S2D4", 3, 6.0, 0., 0., 0.0033, 0., 0.0025, 5.0,
95
96
					0.0025)); // strip 2D, +5 deg stereo angle
	layers.push_back(
Claus Kleinwort's avatar
Claus Kleinwort committed
97
			CreateLayerSit("S2D5", 4, 8.0, 0., 0., 0.0033, 0., 0.0025, -5.,
98
99
					0.0025)); // strip 2D, -5 deg stereo angle
	layers.push_back(
Claus Kleinwort's avatar
Claus Kleinwort committed
100
			CreateLayerSit("S2D6", 5, 10., 0., 0., 0.0033, 0., 0.0025, 5.0,
101
102
					0.0025)); // strip 2D, +5 deg stereo angle
	layers.push_back(
Claus Kleinwort's avatar
Claus Kleinwort committed
103
			CreateLayerSit("S2D7", 6, 12., 0., 0., 0.0033, 0., 0.0025, -5.,
104
					0.0025)); // strip 2D, -5 deg stereo angle
Claus Kleinwort's avatar
Claus Kleinwort committed
105
106
	layers.push_back(
			CreateLayerSit("S1D8", 7, 15., 0., 0., 0.0033, 0., 0.0040)); // strip 1D, no sensitivity to Z
107
108
109
110
111
112

	/* print layers
	 for (unsigned int iLayer = 0; iLayer < layers.size(); ++iLayer) {
	 layers[iLayer].print();
	 } */

Claus Kleinwort's avatar
Claus Kleinwort committed
113
	unsigned int nTry = 10000; //: number of tries
Claus Kleinwort's avatar
Claus Kleinwort committed
114
	std::cout << " GblSit $Id$ " << nTry << ", " << layers.size()
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
			<< std::endl;
	srand(4711);
	clock_t startTime = clock();

	double qbyp = 0.2; // 5 GeV
	const double bfac = 0.003;  // B*c for 1 T
	// const double bfac = 0.;  // B*c for 0 T

	MilleBinary mille; // for producing MillePede-II binary file

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

	for (unsigned int iTry = 0; iTry < nTry; ++iTry) {

		// helix parameter for track generation
		const double genDca = 0.1 * unrm(); // normal
		const double genZ0 = 0.1 * unrm(); // normal
		const double genPhi0 = 0.2 * (2. * unif() - 1.); // uniform
		const double genDzds = 0.3 * (2. * unif() - 1.); // uniform
		const double genCurv = bfac * qbyp * sqrt(1. + genDzds * genDzds);

		//
		// generate hits
		//
		std::vector<Vector2d> hits;
		double curv(genCurv), phi0(genPhi0), dca(genDca), dzds(genDzds), z0(
				genZ0);
		const double cosLambda = 1. / sqrt(1. + dzds * dzds);
		for (unsigned int iLayer = 0; iLayer < layers.size(); ++iLayer) {
			// local constant (Bfield) helix
			GblSimpleHelix hlx = GblSimpleHelix(curv, phi0, dca, dzds, z0);
			// prediction from local helix
			GblHelixPrediction pred = layers[iLayer].intersectWithHelix(hlx);
			// std::cout << " layer " << iLayer << " arc-length " << pred.getArcLength() << std::endl;
			Vector2d meas = pred.getMeasPred();
			// smear according to resolution
			Vector2d sigma = layers[iLayer].getResolution();
			meas[0] += sigma[0] * unrm();
			meas[1] += sigma[1] * unrm();
			// save hit
			hits.push_back(meas);
			// scatter at intersection point
Claus Kleinwort's avatar
Claus Kleinwort committed
160
			Vector3d measPos = pred.getPosition();
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
			double radlen = layers[iLayer].getRadiationLength()
					/ fabs(pred.getCosIncidence());
			double errMs = gblMultipleScatteringError(qbyp, radlen); // simple model
			// move to intersection point
			hlx.moveToXY(measPos[0], measPos[1], phi0, dca, z0); // update phi0, dca, z0
			phi0 += unrm() * errMs / cosLambda; // scattering for phi
			dzds += unrm() * errMs / (cosLambda * cosLambda); // scattering for dzds
			GblSimpleHelix newhlx = GblSimpleHelix(curv, phi0, dca, dzds, z0); // after scattering
			// move back
			newhlx.moveToXY(-measPos[0], -measPos[1], phi0, dca, z0); // update phi0, dca, z0
		}

		//
		// fit track with GBL
		//
		// seed (with true parameters)
		double seedCurv(genCurv), seedPhi0(genPhi0), seedDca(genDca), seedDzds(
				genDzds), seedZ0(genZ0);
		GblSimpleHelix seed = GblSimpleHelix(seedCurv, seedPhi0, seedDca,
				seedDzds, seedZ0);
		// (previous) arc-length
		double sOld = 0.;
		const double cosLambdaSeed = 1. / sqrt(1. + (seedDzds * seedDzds));
		// list of points on trajectory
		std::vector<GblPoint> listOfPoints;
		for (unsigned int iLayer = 0; iLayer < layers.size(); ++iLayer) {
			// std::cout << " hit " << iLayer << " " << hits[iLayer].transpose() << std::endl;
Claus Kleinwort's avatar
Claus Kleinwort committed
188
			GblDetectorLayer& layer = layers[iLayer];
189
			// prediction from seeding helix
Claus Kleinwort's avatar
Claus Kleinwort committed
190
			GblHelixPrediction pred = layer.intersectWithHelix(seed);
191
192
			double sArc = pred.getArcLength(); // arc-length
			Vector2d measPrediction = pred.getMeasPred(); // measurement prediction
Claus Kleinwort's avatar
Claus Kleinwort committed
193
			Vector2d measPrecision = layer.getPrecision(); // measurement precision
194
195
196
			// residuals
			Vector2d res(hits[iLayer][0] - measPrediction[0],
					hits[iLayer][1] - measPrediction[1]);
Claus Kleinwort's avatar
Claus Kleinwort committed
197
198
			// transformation global system to local (curvilinear) (u,v) (matrix from row vectors)
			Matrix<double, 2, 3> transG2l = pred.getCurvilinearDirs();
199
			// transformation measurement system to global system
Claus Kleinwort's avatar
Claus Kleinwort committed
200
			Matrix3d transM2g = layer.getMeasSystemDirs().inverse();
201
202
203
204
205
			// projection matrix (measurement plane to local (u,v))
			Matrix2d proM2l = transG2l * transM2g.block<3, 2>(0, 0); // skip measurement normal
			// projection matrix (local (u,v) to measurement plane)
			Matrix2d proL2m = proM2l.inverse();
			// propagation
Claus Kleinwort's avatar
Claus Kleinwort committed
206
207
			Matrix5d jacPointToPoint = gblSimpleJacobian(
					(sArc - sOld) / cosLambdaSeed, cosLambdaSeed, bfac);
208
209
210
211
212
213
			sOld = sArc;
			// point with (independent) measurements (in measurement system)
			GblPoint point(jacPointToPoint);
			point.addMeasurement(proL2m, res, measPrecision);
			// global labels and parameters for rigid body alignment
			std::vector<int> labGlobal(6);
Claus Kleinwort's avatar
Claus Kleinwort committed
214
			unsigned int layerID = layer.getLayerID();
215
			for (int p = 0; p < 6; p++)
Claus Kleinwort's avatar
Claus Kleinwort committed
216
				labGlobal[p] = layerID * 10 + p + 1;
Claus Kleinwort's avatar
Claus Kleinwort committed
217
218
219
220
			Vector3d pos = pred.getPosition();
			Vector3d dir = pred.getDirection();
			Matrix<double, 2, 6> derGlobal = layer.getRigidBodyDerLocal(pos,
					dir);
221
222
			point.addGlobals(labGlobal, derGlobal);
			// add scatterer to point
Claus Kleinwort's avatar
Claus Kleinwort committed
223
			double radlen = layer.getRadiationLength()
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
					/ fabs(pred.getCosIncidence());
			double errMs = gblMultipleScatteringError(qbyp, radlen); // simple model
			if (errMs > 0.) {
				Vector2d scat(0., 0.);
				Vector2d scatPrec(1. / (errMs * errMs), 1. / (errMs * errMs)); // scattering precision matrix is diagonal in curvilinear system
				point.addScatterer(scat, scatPrec);
			}
			// add point to trajectory
			listOfPoints.push_back(point);
		}
		// create trajectory
		GblTrajectory traj(listOfPoints, bfac != 0.);
		// fit trajectory
		double Chi2;
		int Ndf;
		double lostWeight;
		unsigned int ierr = traj.fit(Chi2, Ndf, lostWeight);
		// std::cout << " Fit " << iTry << ": "<< Chi2 << ", " << Ndf << ", " << lostWeight << std::endl;
		// successfully fitted?
		if (!ierr) {
			// write to MP binary file
			traj.milleOut(mille);
			// update statistics
			Chi2Sum += Chi2;
			NdfSum += Ndf;
			LostSum += lostWeight;
			numFit++;
		}
	}
	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;
	std::cout << " Tracks fitted " << numFit << std::endl;
}

namespace gbl {

/// Create a silicon layer with 1D measurement.
/**
Claus Kleinwort's avatar
Claus Kleinwort committed
265
266
 * Create silicon layer with 1D measurement (u) at fixed X-position.
 *
267
 * \param [in] aName      name
Claus Kleinwort's avatar
Claus Kleinwort committed
268
 * \param [in] layer      layer ID
269
270
271
272
273
274
275
 * \param [in] xPos       X-position (of center)
 * \param [in] yPos       Y-position (of center)
 * \param [in] zPos       Z-position (of center)
 * \param [in] thickness  thickness / radiation_length
 * \param [in] uAngle     angle of u-direction in YZ plane
 * \param [in] uRes       resolution in u-direction
 */
Claus Kleinwort's avatar
Claus Kleinwort committed
276
277
GblDetectorLayer CreateLayerSit(const std::string aName, unsigned int layer,
		double xPos, double yPos, double zPos, double thickness, double uAngle,
Claus Kleinwort's avatar
Claus Kleinwort committed
278
279
280
281
282
283
284
285
286
287
		double uRes) {
	Vector3d aCenter(xPos, yPos, zPos);
	Vector2d aResolution(uRes, 0.);
	Vector2d aPrecision(1. / (uRes * uRes), 0.);
	Matrix3d measTrafo;
	const double cosU = cos(uAngle / 180. * M_PI);
	const double sinU = sin(uAngle / 180. * M_PI);
	measTrafo << 0., cosU, sinU, 0., -sinU, cosU, 1., 0., 0.; // U,V,N
	Matrix3d alignTrafo;
	alignTrafo << 0., 1., 0., 0., 0., 1., 1., 0., 0.; // Y,Z,X
Claus Kleinwort's avatar
Claus Kleinwort committed
288
	return GblDetectorLayer(aName, layer, 1, thickness, aCenter, aResolution,
Claus Kleinwort's avatar
Claus Kleinwort committed
289
			aPrecision, measTrafo, alignTrafo);
290
291
292
293
}

/// Create a silicon layer with 2D measurement.
/**
Claus Kleinwort's avatar
Claus Kleinwort committed
294
295
296
297
 * Create silicon layer with 2D measurement (u,v) at fixed X-position.
 * The measurement directions in the YZ plane can be orthogonal or non-orthogonal
 * (but must be different).
 *
298
 * \param [in] aName      name
Claus Kleinwort's avatar
Claus Kleinwort committed
299
 * \param [in] layer      layer ID
300
301
302
303
304
305
306
307
308
 * \param [in] xPos       X-position (of center)
 * \param [in] yPos       Y-position (of center)
 * \param [in] zPos       Z-position (of center)
 * \param [in] thickness  thickness / radiation_length
 * \param [in] uAngle     angle of u-direction in YZ plane
 * \param [in] uRes       resolution in u-direction
 * \param [in] vAngle     angle of v-direction in YZ plane
 * \param [in] vRes       resolution in v-direction
 */
Claus Kleinwort's avatar
Claus Kleinwort committed
309
310
311
GblDetectorLayer CreateLayerSit(const std::string aName, unsigned int layer,
		double xPos, double yPos, double zPos, double thickness, double uAngle,
		double uRes, double vAngle, double vRes) {
Claus Kleinwort's avatar
Claus Kleinwort committed
312
313
314
315
316
317
318
319
320
321
322
	Vector3d aCenter(xPos, yPos, zPos);
	Vector2d aResolution(uRes, vRes);
	Vector2d aPrecision(1. / (uRes * uRes), 1. / (vRes * vRes));
	Matrix3d measTrafo;
	const double cosU = cos(uAngle / 180. * M_PI);
	const double sinU = sin(uAngle / 180. * M_PI);
	const double cosV = cos(vAngle / 180. * M_PI);
	const double sinV = sin(vAngle / 180. * M_PI);
	measTrafo << 0., cosU, sinU, 0., cosV, sinV, 1., 0., 0.; // U,V,N
	Matrix3d alignTrafo;
	alignTrafo << 0., 1., 0., 0., 0., 1., 1., 0., 0.; // Y,Z,X
Claus Kleinwort's avatar
Claus Kleinwort committed
323
	return GblDetectorLayer(aName, layer, 2, thickness, aCenter, aResolution,
Claus Kleinwort's avatar
Claus Kleinwort committed
324
			aPrecision, measTrafo, alignTrafo);
325
326
327
}

}