GeneralBrokenLines V03-01-01
using EIGEN
exampleDc.cpp
Go to the documentation of this file.
1/*
2 * exampleDc.cpp
3 *
4 * Created on: 6 Nov 2018
5 * Author: kleinwrt
6 */
7
30#include <time.h>
31#include "exampleDc.h"
32#include "GblTrajectory.h"
33
34using namespace gbl;
35using namespace Eigen;
36
38
73void exampleDc() {
74
75 // detector layers (ordered in Z):
76 // name, position (x,y,z), thickness (X/X_0), xz-angle, stereo-angle, resolution
77 double thickness[12];
78 thickness[0] = 0.0025; // gas, wire, wall
79 for (unsigned int iLayer = 1; iLayer < 11; ++iLayer)
80 thickness[iLayer] = 0.0005; // gas, wire
81 thickness[11] = 0.0025; // gas, wire, wall
82 // list of layers
83 std::vector<GblDetectorLayer> layers;
84 const double theta(25.), cosTheta(cos(theta / 180. * M_PI)), sinTheta(
85 sin(theta / 180. * M_PI));
86 // 1. chamber at distance of 200 cm
87 double dist = 200.;
88 for (unsigned int iLayer = 0; iLayer < 6; ++iLayer) {
89 layers.push_back(
90 CreateLayerDc("CH1+", 1, dist * sinTheta, 0.0, dist * cosTheta,
91 thickness[iLayer], theta, 6., 0.030)); // +6 deg stereo layers
92 dist += 2.;
93 }
94 for (unsigned int iLayer = 6; iLayer < 12; ++iLayer) {
95 layers.push_back(
96 CreateLayerDc("CH1-", 1, dist * sinTheta, 0.0, dist * cosTheta,
97 thickness[iLayer], theta, -6., 0.030)); // -6 deg stereo layers
98 dist += 2.;
99 }
100 // 2. chamber at distance of 300 cm
101 dist = 300.;
102 for (unsigned int iLayer = 0; iLayer < 6; ++iLayer) {
103 layers.push_back(
104 CreateLayerDc("CH2+", 2, dist * sinTheta, 0.0, dist * cosTheta,
105 thickness[iLayer], theta, 6., 0.030)); // +6 deg stereo layers
106 dist += 2.;
107 }
108 for (unsigned int iLayer = 6; iLayer < 12; ++iLayer) {
109 layers.push_back(
110 CreateLayerDc("CH2-", 2, dist * sinTheta, 0.0, dist * cosTheta,
111 thickness[iLayer], theta, -6., 0.030)); // -6 deg stereo layers
112 dist += 2.;
113 }
114 // 3. chamber at distance of 400 cm
115 dist = 400.;
116 for (unsigned int iLayer = 0; iLayer < 6; ++iLayer) {
117 layers.push_back(
118 CreateLayerDc("CH3+", 3, dist * sinTheta, 0.0, dist * cosTheta,
119 thickness[iLayer], theta, 6., 0.030)); // +6 deg stereo layers
120 dist += 2.;
121 }
122 for (unsigned int iLayer = 6; iLayer < 12; ++iLayer) {
123 layers.push_back(
124 CreateLayerDc("CH3-", 3, dist * sinTheta, 0.0, dist * cosTheta,
125 thickness[iLayer], theta, -6., 0.030)); // -6 deg stereo layers
126 dist += 2.;
127 }
128
129 /* print layers
130 for (unsigned int iLayer = 0; iLayer < layers.size(); ++iLayer) {
131 layers[iLayer].print();
132 } */
133
134 unsigned int nTry = 10000; //: number of tries
135 std::cout << " GblDc $Id$ " << nTry << ", " << layers.size()
136 << std::endl;
137 srand(4711);
138 clock_t startTime = clock();
139
140 double qbyp = 0.2; // 5 GeV
141 // const double bfac = 0.003; // B*c for 1 T
142 const double bfac = 0.; // B*c for 0 T
143
144 MilleBinary mille; // for producing MillePede-II binary file
145
146 double Chi2Sum = 0.;
147 int NdfSum = 0;
148 double LostSum = 0.;
149 int numFit = 0;
150
151 for (unsigned int iTry = 0; iTry < nTry; ++iTry) {
152
153 // helix parameter for track generation
154 const double genDca = 0.1 * unrm(); // normal
155 const double genZ0 = 0.1 * unrm(); // normal
156 const double genPhi0 = 0.52 * (2. * unif() - 1.); // uniform, [-30..30] deg
157 const double genDzds = 10. * unif() + 1.2; // uniform, lambda ~ [50..85] deg
158 const double genCurv = bfac * qbyp * sqrt(1. + genDzds * genDzds);
159
160 //
161 // generate hits
162 //
163 std::vector<Vector2d> hits;
164 double curv(genCurv), phi0(genPhi0), dca(genDca), dzds(genDzds), z0(
165 genZ0);
166 const double cosLambda = 1. / sqrt(1. + dzds * dzds);
167 for (unsigned int iLayer = 0; iLayer < layers.size(); ++iLayer) {
168 // local constant (Bfield) helix
169 GblSimpleHelix hlx = GblSimpleHelix(curv, phi0, dca, dzds, z0);
170 // prediction from local helix
171 GblHelixPrediction pred = layers[iLayer].intersectWithHelix(hlx);
172 // std::cout << " layer " << iLayer << " arc-length " << pred.getArcLength() << std::endl;
173 Vector2d meas = pred.getMeasPred();
174 // smear according to resolution
175 Vector2d sigma = layers[iLayer].getResolution();
176 meas[0] += sigma[0] * unrm();
177 meas[1] += sigma[1] * unrm();
178 // save hit
179 hits.push_back(meas);
180 // scatter at intersection point
181 Vector3d measPos = pred.getPosition();
182 double radlen = layers[iLayer].getRadiationLength()
183 / fabs(pred.getCosIncidence());
184 double errMs = gblMultipleScatteringError(qbyp, radlen); // simple model
185 // move to intersection point
186 hlx.moveToXY(measPos[0], measPos[1], phi0, dca, z0); // update phi0, dca, z0
187 phi0 += unrm() * errMs / cosLambda; // scattering for phi
188 dzds += unrm() * errMs / (cosLambda * cosLambda); // scattering for dzds
189 GblSimpleHelix newhlx = GblSimpleHelix(curv, phi0, dca, dzds, z0); // after scattering
190 // move back
191 newhlx.moveToXY(-measPos[0], -measPos[1], phi0, dca, z0); // update phi0, dca, z0
192 }
193
194 //
195 // fit track with GBL
196 //
197 // seed (with true parameters)
198 double seedCurv(genCurv), seedPhi0(genPhi0), seedDca(genDca), seedDzds(
199 genDzds), seedZ0(genZ0);
200 GblSimpleHelix seed = GblSimpleHelix(seedCurv, seedPhi0, seedDca,
201 seedDzds, seedZ0);
202 // (previous) arc-length
203 double sOld = 0.;
204 const double cosLambdaSeed = 1. / sqrt(1. + (seedDzds * seedDzds));
205 // list of points on trajectory
206 std::vector<GblPoint> listOfPoints;
207 for (unsigned int iLayer = 0; iLayer < layers.size(); ++iLayer) {
208 // std::cout << " hit " << iLayer << " " << hits[iLayer].transpose() << std::endl;
209 GblDetectorLayer& layer = layers[iLayer];
210 // prediction from seeding helix
211 GblHelixPrediction pred = layer.intersectWithHelix(seed);
212 double sArc = pred.getArcLength(); // arc-length
213 Vector2d measPrediction = pred.getMeasPred(); // measurement prediction
214 Vector2d measPrecision = layer.getPrecision(); // measurement precision
215 // residuals
216 Vector2d res(hits[iLayer][0] - measPrediction[0],
217 hits[iLayer][1] - measPrediction[1]);
218 // transformation global system to local (curvilinear) (u,v) (matrix from row vectors)
219 Matrix<double, 2, 3> transG2l = pred.getCurvilinearDirs();
220 // transformation measurement system to global system
221 Matrix3d transM2g = layer.getMeasSystemDirs().inverse();
222 // projection matrix (measurement plane to local (u,v))
223 Matrix2d proM2l = transG2l * transM2g.block<3, 2>(0, 0);// skip measurement normal
224 // projection matrix (local (u,v) to measurement plane)
225 Matrix2d proL2m = proM2l.inverse();
226 // propagation
227 Matrix5d jacPointToPoint = gblSimpleJacobian(
228 (sArc - sOld) / cosLambdaSeed, cosLambdaSeed, bfac);
229 sOld = sArc;
230 // point with (independent) measurements (in measurement system)
231 GblPoint point(jacPointToPoint);
232 point.addMeasurement(proL2m, res, measPrecision);
233 // global labels and parameters for rigid body alignment
234 std::vector<int> labGlobal(6);
235 Vector3d pos = pred.getPosition();
236 Vector3d dir = pred.getDirection();
237 /* Layer alignment in local (measurement) system
238 for (int p = 0; p < 6; p++)
239 labGlobal[p] = (iLayer + 1) * 10 + p + 1;
240 Matrix<double, 2, 6> derGlobal = layer.getRigidBodyDerLocal(pos,
241 dir); */
242 // Chamber alignment in global system (as common system for both stereo orientations)
243 unsigned int layerID = layer.getLayerID();
244 for (int p = 0; p < 6; p++)
245 labGlobal[p] = layerID * 1000 + p + 1; // chamber alignment
246 Matrix<double, 2, 6> derGlobal = layer.getRigidBodyDerGlobal(pos,
247 dir).block<2, 6>(0, 0);
248 point.addGlobals(labGlobal, derGlobal);
249 // add scatterer to point
250 double radlen = layer.getRadiationLength()
251 / fabs(pred.getCosIncidence());
252 double errMs = gblMultipleScatteringError(qbyp, radlen);// simple model
253 if (errMs > 0.) {
254 Vector2d scat(0., 0.);
255 Vector2d scatPrec(1. / (errMs * errMs), 1. / (errMs * errMs)); // scattering precision matrix is diagonal in curvilinear system
256 point.addScatterer(scat, scatPrec);
257 }
258 // add point to trajectory
259 listOfPoints.push_back(point);
260 }
261 // create trajectory
262 GblTrajectory traj(listOfPoints, bfac != 0.);
263 // fit trajectory
264 double Chi2;
265 int Ndf;
266 double lostWeight;
267 unsigned int ierr = traj.fit(Chi2, Ndf, lostWeight);
268 //std::cout << " Fit " << iTry << ": "<< Chi2 << ", " << Ndf << ", " << lostWeight << std::endl;
269 // successfully fitted?
270 if (!ierr) {
271 // write to MP binary file
272 traj.milleOut(mille);
273 // update statistics
274 Chi2Sum += Chi2;
275 NdfSum += Ndf;
276 LostSum += lostWeight;
277 numFit++;
278 }
279 }
280 clock_t endTime = clock();
281 double diff = endTime - startTime;
282 double cps = CLOCKS_PER_SEC;
283 std::cout << " Time elapsed " << diff / cps << " s" << std::endl;
284 std::cout << " Chi2/Ndf = " << Chi2Sum / NdfSum << std::endl;
285 std::cout << " Tracks fitted " << numFit << std::endl;
286 if (LostSum > 0.)
287 std::cout << " Weight lost " << LostSum << std::endl;
288}
289
290namespace gbl {
291
293
305GblDetectorLayer CreateLayerDc(const std::string aName, unsigned int layer,
306 double xPos, double yPos, double zPos, double thickness, double xzAngle,
307 double stereoAngle, double uRes) {
308 Vector3d aCenter(xPos, yPos, zPos);
309 Vector2d aResolution(uRes, 0.);
310 Vector2d aPrecision(1. / (uRes * uRes), 0.);
311 Matrix3d measTrafo;
312 const double cosXz = cos(xzAngle / 180. * M_PI);
313 const double sinXz = sin(xzAngle / 180. * M_PI);
314 const double cosSt = cos(stereoAngle / 180. * M_PI);
315 const double sinSt = sin(stereoAngle / 180. * M_PI);
316 measTrafo << cosSt * cosXz, sinSt, -cosSt * sinXz, -sinSt * cosXz, cosSt, sinSt
317 * sinXz, sinXz, 0., cosXz; // U,V,N
318 Matrix3d alignTrafo;
319 alignTrafo << cosXz, 0., -sinXz, 0., 1., 0., sinXz, 0., cosXz; // I,J,K
320 return GblDetectorLayer(aName, layer, 1, thickness, aCenter, aResolution,
321 aPrecision, measTrafo, alignTrafo);
322}
323
324}
GblTrajectory definition.
Detector layer.
Definition: GblUtilities.h:107
GblHelixPrediction intersectWithHelix(GblSimpleHelix hlx) const
Intersect with helix.
double getRadiationLength() const
Get radiation length.
unsigned int getLayerID() const
Get layer ID.
Eigen::Matrix3d getMeasSystemDirs() const
Get directions of measurement system.
Eigen::Matrix< double, 3, 6 > getRigidBodyDerGlobal(Eigen::Vector3d &position, Eigen::Vector3d &direction) const
Get rigid body derivatives in global frame.
Eigen::Vector2d getPrecision() const
Get precision.
Prediction on helix.
Definition: GblUtilities.h:49
Eigen::Matrix< double, 2, 3 > getCurvilinearDirs() const
Get curvilinear directions (U,V)
double getArcLength() const
Get arc-length.
const Eigen::Vector2d & getMeasPred() const
Get (measurement) prediction.
double getCosIncidence() const
Get cosine of incidence.
const Eigen::Vector3d & getPosition() const
Get position.
const Eigen::Vector3d & getDirection() const
Get position.
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 addGlobals(const std::vector< int > &aLabels, const Eigen::MatrixBase< Derivative > &aDerivatives)
Add global derivatives to a point.
Definition: GblPoint.h:293
void addScatterer(const Eigen::Vector2d &aResiduals, const Eigen::Matrix2d &aPrecision)
Add a thin scatterer to a point.
Definition: GblPoint.cpp:159
Simple helix.
Definition: GblUtilities.h:78
void moveToXY(double xPos, double yPos, double &newPhi0, double &newDca, double &newZ0) const
Move to new reference point (X,Y)
GBL trajectory.
Definition: GblTrajectory.h:50
void milleOut(MilleBinary &aMille)
Write valid trajectory to Millepede-II binary file.
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight, const std::string &optionList="", unsigned int aLabel=0)
Perform fit of (valid) trajectory.
Millepede-II (binary) record.
Definition: MilleBinary.h:81
void exampleDc()
Drift chamber example.
Definition: exampleDc.cpp:73
Definitions for exampleDc.
Namespace for the general broken lines package.
double unrm()
unit normal distribution, Box-Muller method, polar form
double gblMultipleScatteringError(double qbyp, double xbyx0)
Multiple scattering error.
GblDetectorLayer CreateLayerDc(const std::string aName, unsigned int layer, double xPos, double yPos, double zPos, double thickness, double xzAngle, double stereoAngle, double uRes)
Create a drift chamber layer with 1D measurement.
Definition: exampleDc.cpp:305
double unif()
uniform distribution [0..1]
Matrix5d gblSimpleJacobian(double ds, double cosl, double bfac)
Simple jacobian.
Eigen::Matrix< double, 5, 5 > Matrix5d
Definition: GblData.h:46