GeneralBrokenLines  trunk_rev46
 All Classes Files Functions Variables Typedefs Pages
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 
8 #include <time.h>
9 #include "example1.h"
10 #include "TRandom3.h"
11 #include "GblTrajectory.h"
12 
13 TMatrixD gblSimpleJacobian(double ds, double cosl, double bfac) {
14  /* Simple jacobian: quadratic in arc length difference */
15  TMatrixD jac(5, 5);
16  jac.UnitMatrix();
17  jac[1][0] = -bfac * ds * cosl;
18  jac[3][0] = -0.5 * bfac * ds * ds * cosl;
19  jac[3][1] = ds;
20  jac[4][2] = ds;
21  return jac;
22 }
23 
24 void example1() {
25  /*
26  * Create points on initial trajectory, create trajectory from points,
27  * fit and write trajectory to MP-II binary file,
28  * get track parameter corrections and covariance matrix at points.
29  *
30  * Equidistant measurement layers and thin scatterers, propagation
31  * with simple jacobian (quadratic in arc length differences).
32  */
33 
34 //MP MilleBinary mille; // for producing MillePede-II binary file
35  unsigned int nTry = 10000; //10000; //: number of tries
36  unsigned int nLayer = 10; //: number of detector layers
37  unsigned int nTraj = 1; //: number of trajectories in composed trajectory
38  std::cout << " Gbltst $Rev: 45 $ " << nTry << ", " << nLayer << std::endl;
39 
40  TRandom *r = new TRandom3();
41 
42  clock_t startTime = clock();
43 // track direction
44  double sinLambda = 0.3;
45  double cosLambda = sqrt(1.0 - sinLambda * sinLambda);
46  double sinPhi = 0.;
47  double cosPhi = sqrt(1.0 - sinPhi * sinPhi);
48 // tDir = (cosLambda * cosPhi, cosLambda * sinPhi, sinLambda)
49 // U = Z x T / |Z x T|, V = T x U
50  TMatrixD uvDirT(3, 2);
51  uvDirT[0][0] = -sinPhi;
52  uvDirT[1][0] = cosPhi;
53  uvDirT[2][0] = 0.;
54  uvDirT[0][1] = -sinLambda * cosPhi;
55  uvDirT[1][1] = -sinLambda * sinPhi;
56  uvDirT[2][1] = cosLambda;
57 // measurement resolution
58  TVectorD measErr(2);
59  measErr[0] = 0.001;
60  measErr[1] = 0.001;
61  TVectorD measPrec(2); // (independent) precisions
62  measPrec[0] = 1.0 / (measErr[0] * measErr[0]);
63  measPrec[1] = 1.0 / (measErr[1] * measErr[1]);
64  TMatrixDSym measInvCov(2); // inverse covariance matrix
65  measInvCov.Zero();
66  measInvCov[0][0] = measPrec[0];
67  measInvCov[1][1] = measPrec[1];
68 // scattering error
69  TVectorD scatErr(2);
70  scatErr[0] = 0.001;
71  scatErr[1] = 0.001;
72  TVectorD scatPrec(2);
73  scatPrec[0] = 1.0 / (scatErr[0] * scatErr[0]);
74  scatPrec[1] = 1.0 / (scatErr[1] * scatErr[1]);
75 // (RMS of) CurvLinear track parameters
76  TVectorD clPar(5);
77  TVectorD clErr(5);
78  clErr[0] = 0.001;
79  clErr[1] = -0.1;
80  clErr[2] = 0.2;
81  clErr[3] = -0.15;
82  clErr[4] = 0.25;
83  TMatrixDSym clCov(5), clSeed(5);
84  unsigned int seedLabel = 0;
85 // additional parameters
86  TVectorD addPar(2);
87  addPar[0] = 0.0025;
88  addPar[1] = -0.005;
89 // global labels for MP
90  /*MP std::vector<int> globalLabels(2);
91  globalLabels[0] = 11;
92  globalLabels[1] = 12; */
93 
94  TMatrixD transMatrix(5, 5);
95  transMatrix.UnitMatrix();
96 
97  double bfac = 0.2998; // Bz*c for Bz=1
98  double step = 1.5 / cosLambda; // constant steps in RPhi
99 
100  double Chi2Sum = 0.;
101  int NdfSum = 0;
102  double LostSum = 0.;
103 
104  for (unsigned int iTry = 1; iTry <= nTry; ++iTry) {
105  // create list of list of points and transformation
106  std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > aPointsAndTransList;
107  aPointsAndTransList.reserve(nTraj);
108  std::vector<GblPoint> listOfPoints;
109  listOfPoints.reserve(2 * nLayer);
110  for (unsigned int iTraj = 0; iTraj < nTraj; ++iTraj) {
111  //std::vector<GblPoint> listOfPoints;
112  //listOfPoints.reserve(2 * nLayer);
113  // curvilinear track parameters
114  for (unsigned int i = 0; i < 5; ++i) {
115  clPar[i] = clErr[i] * r->Gaus();
116  }
117  clCov.Zero();
118  for (unsigned int i = 0; i < 5; ++i) {
119  clCov[i][i] = 1.0 * (clErr[i] * clErr[i]);
120  }
121 // std::cout << " Try " << iTry << ":" << clPar << std::endl;
122  TMatrixD addDer(2, 2);
123  addDer.Zero();
124  addDer[0][0] = 1.;
125  addDer[1][1] = 1.;
126 // arclength
127  double s = 0.;
128  std::vector<double> sPoint;
129  TMatrixD jacPointToPoint(5, 5);
130  jacPointToPoint.UnitMatrix();
131  for (unsigned int iLayer = 0; iLayer < nLayer; ++iLayer) {
132  // std::cout << " Layer " << iLayer << ", " << s << std::endl;
133 // measurement directions
134  double sinStereo = (iLayer % 2 == 0) ? 0. : 0.5;
135  double cosStereo = sqrt(1.0 - sinStereo * sinStereo);
136  TMatrixD mDir(2, 3);
137  mDir.Zero();
138  mDir[0][0] = sinStereo;
139  mDir[0][1] = cosStereo;
140  mDir[1][2] = 1.;
141 // projection local (uv) to measurement directions (dm/duv)
142  TMatrixD proL2m = mDir * uvDirT;
143 // point with (independent) measurements (in measurement system)
144  GblPoint point(jacPointToPoint);
145 // measurement - prediction in measurement system with error
146  TVectorD meas = proL2m * clPar.GetSub(3, 4);
147 //MP meas += addDer * addPar; // additional parameters
148  for (unsigned int i = 0; i < 2; ++i) {
149  meas[i] += measErr[i] * r->Gaus();
150  }
151  point.addMeasurement(proL2m, meas, measPrec);
152  /* point with (correlated) measurements (in local system)
153  GblPoint point(jacPointToPoint);
154  // projection measurement to local (uv) directions
155  TMatrixD proM2l = proL2m;
156  proM2l.Invert();
157  // measurement - prediction in local system with error
158  TVectorD meas(2);
159  for (unsigned int i = 0; i < 2; ++i) {
160  meas[i] = measErr[i] * r->Gaus() + addDer(i,0) * 0.0075;
161  }
162  meas = proM2l * meas + clPar.GetSub(3, 4);
163  TMatrixDSym localInvCov = measInvCov;
164  localInvCov.SimilarityT(proL2m);
165  point.addMeasurement(meas, localInvCov);
166 
167  // additional local parameters? */
168 // point.addLocals(addDer);
169 //MP point.addGlobals(globalLabels, addDer);
170  addDer *= -1.; // Der flips sign every measurement
171 // add point to trajectory
172  listOfPoints.push_back(point);
173  unsigned int iLabel = listOfPoints.size();
174  sPoint.push_back(s);
175  if (iLabel == seedLabel) {
176  clSeed = clCov.Invert();
177  }
178 // propagate to scatterer
179  jacPointToPoint = gblSimpleJacobian(step, cosLambda, bfac);
180  clPar = jacPointToPoint * clPar;
181  clCov = clCov.Similarity(jacPointToPoint);
182  s += step;
183  if (iLayer < nLayer - 1) {
184  TVectorD scat(2);
185  scat.Zero();
186  // point with scatterer
187  GblPoint point(jacPointToPoint);
188  point.addScatterer(scat, scatPrec);
189  listOfPoints.push_back(point);
190  iLabel = listOfPoints.size();
191  sPoint.push_back(s);
192  if (iLabel == seedLabel) {
193  clSeed = clCov.Invert();
194  }
195  // scatter a little
196  for (unsigned int i = 0; i < 2; ++i) {
197  clPar[i + 1] += scatErr[i] * r->Gaus();
198  clCov[i + 1] += scatErr[i] * scatErr[i];
199  }
200  // propagate to next measurement layer
201  clPar = jacPointToPoint * clPar;
202  clCov = clCov.Similarity(jacPointToPoint);
203  s += step;
204  }
205  }
206  //aPointsAndTransList.push_back(std::make_pair(listOfPoints, transMatrix));
207 //
208  }
209  // create trajectory
210  GblTrajectory traj(listOfPoints);
211  //GblTrajectory traj(listOfPoints, seedLabel, clSeed); // with external seed
212  //GblTrajectory traj(aPointsAndTransList); // composed
213 // fit trajectory
214  double Chi2;
215  int Ndf;
216  double lostWeight;
217  traj.fit(Chi2, Ndf, lostWeight);
218 // std::cout << " Fit: " << Chi2 << ", " << Ndf << ", " << lostWeight << std::endl;
219  /* TVectorD aCorrection(5);
220  TMatrixDSym aCovariance(5);
221  traj.getResults(1, aCorrection, aCovariance);
222  std::cout << " cor " << std::endl;
223  aCorrection.Print();
224  std::cout << " cov " << std::endl;
225  aCovariance.Print(); */
226 // write to MP binary file
227 //MP traj.milleOut(mille);
228  Chi2Sum += Chi2;
229  NdfSum += Ndf;
230  LostSum += lostWeight;
231  }
232  clock_t endTime = clock();
233  double diff = endTime - startTime;
234  double cps = CLOCKS_PER_SEC;
235  std::cout << " Time elapsed " << diff / cps << " s" << std::endl;
236  std::cout << " Chi2/Ndf = " << Chi2Sum / NdfSum << std::endl;
237 }
238 
Generated on Thu Sep 20 2012 15:45:38 for GeneralBrokenLines by   doxygen 1.8.1