60 unsigned int nTry = 1000;
61 unsigned int nLayer = 10;
62 std::cout <<
" Gbltst-eigen $Id$ " << nTry <<
", " << nLayer << std::endl;
66 clock_t startTime = clock();
68 double sinLambda = 0.3;
69 double cosLambda = sqrt(1.0 - sinLambda * sinLambda);
71 double cosPhi = sqrt(1.0 - sinPhi * sinPhi);
74 Matrix<double, 2, 3> uvDir;
75 uvDir(0, 0) = -sinPhi;
78 uvDir(1, 0) = -sinLambda * cosPhi;
79 uvDir(1, 1) = -sinLambda * sinPhi;
80 uvDir(1, 2) = cosLambda;
83 measErr << 0.001, 0.001;
85 measPrec << 1.0 / (measErr(0) * measErr(0)), 1.0 / (measErr(1) * measErr(1));
88 measInvCov(0, 0) = measPrec[0];
89 measInvCov(1, 1) = measPrec[1];
92 scatErr << 0.001, 0.001;
94 scatPrec << 1.0 / (scatErr(0) * scatErr(0)), 1.0 / (scatErr(1) * scatErr(1));
98 clErr << 0.001, -0.1, 0.2, -0.15, 0.25;
100 unsigned int seedLabel = 99999;
103 addPar << 0.0025, -0.005;
104 std::vector<int> globalLabels;
105 globalLabels.push_back(4711);
106 globalLabels.push_back(4712);
112 double bfac = 0.2998;
113 double step = 1.5 / cosLambda;
120 for (
unsigned int iTry = 1; iTry <= nTry; ++iTry) {
122 for (
unsigned int i = 0; i < 5; ++i) {
123 clPar[i] = clErr[i] *
unrm();
126 for (
unsigned int i = 0; i < 5; ++i) {
127 clCov(i, i) = 1.0 * (clErr[i] * clErr[i]);
137 jacPointToPoint.setIdentity();
139 std::vector<GblPoint> listOfPoints;
140 listOfPoints.reserve(2 * nLayer);
142 for (
unsigned int iLayer = 0; iLayer < nLayer; ++iLayer) {
145 double sinStereo = (iLayer % 2 == 0) ? 0. : 0.1;
146 double cosStereo = sqrt(1.0 - sinStereo * sinStereo);
147 Matrix<double, 3, 2> mDirT;
149 mDirT(1, 0) = cosStereo;
150 mDirT(2, 0) = sinStereo;
151 mDirT(1, 1) = -sinStereo;
152 mDirT(2, 1) = cosStereo;
154 Matrix2d proM2l = uvDir * mDirT;
156 Matrix2d proL2m = proM2l.inverse();
158 GblPoint pointMeas(jacPointToPoint);
160 Vector2d meas = proL2m * clPar.tail(2);
162 for (
unsigned int i = 0; i < 2; ++i) {
163 meas[i] += measErr[i] *
unrm();
182 listOfPoints.push_back(pointMeas);
183 unsigned int iLabel = listOfPoints.size();
184 if (iLabel == seedLabel) {
185 clSeed = clCov.inverse();
190 clPar = jacPointToPoint * clPar;
191 clCov = jacPointToPoint * clCov * jacPointToPoint.transpose();
193 if (iLayer < nLayer - 1) {
194 Vector2d scat(0., 0.);
196 GblPoint pointScat(jacPointToPoint);
198 listOfPoints.push_back(pointScat);
199 iLabel = listOfPoints.size();
200 if (iLabel == seedLabel) {
201 clSeed = clCov.inverse();
204 for (
unsigned int i = 0; i < 2; ++i) {
205 clPar[i + 1] += scatErr[i] *
unrm();
206 clCov(i + 1, i + 1) += scatErr[i] * scatErr[i];
209 clPar = jacPointToPoint * clPar;
210 clCov = jacPointToPoint * clCov * jacPointToPoint.transpose();
228 traj.
fit(Chi2, Ndf, lostWeight);
230 std::cout <<
" Fit: " << Chi2 <<
", " << Ndf <<
", " << lostWeight
231 <<
", " << bandCond << std::endl;
261 LostSum += lostWeight;
265 clock_t endTime = clock();
266 double diff = endTime - startTime;
267 double cps = CLOCKS_PER_SEC;
268 std::cout <<
" Time elapsed " << diff / cps <<
" s" << std::endl;
269 std::cout <<
" Chi2/Ndf = " << Chi2Sum / NdfSum << std::endl;
270 std::cout <<
" Tracks fitted " << numFit << std::endl;
272 std::cout <<
" Weight lost " << LostSum << std::endl;
GblTrajectory definition.
Definitions for GBL utilities.
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.
void addScatterer(const Eigen::Vector2d &aResiduals, const Eigen::Matrix2d &aPrecision)
Add a thin scatterer to a point.
double getBandCondition() const
Get condition from band (decomposition)
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight, const std::string &optionList="", unsigned int aLabel=0)
Perform fit of (valid) trajectory.
Namespace for the general broken lines package.
double unrm()
unit normal distribution, Box-Muller method, polar form
Eigen::Matrix< double, 5, 1 > Vector5d
Matrix5d gblSimpleJacobian(double ds, double cosl, double bfac)
Simple jacobian.
Eigen::Matrix< double, 5, 5 > Matrix5d