ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ZHough_Cylindrical.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ZHough_Cylindrical.cpp
1 #include "ZHough_Cylindrical.h"
2 #include <math.h>
3 #include <iostream>
4 #include <algorithm>
5 #include "CircleHough.h"
6 #include <Eigen/LU>
7 #include <Eigen/Core>
8 
9 
10 using namespace std;
11 using namespace Eigen;
12 
13 
14 ZHough_Cylindrical::ZHough_Cylindrical(unsigned int z0_nbin, unsigned int theta_nbin, ZResolution& min_resolution, ZResolution& max_resolution, ZRange& range, double sxy, double sz) : ZHough(z0_nbin, theta_nbin, min_resolution, max_resolution, range, sxy, sz), nlayers(6)
15 {
16 
17 }
18 
19 
21 {
22 
23 }
24 
25 
26 void ZHough_Cylindrical::findTracks(unsigned int zoomlevel, float xydiffcut, unsigned int max_hits, unsigned int tracks_per_hit, double chi2_cut, float max_kappa_cut, vector<SimpleTrack3D>& tracks, vector<SimpleHit3D>& hits, float min_k, float max_k, float min_phi, float max_phi, float min_d, float max_d, float min_z0, float max_z0, float min_theta, float max_theta, vector<float>& params)
27 {
28  if(_verbosity > 0)
29  {
30  cout<<"findTracks"<<endl;
31  cout<<min_k<<" "<<max_k<<" "<<min_phi<<" "<<max_phi<<" "<<min_d<<" "<<max_d<<" "<<min_z0<<" "<<max_z0<<" "<<min_theta<<" "<<max_theta<<endl;
32  }
33 
34  vector<SimpleHit3D> comb_hits = (*(hits_vec[zoomlevel+1]));
35  if(using_vertex==false)
36  {
37  findTracksCombo_noVertex(comb_hits, zoomlevel, xydiffcut, max_hits, tracks_per_hit, chi2_cut, max_kappa_cut, tracks, hits, min_k, max_k, min_phi, max_phi, min_d, max_d, min_z0, max_z0, min_theta, max_theta, params);
38  }
39  else
40  {
41  findTracksCombo_withVertex(comb_hits, zoomlevel, xydiffcut, max_hits, tracks_per_hit, chi2_cut, max_kappa_cut, tracks, hits, min_k, max_k, min_phi, max_phi, min_d, max_d, min_z0, max_z0, min_theta, max_theta, params);
42  }
43 
44 }
45 
46 
47 void ZHough_Cylindrical::findTracksCombo_noVertex(vector<SimpleHit3D>& comb_hits, unsigned int zoomlevel, float xydiffcut, unsigned int max_hits, unsigned int tracks_per_hit, double chi2_cut, float max_kappa_cut, vector<SimpleTrack3D>& tracks, vector<SimpleHit3D>& hits, float min_k, float max_k, float min_phi, float max_phi, float min_d, float max_d, float min_z0, float max_z0, float min_theta, float max_theta, vector<float>& params)
48 {
49  unsigned int req_hits = (unsigned int)(params[2]);
50 
51  vector<unsigned int> index_holder;
52  index_holder.assign(comb_hits.size(), 0);
53  for(unsigned int i=0;i<comb_hits.size();++i)
54  {
55  index_holder[i]=comb_hits[i].index;
56  comb_hits[i].index = i;
57  }
58 
59  double phidiff;
60  vector<double> chi2_hits;
61  double chi2=0.;
62  vector<bool> use_here;
63  use_here.assign(comb_hits.size(), true);
64  SimpleTrack3D temp_track;
65  temp_track.hits.resize(3, comb_hits[0]);
66 
67  //loop over combinations of 3 hits, each from a different layer
68  for(unsigned int i1=0;i1<comb_hits.size();++i1)
69  {
70  if(use_here[i1]==false){continue;}
71  if(used_vec[index_holder[i1]] >= tracks_per_hit){continue;}
72  temp_track.hits[0] = comb_hits[i1];
73  for(unsigned int i2=(i1+1);i2<comb_hits.size();++i2)
74  {
75  if(use_here[i2]==false){continue;}
76  if(used_vec[index_holder[i1]] >= tracks_per_hit){continue;}
77  if(used_vec[index_holder[i2]] >= tracks_per_hit){continue;}
78  if(comb_hits[i2].layer == comb_hits[i1].layer){continue;}
79  phidiff = temp_track.hits[0].get_phi() - comb_hits[i2].get_phi();
80  if(phidiff < -M_PI){phidiff += 2.*M_PI;}
81  if(phidiff > M_PI){phidiff -= 2.*M_PI;}
82  if(fabs(phidiff) > M_PI/3.){continue;}
83  temp_track.hits[1] = comb_hits[i2];
84  float xdiff2 = (temp_track.hits[0].get_x() - temp_track.hits[1].get_x());
85  xdiff2*=xdiff2;
86  float ydiff2 = (temp_track.hits[0].get_y() - temp_track.hits[1].get_y());
87  ydiff2*=ydiff2;
88  float xydiff2 = xdiff2 + ydiff2;
89  float zdiff2 = (temp_track.hits[0].get_z() - temp_track.hits[1].get_z());
90  zdiff2*=zdiff2;
91  for(unsigned int i3=(i2+1);i3<comb_hits.size();++i3)
92  {
93  if(use_here[i3]==false){continue;}
94  if(used_vec[index_holder[i1]] >= tracks_per_hit){continue;}
95  if(used_vec[index_holder[i2]] >= tracks_per_hit){continue;}
96  if(used_vec[index_holder[i3]] >= tracks_per_hit){continue;}
97  if((comb_hits[i3].layer == comb_hits[i2].layer) || (comb_hits[i3].layer == comb_hits[i1].layer)){continue;}
98  phidiff = temp_track.hits[1].get_phi() - comb_hits[i3].get_phi();
99  if(phidiff < -M_PI){phidiff += 2.*M_PI;}
100  if(phidiff > M_PI){phidiff -= 2.*M_PI;}
101  if(fabs(phidiff) > M_PI/3.){continue;}
102  temp_track.hits[2] = comb_hits[i3];
103  float xdiff2_2 = (temp_track.hits[1].get_x() - temp_track.hits[2].get_x());
104  xdiff2_2*=xdiff2_2;
105  float ydiff2_2 = (temp_track.hits[1].get_y() - temp_track.hits[2].get_y());
106  ydiff2_2*=ydiff2_2;
107  float xydiff2_2 = xdiff2_2 + ydiff2_2;
108  float zdiff2_2 = (temp_track.hits[1].get_z() - temp_track.hits[2].get_z());
109  zdiff2_2*=zdiff2_2;
110 
111  if(fabs(xydiff2*zdiff2_2/(zdiff2*xydiff2_2) - 1.) > 0.2){continue;}
112 
113  temp_track.hits.resize(3, comb_hits[0]);//make sure there are only 3 hits in this track
114  chi2 = fitTrack(temp_track, chi2_hits);
115  if((chi2 < chi2_cut)&&(temp_track.kappa < max_kappa_cut))
116  {
117  //we found a track candidate. see if there are more hits in comb_hits to add to this track
118  vector<SimpleTrack3D> ttracks;
119  SimpleTrack3D ttrack = temp_track;
120  for(unsigned int j=0;j<temp_track.hits.size();++j)
121  {
122  ttrack.hits[j].index = index_holder[ttrack.hits[j].index];
123  ttrack.hits[j].index = index_mapping[ttrack.hits[j].index];
124  }
125  ttracks.push_back(ttrack);
126  unsigned int start_size = tracks.size();
127  chough->addHits((unsigned int)(params[4]), ttracks, tracks, params, tracks_per_hit, params[5]);
128  if(tracks.size() > start_size)
129  {
130  if(tracks.back().hits.size() >= req_hits)
131  {
132  unsigned int count = 0;
133  for(unsigned int j=0;j<temp_track.hits.size();++j)
134  {
135  use_here[temp_track.hits[j].index] = false;
136  }
137  for(unsigned int j=0;j<temp_track.hits.size();++j)
138  {
139  used_vec[index_holder[temp_track.hits[j].index]]++;
140  count++;
141  }
142  for(unsigned int j=temp_track.hits.size();j<tracks.back().hits.size();++j)
143  {
144  for(unsigned int i=0;i<(*(hits_vec[zoomlevel+1])).size();++i)
145  {
146  if(tracks.back().hits[j].index == index_mapping[(*(hits_vec[zoomlevel+1]))[i].index])
147  {
148  used_vec[(*(hits_vec[zoomlevel+1]))[i].index]++;
149  count++;
150  }
151  }
152  }
153  }
154  }
155  }
156  }
157  }
158  }
159  //restore indexes
160  for(unsigned int i=0;i<comb_hits.size();++i)
161  {
162  comb_hits[i].index = index_holder[i];
163  }
164 }
165 
166 
167 void ZHough_Cylindrical::findTracksCombo_withVertex(vector<SimpleHit3D>& comb_hits, unsigned int zoomlevel, float xydiffcut, unsigned int max_hits, unsigned int tracks_per_hit, double chi2_cut, float max_kappa_cut, vector<SimpleTrack3D>& tracks, vector<SimpleHit3D>& hits, float min_k, float max_k, float min_phi, float max_phi, float min_d, float max_d, float min_z0, float max_z0, float min_theta, float max_theta, vector<float>& params)
168 {
169  unsigned int req_hits = (unsigned int)(params[2]);
170 
171  vector<unsigned int> index_holder;
172  index_holder.assign(comb_hits.size(), 0);
173  for(unsigned int i=0;i<comb_hits.size();++i)
174  {
175  index_holder[i]=comb_hits[i].index;
176  comb_hits[i].index = i;
177  }
178 
179  double phidiff;
180  vector<double> chi2_hits;
181  double chi2=0.;
182  vector<bool> use_here;
183  use_here.assign(comb_hits.size(), true);
184  SimpleTrack3D temp_track;
185  temp_track.hits.resize(2, comb_hits[0]);
186 
187  //loop over combinations of 2 hits, each from a different layer
188  for(unsigned int i1=0;i1<comb_hits.size();++i1)
189  {
190  if(use_here[i1]==false){continue;}
191  if(used_vec[index_holder[i1]] >= tracks_per_hit){continue;}
192  temp_track.hits[0] = comb_hits[i1];
193 
194  float xdiff2_v = (temp_track.hits[0].get_x());
195  xdiff2_v*=xdiff2_v;
196  float ydiff2_v = (temp_track.hits[0].get_y());
197  ydiff2_v*=ydiff2_v;
198  float xydiff2_v = xdiff2_v + ydiff2_v;
199  float zdiff2_v = (temp_track.hits[0].get_z());
200  zdiff2_v*=zdiff2_v;
201 
202  for(unsigned int i2=(i1+1);i2<comb_hits.size();++i2)
203  {
204  if(use_here[i1]==false){continue;}
205  if(use_here[i2]==false){continue;}
206  if(comb_hits[i2].layer == comb_hits[i1].layer){continue;}
207  if(used_vec[index_holder[i1]] >= tracks_per_hit){continue;}
208  if(used_vec[index_holder[i2]] >= tracks_per_hit){continue;}
209  phidiff = temp_track.hits[0].get_phi() - comb_hits[i2].get_phi();
210  if(phidiff < -M_PI){phidiff += 2.*M_PI;}
211  if(phidiff > M_PI){phidiff -= 2.*M_PI;}
212  if(fabs(phidiff) > M_PI/3.){continue;}
213  temp_track.hits[1] = comb_hits[i2];
214  float xdiff2 = (temp_track.hits[0].get_x() - temp_track.hits[1].get_x());
215  xdiff2*=xdiff2;
216  float ydiff2 = (temp_track.hits[0].get_y() - temp_track.hits[1].get_y());
217  ydiff2*=ydiff2;
218  float xydiff2 = xdiff2 + ydiff2;
219  float zdiff2 = (temp_track.hits[0].get_z() - temp_track.hits[1].get_z());
220  zdiff2*=zdiff2;
221  if(fabs(xydiff2*zdiff2_v/(zdiff2*xydiff2_v) - 1.) > 0.2){continue;}
222 
223  temp_track.hits.resize(2, comb_hits[0]);//make sure there are only 2 hits in this track
224  chi2 = fitTrack(temp_track, chi2_hits);
225  if((chi2 < chi2_cut)&&(temp_track.kappa < max_kappa_cut))
226  {
227  //we found a track candidate. see if there are more hits in comb_hits to add to this track
228  vector<SimpleTrack3D> ttracks;
229  SimpleTrack3D ttrack = temp_track;
230  for(unsigned int j=0;j<temp_track.hits.size();++j)
231  {
232  ttrack.hits[j].index = index_holder[ttrack.hits[j].index];
233  ttrack.hits[j].index = index_mapping[ttrack.hits[j].index];
234  }
235  ttracks.push_back(ttrack);
236  unsigned int start_size = tracks.size();
237  chough->addHits((unsigned int)(params[4]), ttracks, tracks, params, tracks_per_hit, params[5]);
238  if(tracks.size() > start_size)
239  {
240  if(tracks.back().hits.size() >= req_hits)
241  {
242  unsigned int count = 0;
243  for(unsigned int j=0;j<temp_track.hits.size();++j)
244  {
245  use_here[temp_track.hits[j].index] = false;
246  }
247  for(unsigned int j=0;j<temp_track.hits.size();++j)
248  {
249  used_vec[index_holder[temp_track.hits[j].index]]++;
250  count++;
251  }
252  for(unsigned int j=temp_track.hits.size();j<tracks.back().hits.size();++j)
253  {
254  for(unsigned int i=0;i<(*(hits_vec[zoomlevel+1])).size();++i)
255  {
256  if(tracks.back().hits[j].index == index_mapping[(*(hits_vec[zoomlevel+1]))[i].index])
257  {
258  used_vec[(*(hits_vec[zoomlevel+1]))[i].index]++;
259  count++;
260  }
261  }
262  }
263  }
264  }
265  }
266  }
267  }
268 }
269 
270 
271 
272 double ZHough_Cylindrical::fitTrack(SimpleTrack3D& track, vector<double>& chi2_hit)
273 {
274  if(using_vertex == true)
275  {
276  track.hits.push_back(SimpleHit3D(0.,0.,0., 0, 0));
277  }
278 
279  chi2_hit.clear();
280  chi2_hit.resize(track.hits.size(), 0.);
281 
282  MatrixXd y = MatrixXd::Zero(track.hits.size(), 1);
283  for(unsigned int i=0;i<track.hits.size();i++)
284  {
285  y(i, 0) = ( pow(track.hits[i].get_x(),2) + pow(track.hits[i].get_y(),2) );
286  if((using_vertex==true ) && (i == (track.hits.size() - 1))){y(i, 0) /= vertex_sigma_xy;}
287  else{y(i, 0) /= layer_xy_resolution[track.hits[i].layer];}
288  }
289 
290  MatrixXd X = MatrixXd::Zero(track.hits.size(), 3);
291  for(unsigned int i=0;i<track.hits.size();i++)
292  {
293  X(i, 0) = track.hits[i].get_x();
294  X(i, 1) = track.hits[i].get_y();
295  X(i, 2) = -1.;
296  if((using_vertex==true ) && (i == (track.hits.size() - 1)))
297  {
298  X(i, 0) /= vertex_sigma_xy;
299  X(i, 1) /= vertex_sigma_xy;
300  X(i, 2) /= vertex_sigma_xy;
301  }
302  else
303  {
304  X(i, 0) /= layer_xy_resolution[track.hits[i].layer];
305  X(i, 1) /= layer_xy_resolution[track.hits[i].layer];
306  X(i, 2) /= layer_xy_resolution[track.hits[i].layer];
307  }
308  }
309 
310  MatrixXd Xt = X.transpose();
311 
312  MatrixXd prod = Xt*X;
313  MatrixXd inv = prod.fullPivLu().inverse();
314 
315  MatrixXd beta = inv*Xt*y;
316 
317  double cx = beta(0,0)*0.5;
318  double cy = beta(1,0)*0.5;
319  double r = sqrt(cx*cx + cy*cy - beta(2,0));
320 
321  double phi = atan2(cy, cx);
322  double d = sqrt(cx*cx + cy*cy) - r;
323  double k = 1./r;
324 
325  MatrixXd diff = y - (X*beta);
326  MatrixXd chi2 = (diff.transpose())*diff;
327  for(unsigned int i=0;i<track.hits.size();i++)
328  {
329  chi2_hit[i] += diff(i,0)*diff(i,0);
330  }
331 
332  double dx = d*cos(phi);
333  double dy = d*sin(phi);
334 
335  MatrixXd y2 = MatrixXd::Zero(track.hits.size(), 1);
336  for(unsigned int i=0;i<track.hits.size();i++)
337  {
338  y2(i,0) = track.hits[i].get_z();
339  if((using_vertex==true ) && (i == (track.hits.size() - 1))){y2(i, 0) /= vertex_sigma_z;}
340  else{y2(i, 0) /= layer_z_resolution[track.hits[i].layer];}
341  }
342 
343  MatrixXd X2 = MatrixXd::Zero(track.hits.size(), 2);
344  for(unsigned int i=0;i<track.hits.size();i++)
345  {
346  double D = sqrt( pow(dx - track.hits[i].get_x(), 2) + pow(dy - track.hits[i].get_y(),2));
347  double s = 0.0;
348 
349  if(0.5*k*D > 0.1)
350  {
351  double v = 0.5*k*D;
352  if(v >= 0.999999){v = 0.999999;}
353  s = 2.*asin(v)/k;
354  }
355  else
356  {
357  double temp1 = k*D*0.5;temp1*=temp1;
358  double temp2 = D*0.5;
359  s += 2.*temp2;
360  temp2*=temp1;
361  s += temp2/3.;
362  temp2*=temp1;
363  s += (3./20.)*temp2;
364  temp2*=temp1;
365  s += (5./56.)*temp2;
366  }
367 
368  X2(i,0) = s;
369  X2(i,1) = 1.0;
370 
371  if((using_vertex==true ) && (i == (track.hits.size() - 1)))
372  {
373  X2(i, 0) /= vertex_sigma_z;
374  X2(i, 1) /= vertex_sigma_z;
375  }
376  else
377  {
378  X2(i, 0) /= layer_z_resolution[track.hits[i].layer];
379  X2(i, 1) /= layer_z_resolution[track.hits[i].layer];
380  }
381  }
382 
383  MatrixXd Xt2 = X2.transpose();
384  MatrixXd prod2 = Xt2*X2;
385  MatrixXd inv2 = prod2.fullPivLu().inverse();
386  MatrixXd beta2 = inv2*Xt2*y2;
387 
388  MatrixXd diff2 = y2 - (X2*beta2);
389  MatrixXd chi2_z = (diff2.transpose())*diff2;
390  for(unsigned int i=0;i<track.hits.size();i++)
391  {
392  chi2_hit[i] += diff2(i,0)*diff2(i,0);
393  }
394 
395  double z0 = beta2(1,0);
396  double dzdl = beta2(0,0)/sqrt(1. + beta2(0,0)*beta2(0,0));
397 
398  track.phi = phi;
399  track.d = d;
400  track.kappa = k;
401  track.dzdl = dzdl;
402  track.z0 = z0;
403 // cout<<"fit params = "<<k<<" "<<phi<<" "<<d<<" "<<" "<<z0<<" "<<dzdl<<endl;
404 
405  if(track.kappa!=0.)
406  {
407  r=1./track.kappa;
408  }
409  else
410  {
411  r=1.0e10;
412  }
413 
414  cx = (track.d+r)*cos(track.phi);
415  cy = (track.d+r)*sin(track.phi);
416 
417  double chi2_tot = 0.;
418  for(unsigned int h=0;h<track.hits.size();h++)
419  {
420  double dx1 = track.hits[h].get_x() - cx;
421  double dy1 = track.hits[h].get_y() - cy;
422 
423  double dx2 = track.hits[h].get_x() + cx;
424  double dy2 = track.hits[h].get_y() + cy;
425 
426  double xydiff1 = sqrt(dx1*dx1 + dy1*dy1) - r;
427  double xydiff2 = sqrt(dx2*dx2 + dy2*dy2) - r;
428  double xydiff = xydiff2;
429  if(fabs(xydiff1) < fabs(xydiff2)){ xydiff = xydiff1; }
430 
431  double ls_xy = layer_xy_resolution[track.hits[h].layer];
432  double ls_z = layer_z_resolution[track.hits[h].layer];
433  if((using_vertex == true) && (h == (track.hits.size() - 1)))
434  {
435  ls_xy = vertex_sigma_xy;
436  ls_z = vertex_sigma_z;
437  }
438 
439  chi2_hit[h] = 0.;
440  chi2_hit[h] += xydiff*xydiff/(ls_xy*ls_xy);
441  chi2_hit[h] += diff2(h,0)*diff2(h,0);
442 // chi2_hit[h] += diff2(h,0)*diff2(h,0)/(ls_z*ls_z);
443 
444  chi2_tot += chi2_hit[h];
445  }
446 
447  unsigned int deg_of_freedom = 2*track.hits.size() - 5;
448 
449  if(using_vertex == true)
450  {
451  track.hits.pop_back();
452  chi2_hit.pop_back();
453  }
454 
455 // for(unsigned int i=0;i<track.hits.size();++i)
456 // {
457 // cout<<"chi^2 ( "<<track.hits[i].layer<<" ) = "<<chi2_hit[i]<<endl;
458 // }
459 // cout<<endl;
460 
461  return (chi2_tot)/((double)(deg_of_freedom));
462 }
463