GeneralBrokenLines V03-01-03
using EIGEN
exampleComposedKin.cpp
Go to the documentation of this file.
1/*
2 * exampleCdc.cpp
3 *
4 * Created on: 23 Mar 2023
5 * Author: kleinwrt
6 */
7
30#include <time.h>
31#include "exampleUtilCdc.h"
32#include "GblTrajectory.h"
33
34using namespace gbl;
35using namespace Eigen;
36
38
75
76 // detector setup, ~ Belle-II CDC
77 const unsigned int nSuper = 9; // number of super layers
78 double nLayer[nSuper] = { 8, 6, 6, 6, 6, 6, 6, 6, 6 }; // number of layers per super layer
79 double rInner[nSuper] = { 16.80, 25.70, 36.52, 47.69, 58.41, 69.53, 80.25,
80 91.37, 102.00 }; // inner radius of super layer
81 double rOuter[nSuper] = { 23.80, 34.80, 45.57, 56.69, 67.41, 78.53, 89.25,
82 100.37, 111.14 }; // outer radius of super layer
83 double stereo[nSuper] = { 0., 0.068, 0., -0.060, 0., 0.064, 0., -0.072, 0. }; // stereo angle of super layer
84 double zStart[nSuper] = { -35.9, -51.4, -57.5, -59.6, -61.8, -63.9, -66.0,
85 -68.2, -70.2 }; // -Z end of wires per super layer
86 double zEnd[nSuper] = { 67.9, 98.6, 132.9, 144.7, 146.8, 148.9, 151.0,
87 153.2, 155.3 }; // +Z end of wires per super layer
88
89 unsigned int nTry = 1000; //: number of tries
90 std::cout << " GblComposedKin $Id$ " << nTry << ", " << nSuper << std::endl;
91 srand(4711);
92 clock_t startTime = clock();
93
94 const double bfac = 0.0045; // B*c for 1.5 T
95
96 double beamPos[] = { 0., 0., 0. };
97 double beamSize[] = { 0.005, 0.005, 0.3 };
98 const bool useBeamSpot = true;
99
100 MilleBinary mille; // for producing MillePede-II binary file
101
102 double Chi2Sum = 0.;
103 int NdfSum = 0;
104 double LostSum = 0.;
105 int numFit = 0;
106 unsigned int iEvent = 0;
107
108 // event loop
109 for (unsigned int iTry = 0; iTry < nTry; ++iTry) {
110 unsigned int nTrack;
111 std::vector<std::array<double, 6>> allGenPar; // collect parameters for generation of tracks
112
113 // inline generation of track pair
114 iEvent++;
115 nTrack = 2;
116 // helix parameter for track generation
117 const double qbyp = 0.2; // 5 GeV
118 const double genDca = beamSize[0] * unrm(); // normal
119 const double genZ0 = beamSize[2] * unrm(); // normal
120 const double genPhi0 = M_PI * (2. * unif() - 1.); // uniform, [-30..30] deg
121 const double genDzds = 2. * unif() - 1.; // uniform, lambda ~ [-45..45] deg
122 const double genCurv = bfac * qbyp * sqrt(1. + genDzds * genDzds);
123
124 std::cout << " gen(inline) " << iEvent << ": " << genCurv << " "
125 << genPhi0 << " " << genDca << " " << genDzds << " " << genZ0
126 << std::endl;
127 std::array<double, 6> par = { genCurv, genPhi0, genDca, genDzds, genZ0,
128 1. };
129 allGenPar.push_back(par);
130 // back to back track
131 std::array<double, 6> par2 = { -genCurv, genPhi0 + M_PI, -genDca,
132 -genDzds, genZ0, -1. };
133 allGenPar.push_back(par2);
134
135 std::vector<std::vector<GblPoint>> listsOfPoints;
136 std::vector<Matrix<double, 5, 8>> listOfTrafos;
137 for (unsigned int iTrack = 0; iTrack < nTrack; ++iTrack) {
138 //
139 // generate hits: (virtual) wires as detector layers at constant radii
140 //
141 std::array<double, 6> genPar = allGenPar[iTrack];
142 double curv(genPar[0]), phi0(genPar[1]), dca(genPar[2]), dzds(
143 genPar[3]), z0(genPar[4]);
144 // local constant (Bfield) helix
145 GblSimpleHelix hlx = GblSimpleHelix(curv, phi0, dca, dzds, z0);
146
147 std::vector<GblDetectorLayer> layers;
148
149 // beam spot (required as first (and common) point for composed trajectory)
150 double sArc = hlx.getArcLengthXY(beamPos[0], beamPos[1]);
151 // add virtual layer
152 layers.push_back(
153 CreateImpactPar("impact", 0, beamPos[0], beamPos[1],
154 beamPos[2], phi0 + sArc * curv, dzds, beamSize[0],
155 beamSize[1], beamSize[2]));
156
157 unsigned int cLayer = 0;
158 for (unsigned int iSuper = 0; iSuper < nSuper; ++iSuper) {
159 double radius = rInner[iSuper];
160 double step = (rOuter[iSuper] - radius) / (nLayer[iSuper] - 1);
161
162 for (unsigned int iLayer = 0; iLayer < nLayer[iSuper];
163 ++iLayer) {
164 // check for |dca| < radius
165 if (fabs(dca) < radius) {
166 // arc length to layer
167 double sArc = hlx.getArcLengthR(radius);
168 if (sArc == 0.)
169 goto theEnd;
170 // phi of position at layer
171 double phiPos = hlx.getPhi(radius);
172 // (virtual) wire position
173 double xPos(cos(phiPos) * radius), yPos(
174 sin(phiPos) * radius), zPos(z0 + sArc * dzds);
175 // add virtual wire
176 if (zPos > zStart[iSuper] and zPos < zEnd[iSuper])
177 layers.push_back(
178 CreateWireCdc("wire", ++cLayer, xPos, yPos,
179 zPos, phi0 + sArc * curv, dzds,
180 stereo[iSuper], 0.015));
181 }
182 // next layer
183 radius += step;
184 }
185 }
186 theEnd:
187
188 //
189 // create GBL trajectory (list of GBL points)
190 //
191 // seed (with true parameters)
192 std::array<double, 6> seedPar = allGenPar[iTrack];
193 // transformations to common POSITION parameters (at vertex) of composed trajectory at first point
194 Matrix<double, 5, 8> innerTrafo = Matrix<double, 5, 8>::Zero();
195 if (allGenPar[iTrack][5] < 0.) {
196 // outer track
197 innerTrafo(0, 5) = -1.; // q/p
198 innerTrafo(1, 6) = 1.; // u'
199 innerTrafo(2, 7) = -1.; // v'
200 innerTrafo(3, 3) = -1.; // u
201 innerTrafo(4, 4) = 1.; // v
202 } else {
203 // leading track
204 innerTrafo(0, 0) = 1.; // q/p
205 innerTrafo(1, 1) = 1.; // u'
206 innerTrafo(2, 2) = 1.; // v'
207 innerTrafo(3, 3) = 1.; // u
208 innerTrafo(4, 4) = 1.; // v
209 }
210 listOfTrafos.push_back(innerTrafo);
211 // optionally distort seed
212 //seedPar[0] *= 0.99;
213 //seedPar[1] += 0.01;
214 double seedCurv(seedPar[0]), seedPhi0(seedPar[1]), seedDca(
215 seedPar[2]), seedDzds(seedPar[3]), seedZ0(seedPar[4]);
216 GblSimpleHelix seed = GblSimpleHelix(seedCurv, seedPhi0, seedDca,
217 seedDzds, seedZ0);
218 // (previous) arc-length
219 double sOld = 0.;
220 const double cosLambdaSeed = 1. / sqrt(1. + (seedDzds * seedDzds));
221 // for multiple scattering (not (yet) implemented)
222 const double qbyp = seedCurv / bfac * cosLambdaSeed;
223 // list of points on trajectory
224 std::vector<GblPoint> listOfPoints;
225 for (unsigned int iLayer = 0; iLayer < layers.size(); ++iLayer) {
226 GblDetectorLayer& layer = layers[iLayer];
227 unsigned int lid = layer.getLayerID(); // layer ID (0 = vertex)
228 // std::cout << " hit " << lid << " " << hits[iLayer].transpose() << std::endl;
229 // prediction from seeding helix
230 GblHelixPrediction pred = layer.intersectWithHelix(seed);
231 double sArc = pred.getArcLength(); // arc-length
232 Vector2d measPrediction = pred.getMeasPred(); // measurement prediction
233 Vector2d measPrecision = layer.getPrecision(); // measurement precision
234 // residuals
235 Vector2d res = measPrediction; // (virtual) wire positioned at measurement
236 // smear u according to resolution
237 Vector2d sigma = layer.getResolution();
238 if (lid > 0) // wire layer?
239 res[0] += sigma[0] * unrm();
240 else
241 std::cout << " impact par " << res[0] << " " << res[1]
242 << " " << seedCurv << " " << seedPhi0 << " "
243 << seedDca << " " << seedDzds << " " << seedZ0
244 << " " << sArc << " " << layers.size() << std::endl;
245 // transformation global system to local (curvilinear) (u,v) (matrix from row vectors)
246 Matrix<double, 2, 3> transG2l = pred.getCurvilinearDirs();
247 // transformation measurement system to global system
248 Matrix3d transM2g = layer.getMeasSystemDirs().inverse();
249 // projection matrix (measurement plane to local (u,v))
250 Matrix2d proM2l = transG2l * transM2g.block<3, 2>(0, 0); // skip measurement normal
251 // projection matrix (local (u,v) to measurement plane)
252 Matrix2d proL2m = proM2l.inverse();
253 // propagation
254 Matrix5d jacPointToPoint = gblSimpleJacobian(
255 (sArc - sOld) / cosLambdaSeed, cosLambdaSeed, bfac);
256 sOld = sArc;
257 // point with (independent) measurements (in measurement system)
258 GblPoint point(jacPointToPoint);
259 if (lid > 0 or (iTrack == 0 and useBeamSpot)) // vertex only for first track (large correlations!)
260 point.addMeasurement(proL2m, res, measPrecision);
261 // global labels and parameters for rigid body alignment
262 std::vector<int> labGlobal(6);
263 Vector3d pos = pred.getPosition();
264 Vector3d dir = pred.getDirection();
265 /* Layer alignment in local (measurement) system
266 for (int p = 0; p < 6; p++)
267 labGlobal[p] = (iLayer + 1) * 10 + p + 1;
268 Matrix<double, 2, 6> derGlobal = layer.getRigidBodyDerLocal(pos,
269 dir);
270 */
271 // Chamber alignment in global system
272 unsigned int layerID = layer.getLayerID();
273 for (int p = 0; p < 6; p++)
274 labGlobal[p] = layerID * 1000 + p + 1; // layer alignment
275 Matrix<double, 2, 6> derGlobal = layer.getRigidBodyDerGlobal(
276 pos, dir).block<2, 6>(0, 0);
277 point.addGlobals(labGlobal, derGlobal);
278 // add scatterer to point
279 double radlen = layer.getRadiationLength()
280 / fabs(pred.getCosIncidence());
281 double errMs = gblMultipleScatteringError(qbyp, radlen);// simple model
282 if (errMs > 0.) {
283 Vector2d scat(0., 0.);
284 Vector2d scatPrec(1. / (errMs * errMs),
285 1. / (errMs * errMs)); // scattering precision matrix is diagonal in curvilinear system
286 point.addScatterer(scat, scatPrec);
287 }
288 // add point to trajectory
289 listOfPoints.push_back(point);
290 }
291 listsOfPoints.push_back(listOfPoints);
292 }
293
294 // require 2 decent trajectories for composed trajectory
295 unsigned int nAccepted = 0;
296 for (unsigned int iTrack = 0; iTrack < nTrack; ++iTrack)
297 // sufficient length?
298 if (listsOfPoints[iTrack].size() >= 10)
299 nAccepted++;
300 if (nAccepted != 2)
301 continue;
302
303 // prepare composed trajectory
304 std::vector<std::pair<std::vector<GblPoint>, MatrixXd> > aPointsAndTransList;
305 for (unsigned int iTrack = 0; iTrack < nTrack; ++iTrack) {
306 aPointsAndTransList.push_back(
307 make_pair(listsOfPoints[iTrack], listOfTrafos[iTrack]));
308 }
309 //
310 // fit composed GBL trajectory
311 //
312
313 // create trajectory
314 GblTrajectory traj(aPointsAndTransList);
315 // fit trajectory
316 double Chi2;
317 int Ndf;
318 double lostWeight;
319 unsigned int ierr = traj.fit(Chi2, Ndf, lostWeight);
320 std::cout << " Fit " << iTry << ": " << Chi2 << ", " << Ndf << ", "
321 << lostWeight << std::endl;
322 traj.printTrajectory();
323 // successfully fitted?
324 if (!ierr) {
325 // write to MP binary file
326 traj.milleOut(mille);
327 // update statistics
328 Chi2Sum += Chi2;
329 NdfSum += Ndf;
330 LostSum += lostWeight;
331 numFit++;
332 }
333
334 }
335
336 clock_t endTime = clock();
337 double diff = endTime - startTime;
338 double cps = CLOCKS_PER_SEC;
339 std::cout << " Time elapsed " << diff / cps << " s" << std::endl;
340 std::cout << " Chi2/Ndf = " << Chi2Sum / NdfSum << std::endl;
341 std::cout << " Tracks fitted " << numFit << std::endl;
342 if (LostSum > 0.)
343 std::cout << " Weight lost " << LostSum << std::endl;
344}
345
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 getResolution() const
Get resolution.
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
double getArcLengthR(double aRadius) const
Get (2D) arc length for given radius (to ref. point)
double getPhi(double aRadius) const
Get phi (of point on circle) for given radius (to ref. point)
double getArcLengthXY(double xPos, double yPos) const
Get (2D) arc length for given point.
GBL trajectory.
Definition: GblTrajectory.h:50
void printTrajectory(unsigned int level=0) const
Print GblTrajectory.
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 exampleComposedKin()
Drift chamber example.
Definitions for exampleUtil(ities) for a Cylindrical Drift Chamber.
Namespace for the general broken lines package.
GblDetectorLayer CreateWireCdc(const std::string aName, unsigned int layer, double xPos, double yPos, double zPos, double phi, double tanLambda, double stereoAngle, double uRes)
Create a drift chamber wire with 1D measurement.
double unrm()
unit normal distribution, Box-Muller method, polar form
double gblMultipleScatteringError(double qbyp, double xbyx0)
Multiple scattering error.
GblDetectorLayer CreateImpactPar(const std::string aName, unsigned int layer, double xPos, double yPos, double zPos, double phi, double tanLambda, double xRes, double yRes, double zRes)
Create detector plane for impact parameters as 2D measurement.
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