10 using namespace Eigen;
13 ThreeHitSeedGrower::ThreeHitSeedGrower(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.)
15 unsigned int n_layers = detrad.size();
25 {cout<<
"findTracks::entered with "<<hits.size()<<
" hits"<<endl;
26 vector<double> chi2_hits;
29 temp_track.
hits.resize(3, hits[0]);
30 for(
unsigned int i1=0;i1<hits.size();++i1)
32 temp_track.
hits[0] = hits[i1];
33 for(
unsigned int i2=(i1+1);i2<hits.size();++i2)
35 if( (hits[i2].
layer == hits[i1].
layer)){
continue;}
36 temp_track.
hits[1] = hits[i2];
37 for(
unsigned int i3=(i2+1);i3<hits.size();++i3)
39 if((hits[i3].layer == hits[i2].layer) || (hits[i3].layer == hits[i1].layer)){
continue;}
40 temp_track.
hits[2] = hits[i3];
42 chi2 =
fitTrack(temp_track, chi2_hits);
44 vector<unsigned int> nhit_layer;
48 nhit_layer[temp_track.
hits[
ihit].layer]+=1;
50 if(
GrowTrack(temp_track, nhit_layer, hits, tracks, 0) ==
false)
52 vector<unsigned int> tempcomb;
54 tempcomb[0] = temp_track.
hits[0].index;
55 tempcomb[1] = temp_track.
hits[1].index;
56 tempcomb[2] = temp_track.
hits[2].index;
57 sort(tempcomb.begin(), tempcomb.end());
58 set<vector<unsigned int> >::iterator
it =
combos.find(tempcomb);
59 if(it !=
combos.end()){
continue;}
61 chi2 =
fitTrack(temp_track, chi2_hits);
62 tracks.push_back(temp_track);
63 cout<<
"findTrack::added a 3 hit track"<<endl;
69 cout<<
"leaving findTrack"<<endl;
79 unsigned int init_size = seed_track.
hits.size();
80 vector<double> chi2_hits;
83 while(c_hit < hits.size())
85 hit_added =
addOneHit(seed_track, nhit_layer, c_hit, hits);
87 if(hit_added ==
false)
89 unsigned int final_size = seed_track.
hits.size();
92 if(
GrowTrack(seed_track, nhit_layer, hits, tracks, c_hit) ==
false)
95 chi2 =
fitTrack(seed_track,chi2_hits);
98 vector<unsigned int> tempcomb;
99 tempcomb.assign(seed_track.
hits.size(),0);
102 sort(tempcomb.begin(), tempcomb.end());
103 set<vector<unsigned int> >::iterator
it =
combos.find(tempcomb);
105 seed_track.
hits.pop_back();
109 tracks.push_back(seed_track);
110 cout<<
"GrowTrack::added a track of size: "<<seed_track.
hits.size()<<endl;
121 if(seed_track.
hits[
ihit].index == hits[c_hit].index)
126 if(nhit_layer[hits[c_hit].
layer] == 0 || (nhit_layer[hits[c_hit].
layer] == 1 && hits[c_hit].
layer >1)){
127 vector<double> chi2_hits;
129 seed_track.
hits.push_back(hits[c_hit]);
130 chi2 =
fitTrack(seed_track,chi2_hits);
134 nhit_layer[hits[c_hit].layer]+=1;;
137 seed_track.
hits.pop_back();
144 unsigned int nt = input.size();
145 for(
unsigned int i=0;i<nt;++i)
147 output.push_back(input[i]);
160 chi2_hit.resize(track.
hits.size(), 0.);
162 MatrixXf
y = MatrixXf::Zero(track.
hits.size(), 1);
163 for(
unsigned int i=0;i<track.
hits.size();i++)
165 y(i, 0) = ( pow(track.
hits[i].x,2) + pow(track.
hits[i].y,2) );
170 MatrixXf
X = MatrixXf::Zero(track.
hits.size(), 3);
171 for(
unsigned int i=0;i<track.
hits.size();i++)
173 X(i, 0) = track.
hits[i].x;
174 X(i, 1) = track.
hits[i].y;
190 MatrixXf Xt = X.transpose();
192 MatrixXf prod = Xt*
X;
193 MatrixXf inv = prod.fullPivLu().inverse();
195 MatrixXf beta = inv*Xt*
y;
197 float cx = beta(0,0)*0.5;
198 float cy = beta(1,0)*0.5;
199 float r = sqrt(cx*cx + cy*cy - beta(2,0));
201 float phi = atan2(cy, cx);
202 float d = sqrt(cx*cx + cy*cy) -
r;
205 MatrixXf diff = y - (X*beta);
206 MatrixXf chi2 = (diff.transpose())*diff;
208 float dx = d*cos(phi);
209 float dy = d*sin(phi);
211 MatrixXf
y2 = MatrixXf::Zero(track.
hits.size(), 1);
212 for(
unsigned int i=0;i<track.
hits.size();i++)
214 y2(i,0) = track.
hits[i].z;
219 MatrixXf
X2 = MatrixXf::Zero(track.
hits.size(), 2);
220 for(
unsigned int i=0;i<track.
hits.size();i++)
222 float D = sqrt( pow(dx - track.
hits[i].x, 2) + pow(dy - track.
hits[i].y,2));
228 if(v >= 0.999999){v = 0.999999;}
259 MatrixXf Xt2 = X2.transpose();
260 MatrixXf prod2 = Xt2*
X2;
261 MatrixXf inv2 = prod2.fullPivLu().inverse();
262 MatrixXf beta2 = inv2*Xt2*
y2;
264 MatrixXf diff2 = y2 - (X2*beta2);
265 MatrixXf chi2_z = (diff2.transpose())*diff2;
267 float z0 = beta2(1,0);
268 float dzdl = beta2(0,0)/sqrt(1. + beta2(0,0)*beta2(0,0));
285 cx = (track.
d+
r)*cos(track.
phi);
286 cy = (track.
d+
r)*sin(track.
phi);
289 for(
unsigned int h=0;
h<track.
hits.size();
h++)
291 float dx1 = track.
hits[
h].x - cx;
292 float dy1 = track.
hits[
h].y - cy;
294 float dx2 = track.
hits[
h].x + cx;
295 float dy2 = track.
hits[
h].y + cy;
297 float xydiff1 = sqrt(dx1*dx1 + dy1*dy1) -
r;
298 float xydiff2 = sqrt(dx2*dx2 + dy2*dy2) -
r;
299 float xydiff = xydiff2;
300 if(fabs(xydiff1) < fabs(xydiff2)){ xydiff = xydiff1; }
309 chi2_hit[
h] += xydiff*xydiff/(ls_xy*ls_xy);
310 chi2_hit[
h] += diff2(
h,0)*diff2(
h,0);
312 chi2_tot += chi2_hit[
h];
315 unsigned int deg_of_freedom = 2*track.
hits.size() - 5;
319 track.
hits.pop_back();
323 return (chi2_tot)/((double)(deg_of_freedom));