11 using namespace Eigen;
15 FourHitSeedFinder::FourHitSeedFinder(vector<float>& detrad,
unsigned int n_phi,
unsigned int n_d,
unsigned int n_k,
unsigned int n_dzdl,
unsigned int n_z0,
HelixResolution& min_resolution,
HelixResolution& max_resolution,
HelixRange& range) :
HelixHough(n_phi, n_d, n_k, n_dzdl, n_z0, min_resolution, max_resolution, range), using_vertex(
false), vertex_sigma_xy(0.002), vertex_sigma_z(0.005), chi2_cut(3.)
17 for(
unsigned int i=0;i<detrad.size();++i)
27 cout <<
"findTracks" << endl;
86 vector<double> chi2_hits;
88 temp_track.
hits.resize(6, hits[0]);
90 for(
unsigned int i1=0;i1<hits.size();++i1)
92 temp_track.
hits[0] = hits[i1];
93 for(
unsigned int i2=(i1+1);i2<hits.size();++i2)
95 if( (hits[i2].
layer == hits[i1].
layer)){
continue;}
96 temp_track.
hits[1] = hits[i2];
97 for(
unsigned int i3=(i2+1);i3<hits.size();++i3)
99 if((hits[i3].layer == hits[i2].layer) || (hits[i3].layer == hits[i1].layer)){
continue;}
100 temp_track.
hits[2] = hits[i3];
101 for(
unsigned int i4=(i3+1);i4<hits.size();++i4)
103 if( (hits[i4].layer == hits[i3].layer) || (hits[i4].layer == hits[i2].layer) || (hits[i4].layer == hits[i1].layer)){
continue;}
104 temp_track.
hits[3] = hits[i4];
105 for(
unsigned int i5=(i4+1);i5<hits.size();++i5)
107 if( (hits[i5].layer == hits[i4].layer) || (hits[i5].layer == hits[i3].layer) || (hits[i5].layer == hits[i2].layer) || (hits[i5].layer == hits[i1].layer)){
continue;}
108 temp_track.
hits[4] = hits[i5];
109 for(
unsigned int i6=(i5+1);i6<hits.size();++i6)
111 if( (hits[i6].layer == hits[i5].layer) || (hits[i6].layer == hits[i4].layer) || (hits[i6].layer == hits[i3].layer) || (hits[i6].layer == hits[i2].layer) || (hits[i6].layer == hits[i1].layer)){
continue;}
112 temp_track.
hits[5] = hits[i6];
114 vector<unsigned int> tempcomb;
115 tempcomb.assign(6,0);
116 tempcomb[0] = temp_track.
hits[0].index;
117 tempcomb[1] = temp_track.
hits[1].index;
118 tempcomb[2] = temp_track.
hits[2].index;
119 tempcomb[3] = temp_track.
hits[3].index;
120 tempcomb[4] = temp_track.
hits[4].index;
121 tempcomb[5] = temp_track.
hits[5].index;
123 sort(tempcomb.begin(), tempcomb.end());
124 set<vector<unsigned int> >::iterator
it =
combos.find(tempcomb);
125 if(it !=
combos.end()){
continue;}
136 double chi2 =
fitTrack(temp_track, chi2_hits);
140 if(fabs(chi2) >
chi2_cut){
continue;}
141 tracks.push_back(temp_track);
153 unsigned int nt = input.size();
154 for(
unsigned int i=0;i<nt;++i)
156 output.push_back(input[i]);
158 cout<<
"# combinations = "<<
combos.size()<<endl;
169 bool swapped =
false;
170 if( fabs(track.
hits[0].x - track.
hits[1].x) < fabs(track.
hits[0].y - track.
hits[1].y) )
172 for(
unsigned int i=0;i<track.
hits.size();i++)
183 chi2_hit.resize(track.
hits.size(), 0.);
185 MatrixXf
y = MatrixXf::Zero(track.
hits.size(), 1);
186 for(
unsigned int i=0;i<track.
hits.size();i++)
188 y(i, 0) = track.
hits[i].y;
193 MatrixXf
X = MatrixXf::Zero(track.
hits.size(), 2);
194 for(
unsigned int i=0;i<track.
hits.size();i++)
197 X(i, 1) = track.
hits[i].x;
210 MatrixXf Xt = X.transpose();
212 MatrixXf prod = Xt*
X;
213 MatrixXf inv = prod.fullPivLu().inverse();
215 MatrixXf beta = inv*Xt*
y;
217 float a1 = beta(1,0);
218 float phi = 0.5*
M_PI - atan(-a1);
221 for(
unsigned int i=0;i<track.
hits.size();i++)
229 float d = track.
hits[0].x*cos(phi) + track.
hits[0].y*sin(phi);
232 MatrixXf diff = y - (X*beta);
233 MatrixXf chi2 = (diff.transpose())*diff;
235 float dx = d*cos(phi);
236 float dy = d*sin(phi);
238 MatrixXf
y2 = MatrixXf::Zero(track.
hits.size(), 1);
239 for(
unsigned int i=0;i<track.
hits.size();i++)
241 y2(i,0) = track.
hits[i].z;
246 MatrixXf
X2 = MatrixXf::Zero(track.
hits.size(), 2);
247 for(
unsigned int i=0;i<track.
hits.size();i++)
249 float D = sqrt( pow(dx - track.
hits[i].x, 2) + pow(dy - track.
hits[i].y,2));
266 MatrixXf Xt2 = X2.transpose();
267 MatrixXf prod2 = Xt2*
X2;
268 MatrixXf inv2 = prod2.fullPivLu().inverse();
269 MatrixXf beta2 = inv2*Xt2*
y2;
271 MatrixXf diff2 = y2 - (X2*beta2);
272 MatrixXf chi2_z = (diff2.transpose())*diff2;
274 float z0 = beta2(1,0);
275 float dzdl = beta2(0,0)/sqrt(1. + beta2(0,0)*beta2(0,0));
285 for(
unsigned int h=0;
h<track.
hits.size();
h++)
288 chi2_hit[
h] += diff(
h,0)*diff(
h,0);
289 chi2_hit[
h] += diff2(
h,0)*diff2(
h,0);
291 chi2_tot += chi2_hit[
h];
294 unsigned int deg_of_freedom = 2*track.
hits.size() - 5;
298 track.
hits.pop_back();
302 return (chi2_tot)/((double)(deg_of_freedom));
314 chi2_hit.resize(track.
hits.size(), 0.);
316 MatrixXf
y = MatrixXf::Zero(track.
hits.size(), 1);
317 for(
unsigned int i=0;i<track.
hits.size();i++)
319 y(i, 0) = ( pow(track.
hits[i].x,2) + pow(track.
hits[i].y,2) );
324 MatrixXf
X = MatrixXf::Zero(track.
hits.size(), 3);
325 for(
unsigned int i=0;i<track.
hits.size();i++)
327 X(i, 0) = track.
hits[i].x;
328 X(i, 1) = track.
hits[i].y;
344 MatrixXf Xt = X.transpose();
346 MatrixXf prod = Xt*
X;
349 MatrixXf beta = prod.ldlt().solve(Xty);
351 float cx = beta(0,0)*0.5;
352 float cy = beta(1,0)*0.5;
353 float r = sqrt(cx*cx + cy*cy - beta(2,0));
355 float phi = atan2(cy, cx);
356 float d = sqrt(cx*cx + cy*cy) -
r;
359 MatrixXf diff = y - (X*beta);
360 MatrixXf chi2 = (diff.transpose())*diff;
362 float dx = d*cos(phi);
363 float dy = d*sin(phi);
365 MatrixXf
y2 = MatrixXf::Zero(track.
hits.size(), 1);
366 for(
unsigned int i=0;i<track.
hits.size();i++)
368 y2(i,0) = track.
hits[i].z;
373 MatrixXf
X2 = MatrixXf::Zero(track.
hits.size(), 2);
374 for(
unsigned int i=0;i<track.
hits.size();i++)
376 float D = sqrt( pow(dx - track.
hits[i].x, 2) + pow(dy - track.
hits[i].y,2));
382 if(v >= 0.999999){v = 0.999999;}
413 MatrixXf Xt2 = X2.transpose();
414 MatrixXf prod2 = Xt2*
X2;
416 MatrixXf Xty2 = Xt2*
y2;
417 MatrixXf beta2 = prod2.ldlt().solve(Xty2);
420 MatrixXf diff2 = y2 - (X2*beta2);
421 MatrixXf chi2_z = (diff2.transpose())*diff2;
423 float z0 = beta2(1,0);
424 float dzdl = beta2(0,0)/sqrt(1. + beta2(0,0)*beta2(0,0));
441 cx = (track.
d+
r)*cos(track.
phi);
442 cy = (track.
d+
r)*sin(track.
phi);
445 for(
unsigned int h=0;
h<track.
hits.size();
h++)
447 float dx1 = track.
hits[
h].x - cx;
448 float dy1 = track.
hits[
h].y - cy;
450 float dx2 = track.
hits[
h].x + cx;
451 float dy2 = track.
hits[
h].y + cy;
453 float xydiff1 = sqrt(dx1*dx1 + dy1*dy1) -
r;
454 float xydiff2 = sqrt(dx2*dx2 + dy2*dy2) -
r;
455 float xydiff = xydiff2;
456 if(fabs(xydiff1) < fabs(xydiff2)){ xydiff = xydiff1; }
465 chi2_hit[
h] += xydiff*xydiff/(ls_xy*ls_xy);
466 chi2_hit[
h] += diff2(
h,0)*diff2(
h,0);
468 chi2_tot += chi2_hit[
h];
471 unsigned int deg_of_freedom = 2*track.
hits.size() - 5;
475 track.
hits.pop_back();
479 return (chi2_tot)/((double)(deg_of_freedom));