ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CylindricalHough.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CylindricalHough.cpp
1 #include "CylindricalHough.h"
2 #include "ZHough_Cylindrical.h"
3 
4 
5 using namespace std;
6 
7 
8 bool CylindricalHough::intersect_circles(bool hel, double startx, double starty, double rad_det, double rad_trk, double cx, double cy, double& x, double& y)
9 {
10  double cx_det = -vertex_x;
11  double cy_det = -vertex_y;
12 
13  double d2 = ((cx-cx_det)*(cx-cx_det) + (cy-cy_det)*(cy-cy_det));
14  double d = sqrt(d2);
15  if(d > (rad_det + rad_trk))
16  {
17  return false;
18  }
19  if(d < fabs(rad_det - rad_trk))
20  {
21  return false;
22  }
23 
24  double r2 = rad_trk*rad_trk;
25 
26  double d_inv = 1./d;
27  double R2 = rad_det*rad_det;
28  double a = 0.5*(R2 - r2 + d2)*d_inv;
29  double h = a*d_inv;
30  double P2x = cx_det + (cx-cx_det)*h;
31  double P2y = cy_det + (cy-cy_det)*h;;
32  h = sqrt(R2 - a*a);
33 
34  double ux = -(cy-cy_det)*d_inv;
35  double uy = (cx-cx_det)*d_inv;
36  double P3x1 = P2x + ux*h;
37  double P3y1 = P2y + uy*h;
38  ux = -ux;
39  uy = -uy;
40  double P3x2 = P2x + ux*h;
41  double P3y2 = P2y + uy*h;
42 
43  double d1_2 = (startx - P3x1)*(startx - P3x1) + (starty - P3y1)*(starty - P3y1);
44  double d2_2 = (startx - P3x2)*(startx - P3x2) + (starty - P3y2)*(starty - P3y2);
45 
46  if(d1_2 < d2_2)
47  {
48  x = P3x1;
49  y = P3y1;
50  }
51  else
52  {
53  x = P3x2;
54  y = P3y2;
55  }
56 
57  return true;
58 }
59 
60 
61 CylindricalHough::CylindricalHough(vector<float>& detrad, unsigned int inv_radius_nbin, unsigned int center_angle_nbin, unsigned int dca_origin_nbin, CircleResolution& min_resolution, CircleResolution& max_resolution, CircleRange& range, unsigned int z0_nbin, unsigned int theta_nbin, ZResolution& minzres, ZResolution& maxzres, ZRange& zrange, double sxy, double sz) : CircleHough(inv_radius_nbin, center_angle_nbin, dca_origin_nbin, min_resolution, max_resolution, range, z0_nbin, theta_nbin, minzres, maxzres, zrange, sxy, sz, false), vertex_x(0.), vertex_y(0.), vertex_z(0.), phicut(0.1)
62 {
63  for(unsigned int i=0;i<detrad.size();++i)
64  {
65  detector_radii.push_back(detrad[i]);
66  }
67  init_ZHough(z0_nbin, theta_nbin, minzres, maxzres, zrange);
68 }
69 
70 
72 {
73 
74 }
75 
76 
77 void CylindricalHough::init_ZHough(int z0_nbin, unsigned int theta_nbin, ZResolution& minzres, ZResolution& maxzres, ZRange& zrange)
78 {
79  if(zhough!=NULL){delete zhough;}
80  zhough = new ZHough_Cylindrical(z0_nbin, theta_nbin, minzres, maxzres, zrange, sigma_xy, sigma_z);
81  ((ZHough_Cylindrical*)(zhough))->setNLayers(detector_radii.size());
82  vector<double> lxy, lz;
83  for(unsigned int i=0;i<detector_radii.size();++i)
84  {
85  lxy.push_back(0.005/3.);
86  lz.push_back(0.0425/3.);
87  }
88  ((ZHough_Cylindrical*)(zhough))->setLayerResolution(lxy, lz);
89  ((ZHough_Cylindrical*)(zhough))->setVertexResolution(0.002, 0.01);
90 
91  zhough->setCircleHough(this);
92 }
93 
94 
95 void CylindricalHough::setLayerResolution(vector<double>& lxy, vector<double>& lz)
96 {
97  ((ZHough_Cylindrical*)(zhough))->setLayerResolution(lxy, lz);
98 }
99 
100 
101 void CylindricalHough::setVertexResolution(double vxy, double vz)
102 {
103  ((ZHough_Cylindrical*)(zhough))->setVertexResolution(vxy, vz);
104 }
105 
106 
107 void CylindricalHough::customFindHelicesInit(vector<SimpleHit3D>& hits, unsigned int min_hits, unsigned int max_hits, unsigned int min_zhits, unsigned int max_zhits, double chi2_cut, float xydiffcut, vector<SimpleTrack3D>& tracks, unsigned int maxtracks)
108 {
109  AngleIndexList templist;
110  angle_list.clear();
111  angle_list.assign(detector_radii.size(), templist);
112  for(unsigned int i=0;i<hits.size();i++)
113  {
114  AngleIndexPair temppair(atan2(hits[i].get_y(), hits[i].get_x()),hits[i].index);
115  angle_list[hits[i].layer].addPair( temppair );
116  }
117 }
118 
119 
120 //params:
121 //0 <--> xydiffcut
122 //1 <--> chi2_cut
123 //2 <--> min_hits
124 //3 <--> max_kappa_cut
125 void CylindricalHough::addHits(unsigned int zlevel, vector<SimpleTrack3D>& temptracks, vector<SimpleTrack3D>& tracks, vector<float>& params, int tracks_per_hit, float z_cut)
126 {
127  vector<double> chi2_hit;
128 // float xydiffcut = params[0];
129  float chi2_cut = params[1];
130  unsigned int min_hits = (unsigned int)(params[2]);
131  float max_kappa_cut = params[3];
132 
133  for(unsigned int trk=0;trk<temptracks.size();++trk)
134  {
135  if(temptracks[trk].hits.size()<3 && using_vertex==false){continue;}
136  else if(temptracks[trk].hits.size()<2 && using_vertex==true){continue;}
137  zhough->fitTrack(temptracks[trk], chi2_hit);
138  SimpleTrack3D temptrack = temptracks[trk];
139 
140  float cx =0.0; // center of rotation, x
141  float cy =0.0; // center of rotation, y
142  float r =0.0; // radius of rotation
143 
144  if(fabs(temptracks[trk].kappa) < 1.0e-8){temptracks[trk].kappa = 1.0e-8;}
145  // radius or something very straight
146  r = 1./temptracks[trk].kappa;
147 
148  // center of rotation
149  cx = (temptracks[trk].d+r)*cos(temptracks[trk].phi);
150  cy = (temptracks[trk].d+r)*sin(temptracks[trk].phi);
151 
152  // counter-clockwise is represented by hel=true
153  double v1x = temptracks[trk].hits[0].get_x();
154  double v1y = temptracks[trk].hits[0].get_y();
155  double v1z = temptracks[trk].hits[0].get_z();
156  double v2x = temptracks[trk].hits.back().get_x();
157  double v2y = temptracks[trk].hits.back().get_y();
158  double diffx = v2x-v1x;
159  double diffy = v2y-v1y;
160  double cross = diffx*cy - diffy*cx;
161  bool helicity=true;
162  if(cross<0.){helicity=false;}
163 
164 
165  vector<bool> layer_used;
166  layer_used.assign(angle_list.size(), false);
167  for(unsigned int h=0;h<temptracks[trk].hits.size();++h)
168  {
169  layer_used[temptracks[trk].hits[h].layer]=true;
170  }
171  vector<int> unused_layers;
172  for(unsigned int ll=0;ll<layer_used.size();++ll)
173  {
174  if(layer_used[ll]==false){unused_layers.push_back(ll);}
175  }
176 
177  for(unsigned int ll=0;ll<unused_layers.size();ll++)
178  {
179  if(((unused_layers.size() - ll) + temptracks[trk].hits.size()) < min_hits){break;}
180 
181  int layer = unused_layers[ll];
182  int closest_layer = temptracks[trk].hits[0].layer;
183  double startx = v1x;
184  double starty = v1y;
185  double startz = v1z;
186  for(int cl=1;cl<(int)(temptracks[trk].hits.size());cl++)
187  {
188  if( fabs(layer - temptracks[trk].hits[cl].layer) < fabs(layer - closest_layer) )
189  {
190  closest_layer = temptracks[trk].hits[cl].layer;
191  startx = temptracks[trk].hits[cl].get_x();
192  starty = temptracks[trk].hits[cl].get_y();
193  startz = temptracks[trk].hits[cl].get_z();
194  }
195  }
196 
197  double x_intersect=0.;
198  double y_intersect=0.;
199  double phi_error = phicut;
200 
201  double xy_d_cut = 3.*((ZHough_Cylindrical*)(zhough))->getLayerResolution_xy(layer);
202  phi_error += xy_d_cut/(2.*M_PI*detector_radii[layer]);
203 // if(phi_error < 0.02){phi_error = 0.02;}
204 
205  bool try_helicity = helicity;
206 // if(layer < temptracks[trk].hits[0].layer){try_helicity = !helicity;}
207 
208  bool intersected = intersect_circles(try_helicity, startx, starty, detector_radii[layer], r, cx, cy, x_intersect, y_intersect);
209  if( ( intersected==false ) || ( x_intersect != x_intersect ) || ( y_intersect != y_intersect ) )
210  {
211  continue;
212  }
213  double intersect_phi = atan2(y_intersect, x_intersect);
214  if(intersect_phi < 0.){intersect_phi += 2.*M_PI;}
215  double k = temptracks[trk].kappa;
216  double D = sqrt((startx-x_intersect)*(startx-x_intersect) + (starty-y_intersect)*(starty-y_intersect));
217  double dzdl = temptracks[trk].dzdl;
218  double s=0.;
219  double z_intersect=0.;
220  if(0.5*k*D > 0.01)
221  {
222  double v = 0.5*k*D;
223  if(v >= 0.999999){v = 0.999999;}
224  s = 2.*asin(v)/k;
225  }
226  else
227  {
228  double temp1 = k*D*0.5;temp1*=temp1;
229  double temp2 = D*0.5;
230  s += 2.*temp2;
231  temp2*=temp1;
232  s += temp2/3.;
233  temp2*=temp1;
234  s += (3./20.)*temp2;
235  temp2*=temp1;
236  s += (5./56.)*temp2;
237  }
238  double dz = sqrt(s*s*dzdl*dzdl/(1. - dzdl*dzdl));
239  if(layer < closest_layer){dz=-dz;}
240  if(dzdl>0.){z_intersect = startz + dz;}
241  else{z_intersect = startz - dz;}
242  if(z_intersect != z_intersect)
243  {
244  continue;
245  }
246 
247  vector<AngleIndexPair*> hit_candidates;
248  angle_list[layer].getRangeList(intersect_phi, phi_error, hit_candidates);
249 
250  for(unsigned int hit_cand=0;hit_cand<hit_candidates.size();hit_cand++)
251  {
252 // //has hit been used in too many tracks?
253 // if(used_vec[hit_candidates[hit_cand]->index] >= tracks_per_hit){continue;}
254  //is this hit close enough in z?
255 
256  double z_error = z_cut + 3.*((ZHough_Cylindrical*)(zhough))->getLayerResolution_z(layer);
257  if(fabs( ((*(hits_vec[0]))[hit_candidates[hit_cand]->index]).get_z() - z_intersect ) >= z_error )
258  {
259  continue;
260  }
261 
262  temptracks[trk].hits.push_back((*(hits_vec[0]))[hit_candidates[hit_cand]->index]);
263  //fit the track with the hit candidate added on
264  double chi2 = zhough->fitTrack(temptracks[trk], chi2_hit);
265  //if this hit doesn't belong, then get rid of it
266  if( (chi2 != chi2) || (chi2 > chi2_cut) || (temptracks[trk].kappa > max_kappa_cut))
267  {
268  temptracks[trk] = temptrack;
269  }
270  else
271  {
272  temptrack = temptracks[trk];
273  break;
274  }
275  }
276  }
277  if(temptracks[trk].hits.size() >= min_hits)
278  {
279  // loop over all tracks
280  for(unsigned int itrack=0; itrack<tracks.size();itrack++)
281  {
282  // sort all hit indexes within the track
283  sort (tracks[itrack].hits.begin(), tracks[itrack].hits.end(),SimpleHit3D_LessThan);
284  }
285 
286  // now sort the tracks themselves
287  sort (tracks.begin(), tracks.end(), SimpleTrack3D_LessThan);
288 
289  sort (temptracks[trk].hits.begin(), temptracks[trk].hits.end(),SimpleHit3D_LessThan);
290 
291 
292  if(binary_search(tracks.begin(), tracks.end(), temptracks[trk], SimpleTrack3D_LessThan) == false)
293  {
294  // add to the output list
295  tracks.push_back(temptracks[trk]);
296  // add the hits to the global usage list
297  for(unsigned int h=0;h<temptracks[trk].hits.size();h++)
298  {
299  // doesn't this double count the already used hits??? -MPM
300  // no, actually this is the first time the hits get marked -Theo K
301  // so the usage in ZHough during the track construction is not passed up??? -MPM
302  used_vec[temptracks[trk].hits[h].index]++;
303  }
304  }
305  }
306 
307  }
308 
309 // //---------------------------
310 // // Eliminate duplicate tracks
311 // //---------------------------
312 //
313 // // loop over all tracks
314 // for(unsigned int itrack=0; itrack<tracks.size();itrack++)
315 // {
316 // // sort all hit indexes within the track
317 // sort (tracks[itrack].hits.begin(), tracks[itrack].hits.end(),SimpleHit3D_LessThan);
318 // }
319 //
320 // // now sort the tracks themselves
321 // sort (tracks.begin(), tracks.end(), SimpleTrack3D_LessThan);
322 //
323 // // remove the duplicate tracks which will be neighbors after the sorting above
324 // vector<SimpleTrack3D>::iterator it = unique (tracks.begin(), tracks.end(), SimpleTrack3D_Equality);
325 //
326 // // resize the vector to eliminate the duplicates
327 // tracks.resize( it - tracks.begin() );
328 
329 }
330 
331