12 unsigned int n_d,
unsigned int n_k,
13 unsigned int n_dzdl,
unsigned int n_z0,
16 HelixHough(n_phi, n_d, n_k, n_dzdl, n_z0, min_resolution, max_resolution, range),
18 vertex_sigma_xy(0.002),
19 vertex_sigma_z(0.005),
23 for(
unsigned int i=0;i<detrad.size();++i)
43 for(
unsigned int i=0;i<hits.size();i++)
46 if(
phis.find(hits[i].index) ==
phis.end() )
49 phis.insert(make_pair(hits[i].index,atan2(hits[i].
y,hits[i].
x)));
53 vector<double> chi2_hits;
55 temp_track.
hits.resize(4, hits[0]);
57 for(
unsigned int i1=0;i1<hits.size();++i1)
59 temp_track.
hits[0] = hits[i1];
61 for(
unsigned int i2=(i1+1);i2<hits.size();++i2)
63 if( hits[i2].
layer == hits[i1].
layer )
continue;
64 if( fabs(
phis[hits[i2].index] -
phis[hits[i1].index] ) >
M_PI_2 )
continue;
66 temp_track.
hits[1] = hits[i2];
68 for(
unsigned int i3=(i2+1);i3<hits.size();++i3)
70 if( hits[i3].layer == hits[i2].layer )
continue;
71 if( hits[i3].layer == hits[i1].layer )
continue;
72 if( fabs(
phis[hits[i3].index] -
phis[hits[i2].index] ) >
M_PI_2 )
continue;
73 if( fabs(
phis[hits[i3].index] -
phis[hits[i1].index] ) >
M_PI_2 )
continue;
75 temp_track.
hits[2] = hits[i3];
77 for(
unsigned int i4=(i3+1);i4<hits.size();++i4)
79 if( hits[i4].layer == hits[i3].layer )
continue;
80 if( hits[i4].layer == hits[i2].layer )
continue;
81 if( hits[i4].layer == hits[i1].layer )
continue;
82 if( fabs(
phis[hits[i4].index] -
phis[hits[i3].index] ) >
M_PI_2 )
continue;
83 if( fabs(
phis[hits[i4].index] -
phis[hits[i2].index] ) >
M_PI_2 )
continue;
84 if( fabs(
phis[hits[i4].index] -
phis[hits[i1].index] ) >
M_PI_2 )
continue;
86 temp_track.
hits[3] = hits[i4];
89 vector<unsigned int> tempcomb;
91 tempcomb[0] = temp_track.
hits[0].index;
92 tempcomb[1] = temp_track.
hits[1].index;
93 tempcomb[2] = temp_track.
hits[2].index;
94 tempcomb[3] = temp_track.
hits[3].index;
97 sort(tempcomb.begin(), tempcomb.end());
98 set<vector<unsigned int> >::iterator
it =
combos.find(tempcomb);
99 if(it !=
combos.end())
continue;
105 double chi2 =
fitTrack(temp_track, chi2_hits);
109 tracks.push_back(temp_track);
120 for(
unsigned int i=0;i<hits.size();i++)
123 if(
phis.find(hits[i].index) ==
phis.end() )
126 phis.insert(make_pair(hits[i].index,atan2(hits[i].
y,hits[i].
x)));
130 vector<double> chi2_hits;
132 temp_track.
hits.resize(5, hits[0]);
134 for(
unsigned int i1=0;i1<hits.size();++i1)
136 temp_track.
hits[0] = hits[i1];
138 for(
unsigned int i2=(i1+1);i2<hits.size();++i2)
140 if( hits[i2].
layer == hits[i1].
layer )
continue;
144 temp_track.
hits[1] = hits[i2];
146 for(
unsigned int i3=(i2+1);i3<hits.size();++i3)
148 if( hits[i3].layer == hits[i2].layer )
continue;
149 if( hits[i3].layer == hits[i1].layer )
continue;
154 temp_track.
hits[2] = hits[i3];
156 for(
unsigned int i4=(i3+1);i4<hits.size();++i4)
158 if( hits[i4].layer == hits[i3].layer )
continue;
159 if( hits[i4].layer == hits[i2].layer )
continue;
160 if( hits[i4].layer == hits[i1].layer )
continue;
166 temp_track.
hits[3] = hits[i4];
167 for(
unsigned int i5=(i4+1);i5<hits.size();++i5)
169 if( hits[i5].layer == hits[i4].layer )
continue;
170 if( hits[i5].layer == hits[i3].layer )
continue;
171 if( hits[i5].layer == hits[i2].layer )
continue;
172 if( hits[i5].layer == hits[i1].layer )
continue;
174 temp_track.
hits[4]= hits[i5];
177 vector<unsigned int> tempcomb;
178 tempcomb.assign(5,0);
179 tempcomb[0] = temp_track.
hits[0].index;
180 tempcomb[1] = temp_track.
hits[1].index;
181 tempcomb[2] = temp_track.
hits[2].index;
182 tempcomb[3] = temp_track.
hits[3].index;
183 tempcomb[4] = temp_track.
hits[4].index;
187 sort(tempcomb.begin(), tempcomb.end());
188 set<vector<unsigned int> >::iterator
it =
combos.find(tempcomb);
189 if(it !=
combos.end())
continue;
197 double chi2 =
fitTrack(temp_track, chi2_hits);
202 tracks.push_back(temp_track);
218 for(
unsigned int i=0;i<hits.size();i++)
221 if(
phis.find(hits[i].index) ==
phis.end() )
224 phis.insert(make_pair(hits[i].index,atan2(hits[i].
y,hits[i].
x)));
228 vector<double> chi2_hits;
230 temp_track.
hits.resize(6, hits[0]);
232 for(
unsigned int i1=0;i1<hits.size();++i1)
234 temp_track.
hits[0] = hits[i1];
236 for(
unsigned int i2=(i1+1);i2<hits.size();++i2)
238 if( hits[i2].
layer == hits[i1].
layer )
continue;
242 temp_track.
hits[1] = hits[i2];
244 for(
unsigned int i3=(i2+1);i3<hits.size();++i3)
246 if( hits[i3].layer == hits[i2].layer )
continue;
247 if( hits[i3].layer == hits[i1].layer )
continue;
251 temp_track.
hits[2] = hits[i3];
253 for(
unsigned int i4=(i3+1);i4<hits.size();++i4)
255 if( hits[i4].layer == hits[i3].layer )
continue;
256 if( hits[i4].layer == hits[i2].layer )
continue;
257 if( hits[i4].layer == hits[i1].layer )
continue;
262 temp_track.
hits[3] = hits[i4];
264 for(
unsigned int i5=(i4+1);i5<hits.size();++i5)
266 if( hits[i5].layer == hits[i4].layer )
continue;
267 if( hits[i5].layer == hits[i3].layer )
continue;
268 if( hits[i5].layer == hits[i2].layer )
continue;
269 if( hits[i5].layer == hits[i1].layer )
continue;
275 temp_track.
hits[4] = hits[i5];
277 for(
unsigned int i6=(i5+1);i6<hits.size();++i6)
279 if( hits[i6].layer == hits[i5].layer )
continue;
280 if( hits[i6].layer == hits[i4].layer )
continue;
281 if( hits[i6].layer == hits[i3].layer )
continue;
282 if( hits[i6].layer == hits[i2].layer )
continue;
283 if( hits[i6].layer == hits[i1].layer )
continue;
290 temp_track.
hits[5] = hits[i6];
293 vector<unsigned int> tempcomb;
294 tempcomb.assign(6,0);
295 tempcomb[0] = temp_track.
hits[0].index;
296 tempcomb[1] = temp_track.
hits[1].index;
297 tempcomb[2] = temp_track.
hits[2].index;
298 tempcomb[3] = temp_track.
hits[3].index;
299 tempcomb[4] = temp_track.
hits[4].index;
300 tempcomb[5] = temp_track.
hits[5].index;
303 sort(tempcomb.begin(), tempcomb.end());
304 set<vector<unsigned int> >::iterator
it =
combos.find(tempcomb);
314 double chi2 =
fitTrack(temp_track, chi2_hits);
320 tracks.push_back(temp_track);
332 unsigned int nt = input.size();
333 for(
unsigned int i=0;i<nt;++i)
335 output.push_back(input[i]);
348 chi2_hit.resize(track.
hits.size(), 0.);
350 MatrixXf
y = MatrixXf::Zero(track.
hits.size(), 1);
351 for(
unsigned int i=0;i<track.
hits.size();i++)
353 y(i, 0) = ( pow(track.
hits[i].x,2) + pow(track.
hits[i].y,2) );
358 MatrixXf
X = MatrixXf::Zero(track.
hits.size(), 3);
359 for(
unsigned int i=0;i<track.
hits.size();i++)
361 X(i, 0) = track.
hits[i].x;
362 X(i, 1) = track.
hits[i].y;
378 MatrixXf Xt = X.transpose();
380 MatrixXf prod = Xt*
X;
381 MatrixXf inv = prod.fullPivLu().inverse();
383 MatrixXf beta = inv*Xt*
y;
385 float cx = beta(0,0)*0.5;
386 float cy = beta(1,0)*0.5;
387 float r = sqrt(cx*cx + cy*cy - beta(2,0));
389 float phi = atan2(cy, cx);
390 float d = sqrt(cx*cx + cy*cy) -
r;
393 MatrixXf diff = y - (X*beta);
394 MatrixXf chi2 = (diff.transpose())*diff;
396 float dx = d*cos(phi);
397 float dy = d*sin(phi);
399 MatrixXf
y2 = MatrixXf::Zero(track.
hits.size(), 1);
400 for(
unsigned int i=0;i<track.
hits.size();i++)
402 y2(i,0) = track.
hits[i].z;
407 MatrixXf
X2 = MatrixXf::Zero(track.
hits.size(), 2);
408 for(
unsigned int i=0;i<track.
hits.size();i++)
410 float D = sqrt( pow(dx - track.
hits[i].x, 2) + pow(dy - track.
hits[i].y,2));
416 if(v >= 0.999999){v = 0.999999;}
447 MatrixXf Xt2 = X2.transpose();
448 MatrixXf prod2 = Xt2*
X2;
449 MatrixXf inv2 = prod2.fullPivLu().inverse();
450 MatrixXf beta2 = inv2*Xt2*
y2;
452 MatrixXf diff2 = y2 - (X2*beta2);
453 MatrixXf chi2_z = (diff2.transpose())*diff2;
455 float z0 = beta2(1,0);
456 float dzdl = beta2(0,0)/sqrt(1. + beta2(0,0)*beta2(0,0));
473 cx = (track.
d+
r)*cos(track.
phi);
474 cy = (track.
d+
r)*sin(track.
phi);
477 for(
unsigned int h=0;
h<track.
hits.size();
h++)
479 float dx1 = track.
hits[
h].x - cx;
480 float dy1 = track.
hits[
h].y - cy;
482 float dx2 = track.
hits[
h].x + cx;
483 float dy2 = track.
hits[
h].y + cy;
485 float xydiff1 = sqrt(dx1*dx1 + dy1*dy1) -
r;
486 float xydiff2 = sqrt(dx2*dx2 + dy2*dy2) -
r;
487 float xydiff = xydiff2;
488 if(fabs(xydiff1) < fabs(xydiff2)){ xydiff = xydiff1; }
497 chi2_hit[
h] += xydiff*xydiff/(ls_xy*ls_xy);
498 chi2_hit[
h] += diff2(
h,0)*diff2(
h,0);
500 chi2_tot += chi2_hit[
h];
503 unsigned int deg_of_freedom = 2*track.
hits.size() - 5;
507 track.
hits.pop_back();
511 return (chi2_tot)/((double)(deg_of_freedom));