ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ThreeHitSeedGrower.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ThreeHitSeedGrower.cpp
1 #include "ThreeHitSeedGrower.h"
2 #include <cmath>
3 #include <iostream>
4 #include <algorithm>
5 #include <Eigen/LU>
6 #include <Eigen/Core>
7 
8 
9 using namespace std;
10 using namespace Eigen;
11 
12 
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.)
14 {
15  unsigned int n_layers = detrad.size();
16  _max_hits = n_layers + 2;
17  for(unsigned int i=0;i<n_layers;++i)
18  {
19  detector_radii.push_back(detrad[i]);
20  }
21 }
22 
23 
24 void ThreeHitSeedGrower::findTracks(vector<SimpleHit3D>& hits, vector<SimpleTrack3D>& tracks)
25 {cout<<"findTracks::entered with "<<hits.size()<<" hits"<<endl;
26  vector<double> chi2_hits;
27  double chi2;
28  SimpleTrack3D temp_track;
29  temp_track.hits.resize(3, hits[0]);
30  for(unsigned int i1=0;i1<hits.size();++i1)
31  {
32  temp_track.hits[0] = hits[i1];
33  for(unsigned int i2=(i1+1);i2<hits.size();++i2)
34  {//TO DO: ADD PHICUT
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)
38  {
39  if((hits[i3].layer == hits[i2].layer) || (hits[i3].layer == hits[i1].layer)){continue;}
40  temp_track.hits[2] = hits[i3];
41 
42  chi2 = fitTrack(temp_track, chi2_hits);
43  if(chi2 > chi2_cut){continue;}
44  vector<unsigned int> nhit_layer;
45  nhit_layer.assign(detector_radii.size(),0);
46  for(unsigned int ihit = 0; ihit<3;ihit++)
47  {
48  nhit_layer[temp_track.hits[ihit].layer]+=1;
49  }
50  if(GrowTrack(temp_track, nhit_layer, hits, tracks, 0) == false)
51  {//refit if track wasn't grown.
52  vector<unsigned int> tempcomb;
53  tempcomb.assign(3,0);
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;}
60  combos.insert(tempcomb);
61  chi2 = fitTrack(temp_track, chi2_hits);
62  tracks.push_back(temp_track);
63  cout<<"findTrack::added a 3 hit track"<<endl;
64  }
65 
66  }
67  }
68  }
69  cout<<"leaving findTrack"<<endl;
70 }
71 //this is my attempt at implementation of 3 hit seeding adding hits up to 8 total
72 //plan for this function:
73 //will figure out which layers need new hits the most from nhit_layer
74 //incrementally adds a hit from hits and fits it
75 //tests to see if the chi squared is still in bound
76 //recursively calls grow track on newly accepted track and adds it to tracks.
77 //
78 bool ThreeHitSeedGrower::GrowTrack(SimpleTrack3D seed_track, vector<unsigned int>& nhit_layer, vector<SimpleHit3D>& hits, vector<SimpleTrack3D>& tracks, unsigned int c_hit){
79  unsigned int init_size = seed_track.hits.size();
80  vector<double> chi2_hits;
81  double chi2;
82  bool hit_added;
83  while(c_hit < hits.size())
84  {
85  hit_added = addOneHit(seed_track, nhit_layer, c_hit, hits);
86  c_hit+=1;
87  if(hit_added == false)
88  continue;
89  unsigned int final_size = seed_track.hits.size();
90  if(final_size < _max_hits)
91  {
92  if(GrowTrack(seed_track, nhit_layer, hits, tracks, c_hit) == false)
93  //attempt to continue growing track. if fails, reset fitting parameters
94  {
95  chi2 = fitTrack(seed_track,chi2_hits);
96  }
97  }
98  vector<unsigned int> tempcomb;
99  tempcomb.assign(seed_track.hits.size(),0);
100  for(unsigned int ihit = 0; ihit < seed_track.hits.size();ihit++)
101  tempcomb[ihit] = seed_track.hits[ihit].index;
102  sort(tempcomb.begin(), tempcomb.end());
103  set<vector<unsigned int> >::iterator it = combos.find(tempcomb);
104  if(it != combos.end()){
105  seed_track.hits.pop_back();
106  continue;
107  }
108  combos.insert(tempcomb);
109  tracks.push_back(seed_track);
110  cout<<"GrowTrack::added a track of size: "<<seed_track.hits.size()<<endl;
111  return true;
112  }
113  return false;
114 }
115 
116 bool ThreeHitSeedGrower::addOneHit(SimpleTrack3D & seed_track, vector<unsigned int> & nhit_layer, unsigned int c_hit, vector<SimpleHit3D>& hits)
117 {//cout<<"addOneHit::newly entered"<<endl;
118  for(unsigned int ihit =0; ihit<seed_track.hits.size();ihit++)
119  {
120  // cout<<"addOneHit::hit: "<<ihit<<" index = "<<seed_track.hits[ihit].index<<endl;
121  if(seed_track.hits[ihit].index == hits[c_hit].index)
122  return false;
123  }
124  // cout<<"addOneHit::hit: "<<c_hit<<" passed with index = "<<hits[c_hit].index<<endl;
125  //this conditional logic specifically applies to the VTX, maybe modularize this later.
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;
128  double chi2;
129  seed_track.hits.push_back(hits[c_hit]);
130  chi2 = fitTrack(seed_track,chi2_hits);
131  if(chi2 <= chi2_cut)
132  {
133  return true;
134  nhit_layer[hits[c_hit].layer]+=1;;
135  }
136  else
137  seed_track.hits.pop_back();
138  }
139  return false;
140 }
141 
142 void ThreeHitSeedGrower::finalize(vector<SimpleTrack3D>& input, vector<SimpleTrack3D>& output)
143 {
144  unsigned int nt = input.size();
145  for(unsigned int i=0;i<nt;++i)
146  {
147  output.push_back(input[i]);
148  }
149 }
150 
151 
152 double ThreeHitSeedGrower::fitTrack(SimpleTrack3D & track, vector<double>& chi2_hit)
153 {
154  if(using_vertex == true)
155  {
156  track.hits.push_back(SimpleHit3D(0.,0., 0.,0., 0.,0., 0, 0));
157  }
158 
159  chi2_hit.clear();
160  chi2_hit.resize(track.hits.size(), 0.);
161 
162  MatrixXf y = MatrixXf::Zero(track.hits.size(), 1);
163  for(unsigned int i=0;i<track.hits.size();i++)
164  {
165  y(i, 0) = ( pow(track.hits[i].x,2) + pow(track.hits[i].y,2) );
166  if((using_vertex==true ) && (i == (track.hits.size() - 1))){y(i, 0) /= vertex_sigma_xy;}
167  else{y(i, 0) /= layer_xy_resolution[track.hits[i].layer];}
168  }
169 
170  MatrixXf X = MatrixXf::Zero(track.hits.size(), 3);
171  for(unsigned int i=0;i<track.hits.size();i++)
172  {
173  X(i, 0) = track.hits[i].x;
174  X(i, 1) = track.hits[i].y;
175  X(i, 2) = -1.;
176  if((using_vertex==true ) && (i == (track.hits.size() - 1)))
177  {
178  X(i, 0) /= vertex_sigma_xy;
179  X(i, 1) /= vertex_sigma_xy;
180  X(i, 2) /= vertex_sigma_xy;
181  }
182  else
183  {
184  X(i, 0) /= layer_xy_resolution[track.hits[i].layer];
185  X(i, 1) /= layer_xy_resolution[track.hits[i].layer];
186  X(i, 2) /= layer_xy_resolution[track.hits[i].layer];
187  }
188  }
189 
190  MatrixXf Xt = X.transpose();
191 
192  MatrixXf prod = Xt*X;
193  MatrixXf inv = prod.fullPivLu().inverse();
194 
195  MatrixXf beta = inv*Xt*y;
196 
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));
200 
201  float phi = atan2(cy, cx);
202  float d = sqrt(cx*cx + cy*cy) - r;
203  float k = 1./r;
204 
205  MatrixXf diff = y - (X*beta);
206  MatrixXf chi2 = (diff.transpose())*diff;
207 
208  float dx = d*cos(phi);
209  float dy = d*sin(phi);
210 
211  MatrixXf y2 = MatrixXf::Zero(track.hits.size(), 1);
212  for(unsigned int i=0;i<track.hits.size();i++)
213  {
214  y2(i,0) = track.hits[i].z;
215  if((using_vertex==true ) && (i == (track.hits.size() - 1))){y2(i, 0) /= vertex_sigma_z;}
216  else{y2(i, 0) /= layer_z_resolution[track.hits[i].layer];}
217  }
218 
219  MatrixXf X2 = MatrixXf::Zero(track.hits.size(), 2);
220  for(unsigned int i=0;i<track.hits.size();i++)
221  {
222  float D = sqrt( pow(dx - track.hits[i].x, 2) + pow(dy - track.hits[i].y,2));
223  float s = 0.0;
224 
225  if(0.5*k*D > 0.1)
226  {
227  float v = 0.5*k*D;
228  if(v >= 0.999999){v = 0.999999;}
229  s = 2.*asin(v)/k;
230  }
231  else
232  {
233  float temp1 = k*D*0.5;temp1*=temp1;
234  float temp2 = D*0.5;
235  s += 2.*temp2;
236  temp2*=temp1;
237  s += temp2/3.;
238  temp2*=temp1;
239  s += (3./20.)*temp2;
240  temp2*=temp1;
241  s += (5./56.)*temp2;
242  }
243 
244  X2(i,0) = s;
245  X2(i,1) = 1.0;
246 
247  if((using_vertex==true ) && (i == (track.hits.size() - 1)))
248  {
249  X2(i, 0) /= vertex_sigma_z;
250  X2(i, 1) /= vertex_sigma_z;
251  }
252  else
253  {
254  X2(i, 0) /= layer_z_resolution[track.hits[i].layer];
255  X2(i, 1) /= layer_z_resolution[track.hits[i].layer];
256  }
257  }
258 
259  MatrixXf Xt2 = X2.transpose();
260  MatrixXf prod2 = Xt2*X2;
261  MatrixXf inv2 = prod2.fullPivLu().inverse();
262  MatrixXf beta2 = inv2*Xt2*y2;
263 
264  MatrixXf diff2 = y2 - (X2*beta2);
265  MatrixXf chi2_z = (diff2.transpose())*diff2;
266 
267  float z0 = beta2(1,0);
268  float dzdl = beta2(0,0)/sqrt(1. + beta2(0,0)*beta2(0,0));
269 
270  track.phi = phi;
271  track.d = d;
272  track.kappa = k;
273  track.dzdl = dzdl;
274  track.z0 = z0;
275 
276  if(track.kappa!=0.)
277  {
278  r=1./track.kappa;
279  }
280  else
281  {
282  r=1.0e10;
283  }
284 
285  cx = (track.d+r)*cos(track.phi);
286  cy = (track.d+r)*sin(track.phi);
287 
288  float chi2_tot = 0.;
289  for(unsigned int h=0;h<track.hits.size();h++)
290  {
291  float dx1 = track.hits[h].x - cx;
292  float dy1 = track.hits[h].y - cy;
293 
294  float dx2 = track.hits[h].x + cx;
295  float dy2 = track.hits[h].y + cy;
296 
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; }
301 
302  float ls_xy = layer_xy_resolution[track.hits[h].layer];
303  if((using_vertex == true) && (h == (track.hits.size() - 1)))
304  {
305  ls_xy = vertex_sigma_xy;
306  }
307 
308  chi2_hit[h] = 0.;
309  chi2_hit[h] += xydiff*xydiff/(ls_xy*ls_xy);
310  chi2_hit[h] += diff2(h,0)*diff2(h,0);
311 
312  chi2_tot += chi2_hit[h];
313  }
314 
315  unsigned int deg_of_freedom = 2*track.hits.size() - 5;
316 
317  if(using_vertex == true)
318  {
319  track.hits.pop_back();
320  chi2_hit.pop_back();
321  }
322 
323  return (chi2_tot)/((double)(deg_of_freedom));
324 }
325 
326