EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HelixHoughFuncs_v1.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HelixHoughFuncs_v1.cc
1 #include "HelixHoughFuncs_v1.h"
2 
3 #include <cmath>
4 #include <vector>
5 
6 using namespace std;
7 
8 HelixHoughFuncs_v1::HelixHoughFuncs_v1()
9  : _hough_space(nullptr),_cur_zoom(0)
10 {
11 
12 }
13 
14 HelixHoughFuncs_v1::HelixHoughFuncs_v1(const HelixHoughFuncs_v1& hough_funcs) {
15  *this = hough_funcs;
16  return;
17 }
18 
19 void HelixHoughFuncs_v1::set_hough_space(HelixHoughSpace* hough_space) {
20 
21  _hough_space = hough_space;
22 }
23 
24 
25 void HelixHoughFuncs_v1::calculate_dzdl_range(float* hitpos3d, std::vector<float>& z0_range, std::vector<float>& kappa_phi_d_ranges, float* dzdl_range) {
26 
27  float x = hitpos3d[0];
28  float y = hitpos3d[1];
29  float z = hitpos3d[2];
30 // _hough_space->print_para_range();
31  // iz has to be an array of iz[zoomlevel], sot that it contains history
32  float z0_min = z0_range[0];
33  float z0_max = z0_range[1];
34  float kappa_min = kappa_phi_d_ranges[0];
35  float kappa_max = kappa_phi_d_ranges[1];
36  float phi_min = kappa_phi_d_ranges[2];
37  float phi_max = kappa_phi_d_ranges[3];
38 
39  float cosphi_min = cos(phi_min);
40  float cosphi_max = cos(phi_max);
41  float sinphi_min = sin(phi_min);
42  float sinphi_max = sin(phi_max);
43 
44  if (phi_max -phi_min > M_PI) { // for voting z only
45  cosphi_min = 1.;
46  cosphi_max = -1;
47  sinphi_min = 0.;
48  sinphi_max = 0.;
49  }
50 
51  float d_min = kappa_phi_d_ranges[4];
52  float d_max = kappa_phi_d_ranges[5];
53 
54 // cout<<"x "<<x <<" z0_min "<<z0_min<<" cosphi_min "<<cosphi_min<<" cosphi_max "<<cosphi_max<<endl;
55 // cout<<"d_min "<<d_min<<" d_max "<<d_max<<endl;
56 // cout<<"kappa_min "<<kappa_min<<" kappa_max "<<kappa_max<<endl;
57 // cout<<"phi_min "<<phi_min<<" phi_max "<<phi_max<<endl;
58  float d, k, dx, dy;
59  float dzdl[4]={-999,-999,-999,-999};
60  for (int icase=0; icase<2; icase++){
61  switch (icase){
62  case 0:
63  d= d_min; dx=d*cosphi_min; dy=d*sinphi_min;k=kappa_min;
64  break;
65  case 1:
66  d= d_max; dx=d*cosphi_max; dy=d*sinphi_max;k=kappa_max;
67  break;
68  }
69 
70  float D = sqrt(pow(x-dx,2)+pow(y-dy,2));
71  float a = k*D/2.;float a2 = a*a;
72  float s = -999.;;
73  if (a > .999) {a = 0.999; s = 2*asin(a)/k; }
74  else if (a < 0.01){ s = D + a2*D/6. + 3/40.*a2*a2*D + 5/112.*a2*a2*a2;}
75  else { s = 2*asin(a)/k;}
76 
77  switch (icase){
78  case 0:
79  dzdl[0] = pow(z-z0_max,2)/(pow(z-z0_max,2)+s*s);
80  dzdl[0] = sqrt(dzdl[0]);
81  if (z < z0_max) dzdl[0] *= -1;
82  dzdl[1] = pow(z-z0_min,2)/(pow(z-z0_min,2)+s*s);
83  dzdl[1] = sqrt(dzdl[1]);
84  if (z < z0_min) dzdl[1] *= -1;
85 // cout<<"dzdl[0] "<< dzdl[0]<<endl;
86 // cout<<"dzdl[1] "<<dzdl[1]<<endl;
87  break;
88 
89  case 1:
90  dzdl[2] = pow(z-z0_max,2)/(pow(z-z0_max,2)+s*s);
91  dzdl[2] = sqrt(dzdl[2]);
92  if (z < z0_max) dzdl[2] *= -1;
93  dzdl[3] = pow(z-z0_min,2)/(pow(z-z0_min,2)+s*s);
94  dzdl[3] = sqrt(dzdl[3]);
95  if (z < z0_min) dzdl[3] *= -1;
96 // cout<<"dzdl[2] "<<dzdl[2]<<endl;
97 // cout<<"dzdl[3] "<<dzdl[3]<<endl;
98  break;
99  }
100 // if (z0_min<5. && z0_max>5.){
101 // cout<<"x "<<x<<" y "<<y<< " z "<<z<<endl;
102 // cout<<"d "<<d<<" k "<<k<<" dx "<<dx<<" dy "<<dy<<" z-z0_max "<<z-z0_max<<" z-z0_min "<<z-z0_min<<endl;
103 // cout<<"D "<<D<<" s "<<s<<endl;
104 // }
105  }
106  float dzdl_min = -999;
107  float dzdl_max = -999;
108 
109 // if (z0_min<5. && z0_max>5.){
110 // cout<<"dzdl (dmin,kmin, phimin) max "<<dzdl[0]<<" dzdl (dmin,kmin,phimin) min "<<dzdl[1]<<endl;
111 // cout<<"dzdl (dmax,kmax, phimax) max "<<dzdl[2]<<" dzdl (dmax,kmax,phimax) min "<<dzdl[3]<<endl;
112 // }
113 
114  dzdl_min = dzdl[0];
115  for (int i = 1; i<4; i++) {
116  if (dzdl_min > dzdl[i])
117  dzdl_min = dzdl[i];
118  }
119  dzdl_max = dzdl[0];
120  for (int j = 1; j<4; j++) {
121  if (dzdl_max < dzdl[j])
122  dzdl_max = dzdl[j];
123  }
124 // if (z0_min<5. && z0_max>5.)
125 // if (_cur_zoom==1)
126 // cout<<"dzdl_min "<<dzdl_min<<" dzdl_max "<<dzdl_max<<endl;
127 
128 // dzdl_range[0]= _hough_space->get_dzdl_bin(_cur_zoom,dzdl_min); // min bin
129 // dzdl_range[1]= _hough_space->get_dzdl_bin(_cur_zoom,dzdl_max); // max bin
130  dzdl_range[0] = dzdl_min;
131  dzdl_range[1] = dzdl_max;
132 
133 }
134 
135 // if not separate by helicity
136 void HelixHoughFuncs_v1::calculate_phi_range(float* hitpos2d, std::vector<float>& kappa_d_ranges, float* phi_r_range, float* phi_l_range){
137 
138  float x = hitpos2d[0];
139  float y = hitpos2d[1];
140  float hitphi= atan2(y,x);
141  if (hitphi < 0.) hitphi += 2.*M_PI;
142 
143  float kmin = kappa_d_ranges[0];
144  float kmax = kappa_d_ranges[1];
145  float dmin = kappa_d_ranges[2];
146  float dmax = kappa_d_ranges[3];
147 
148 // cout<<"kmin "<<kmin<<" kmax "<<kmax<<" dmin "<< dmin<<" dmax "<<dmax<<endl;
149 
150  float d, k;
151  float D, Dinv, ak, hk, hksq, xk1, xk2, yk1,yk2;
152  float phi[8] = {-999,-999,-999,-999,-999,-999,-999,-999};
153 
154  for (int icase=0; icase<4; icase++ ){
155  switch (icase){
156  case 0:
157  d= dmin; k=kmin;
158  break;
159  case 1:
160  d= dmin; k=kmax;
161  break;
162  case 2:
163  d= dmax; k=kmin;
164  break;
165  case 3:
166  d= dmax; k=kmax;
167  break;
168  }
169 
170  // calc phi business
171  D = sqrt(x*x+y*y);
172  Dinv = 1./D;
173  ak = (2.*d + d*d*k + D*D*k)/2.*Dinv;
174  hksq = pow(d*k+1.,2) - ak*ak;
175  if (hksq>=0.) {
176  hk=sqrt(hksq);
177  xk1 = (ak*x + hk*y)*Dinv;
178  xk2 = (ak*x - hk*y)*Dinv;
179  yk1 = (ak*y - hk*x)*Dinv;
180  yk2 = (ak*y + hk*x)*Dinv;
181  }
182 
183  switch (icase){
184  case 0:
185  if (hksq<0.){
186  phi[0]= hitphi;//right
187  phi[1]= hitphi;//left
188  }else {
189  phi[0]= atan2(yk1,xk1);// right
190  if (phi[0]<0.) phi[0]+= 2.*M_PI;
191  phi[1]= atan2(yk2,xk2);// left
192  if (phi[1]<0.) phi[1]+= 2.*M_PI;
193  }
194  break;
195 
196  case 1:
197  if (hksq<0.){
198  phi[2]= hitphi;
199  phi[3]= hitphi;
200  }else {
201  phi[2]=atan2(yk1,xk1);// right
202  if (phi[2]<0.) phi[2]+= 2.*M_PI;
203  phi[3]= atan2(yk2,xk2);// left
204  if (phi[3]<0.) phi[3]+= 2.*M_PI;
205  }
206  break;
207 
208  case 2:
209  if (hksq<0.){
210  phi[4]= hitphi;
211  phi[5]= hitphi;
212  }else{
213  phi[4]=atan2(yk1,xk1);//right
214  if (phi[4]<0.) phi[4]+= 2.*M_PI;
215  phi[5]=atan2(yk2,xk2);//left
216  if (phi[5]<0.) phi[5]+= 2.*M_PI;
217  }
218  break;
219 
220  case 3:
221  if (hksq<0.){
222  phi[6]= hitphi;
223  phi[7]= hitphi;
224  } else{
225  phi[6]=atan2(yk1,xk1);//right
226  if (phi[6]<0.) phi[6]+= 2.*M_PI;
227  phi[7]=atan2(yk2,xk2);//left
228  if (phi[7]<0.) phi[7]+= 2.*M_PI;
229  }
230  break;
231  }
232 
233  }
234 // for (int i = 0; i<8;i++)
235 // cout<<"phi["<<i<<"]="<<phi[i]<<endl;
236  float twopi = 2.*M_PI;
237  float piover2=M_PI/2.;
238  float threepiover2=3.*M_PI/2.;
239  bool first = (phi[0]<piover2) || (phi[2]<piover2) || (phi[4]<piover2 || (phi[6]<piover2));
240  bool fourth = (phi[0]>threepiover2) || (phi[2]>threepiover2) || (phi[4]>threepiover2) || (phi[6]>threepiover2);
241 
242  if (first && fourth){
243  if (phi[0]>threepiover2) phi[0] -= twopi;
244  if (phi[2]>threepiover2) phi[2] -= twopi;
245  if (phi[4]>threepiover2) phi[4] -= twopi;
246  if (phi[6]>threepiover2) phi[6] -= twopi;
247  }
248 
249  float phi_r_min = phi[0];
250  for (int i =1; i<4; i++){
251  if(phi[2*i]<phi_r_min) phi_r_min= phi[2*i];
252  }
253  float phi_r_max = phi[0];
254  for (int i =1; i<4; i++){
255  if(phi[2*i]>phi_r_max) phi_r_max=phi[2*i];
256  }
257 
258  phi_r_range[0] = phi_r_min;
259  phi_r_range[1] = phi_r_max;
260 
261 
262  first = (phi[1]<piover2) || (phi[3]<piover2) || (phi[5]<piover2 || (phi[7]<piover2));
263  fourth = (phi[1]>threepiover2) || (phi[3]>threepiover2) || (phi[5]>threepiover2) || (phi[7]>threepiover2);
264 
265  if (first && fourth){
266  if (phi[1]>threepiover2) phi[1] -= twopi;
267  if (phi[3]>threepiover2) phi[3] -= twopi;
268  if (phi[5]>threepiover2) phi[5] -= twopi;
269  if (phi[7]>threepiover2) phi[7] -= twopi;
270  }
271 
272  float phi_l_min = phi[1];
273  for (int i =1; i<4; i++){
274  if(phi[2*i+1]<phi_l_min) phi_l_min= phi[2*i+1];
275  }
276  float phi_l_max = phi[1];
277  for (int i =1; i<4; i++){
278  if(phi[2*i+1]>phi_l_max) phi_l_max=phi[2*i+1];
279  }
280 
281  phi_l_range[0] = phi_l_min;
282  phi_l_range[1] = phi_l_max;
283 
284 // if (dmin<.02 && dmax>0.02)
285 // cout<<"phi_l_min "<<phi_l_min<<" phi_l_max "<<phi_l_max<<" phi_r_min "<<phi_r_min<<" phi_r_max "<<phi_r_max<<endl;
286 
287 }
288 
289 // if separage by helicity
290 // initial k
291 void HelixHoughFuncs_v1::calculate_phi_range(float* hitpos2d, std::vector<float>& kappa_d_ranges, int helicity, float* phi_range, float* phi_next_range){
292 
293  float x = hitpos2d[0];
294  float y = hitpos2d[1];
295  float hitphi= atan2(y,x);
296  if (hitphi < 0.) hitphi += 2.*M_PI;
297 
298  float kmin = kappa_d_ranges[0];
299  float kmax = kappa_d_ranges[1];
300  float dmin = kappa_d_ranges[2];
301  float dmax = kappa_d_ranges[3];
302 
303  float d, k;
304  float D, Dinv, ak, hk=0., hksq, xk1, xk2, xk, yk1,yk2, yk;
305  float cross, phi[4];
306 
307  for (int icase=0; icase<4; icase++ ){
308  switch (icase){
309  case 0:
310  d= dmin; k=kmin;
311  break;
312  case 1:
313  d= dmin; k=kmax;
314  break;
315  case 2:
316  d= dmax; k=kmin;
317  break;
318  case 3:
319  d= dmax; k=kmax;
320  break;
321  }
322 
323  // calc phi business
324  D = sqrt(x*x+y*y);
325  Dinv = 1./D;
326  ak = (2.*d + d*d*k + D*D*k)/(2.*D);
327  hksq = pow(d*k+1.,2) - ak*ak;
328  if (hksq>=0.)
329  hk=sqrt(hksq);
330  xk1 = (ak*x + hk*y)*Dinv;
331  xk2 = (ak*x - hk*y)*Dinv;
332  yk1 = (ak*y - hk*x)*Dinv;
333  yk2 = (ak*y + hk*x)*Dinv;
334 
335 
336  cross = x*yk1 - y*xk1;
337  if (cross*helicity>0){
338  xk=xk1;
339  yk=yk1;
340  } else {
341  xk=xk2;
342  yk=yk2;
343  }
344 
345  switch (icase){
346  case 0:
347  if (hksq<0.){
348  phi[0]= hitphi;//left
349  }else {
350  phi[0]= atan2(yk,xk);// left
351  if (phi[0]<0.) phi[0]+= 2.*M_PI;
352  }
353  break;
354 
355  case 1: // min d, max k next_range[0]
356  if (hksq<0.){
357  phi[1]= hitphi;
358  }else {
359  phi[1]= atan2(yk,xk);// left
360  if (phi[1]<0.) phi[1]+= 2.*M_PI;
361  }
362  break;
363 
364  case 2:
365  if (hksq<0.){
366  phi[2]= hitphi;
367  }else{
368  phi[2]=atan2(yk,xk);
369  if (phi[2]<0.) phi[2]+= 2.*M_PI;
370  }
371  break;
372 
373  case 3: // max d, max k phi_next_range[1]
374  if (hksq<0.){
375  phi[3]= hitphi;
376  } else{
377  phi[3]=atan2(yk,xk);
378  if (phi[3]<0.) phi[3]+= 2.*M_PI;
379  }
380  break;
381  }
382 
383  }// cases
384 
385  phi_next_range[0] = phi[1];
386  phi_next_range[1] = phi[3];
387 
388  float twopi = 2.*M_PI;
389  float piover2=M_PI/2.;
390  float threepiover2=3.*M_PI/2.;
391  bool first = (phi[0]<piover2) || (phi[1]<piover2) || (phi[2]<piover2 || (phi[3]<piover2));
392  bool fourth = (phi[0]>threepiover2) || (phi[1]>threepiover2) || (phi[2]>threepiover2) || (phi[3]>threepiover2);
393 
394  if (first && fourth){
395  if (phi[0]>threepiover2) phi[0] -= twopi;
396  if (phi[1]>threepiover2) phi[1] -= twopi;
397  if (phi[2]>threepiover2) phi[2] -= twopi;
398  if (phi[3]>threepiover2) phi[3] -= twopi;
399  }
400 
401  float phi_min = phi[0];
402  for (int i =0; i<4; i++){
403  if(phi[i]<phi_min) phi_min= phi[i];
404  }
405  float phi_max = phi[0];
406  for (int i =0; i<4; i++){
407  if(phi[i]>phi_max) phi_max=phi[i];
408  }
409 
410  phi_range[0] = phi_min;
411  phi_range[1] = phi_max;
412 
413 }
414 
415 
416 void HelixHoughFuncs_v1::calculate_phi_range(float* hitpos2d, std::vector<float>& kappa_d_ranges, int helicity, float* phi_range, float* phi_prev_range, float* phi_next_range){
417 
418  float x = hitpos2d[0];
419  float y = hitpos2d[1];
420  float hitphi= atan2(y,x);
421  if (hitphi < 0.) hitphi += 2.*M_PI;
422 
423 // float kmin = kappa_d_ranges[0];
424  float kmax = kappa_d_ranges[1];
425  float dmin = kappa_d_ranges[2];
426  float dmax = kappa_d_ranges[3];
427 
428  float d, k;
429  float D, Dinv, ak, hk=0., hksq, xk1, xk2, xk, yk1,yk2, yk;
430  float cross, phi[4];
431 
432 
433  for (int icase=0; icase<2; icase++ ){
434  switch (icase){
435  case 0:
436  d= dmin; k=kmax;
437  break;
438  case 1:
439  d= dmax; k=kmax;
440  break;
441  }
442 
443  // calc phi business
444  D = sqrt(x*x+y*y);
445  Dinv = 1./D;
446  ak = (2.*d + d*d*k + D*D*k)/(2.*D);
447  hksq = pow(d*k+1.,2) - ak*ak;
448  if (hksq>=0.)
449  hk=sqrt(hksq);
450  xk1 = (ak*x + hk*y)*Dinv;
451  xk2 = (ak*x - hk*y)*Dinv;
452  yk1 = (ak*y - hk*x)*Dinv;
453  yk2 = (ak*y + hk*x)*Dinv;
454 
455 
456  cross = x*yk1 - y*xk1;
457  if (cross*helicity>0){
458  xk=xk1;
459  yk=yk1;
460  } else {
461  xk=xk2;
462  yk=yk2;
463  }
464 
465  switch (icase){
466  case 0:
467  if (hksq<0.){
468  phi[0]= hitphi;//left
469  }else {
470  phi[0]= atan2(yk,xk);// left
471  if (phi[0]<0.) phi[0]+= 2.*M_PI;
472  }
473  break;
474 
475  case 1: // min d, max k next_range[0]
476  if (hksq<0.){
477  phi[1]= hitphi;
478  }else {
479  phi[1]= atan2(yk,xk);// left
480  if (phi[1]<0.) phi[1]+= 2.*M_PI;
481  }
482  break;
483 
484  }
485 
486  }// cases
487 
488  phi_next_range[0] = phi[0];
489  phi_next_range[1] = phi[1];
490  phi[2] = phi_prev_range[0];
491  phi[3] = phi_prev_range[1];
492 
493  float twopi = 2.*M_PI;
494  float piover2=M_PI/2.;
495  float threepiover2=3.*M_PI/2.;
496  bool first = (phi[0]<piover2) || (phi[1]<piover2) || (phi[2]<piover2 || (phi[3]<piover2));
497  bool fourth = (phi[0]>threepiover2) || (phi[1]>threepiover2) || (phi[2]>threepiover2) || (phi[3]>threepiover2);
498 
499  if (first && fourth){
500  if (phi[0]>threepiover2) phi[0] -= twopi;
501  if (phi[1]>threepiover2) phi[1] -= twopi;
502  if (phi[2]>threepiover2) phi[2] -= twopi;
503  if (phi[3]>threepiover2) phi[3] -= twopi;
504  }
505 
506  float phi_min = phi[0];
507  for (int i =0; i<4; i++){
508  if(phi[i]<phi_min) phi_min= phi[i];
509  }
510  float phi_max = phi[0];
511  for (int i =0; i<4; i++){
512  if(phi[i]>phi_max) phi_max=phi[i];
513  }
514 
515  phi_range[0] = phi_min;
516  phi_range[1] = phi_max;
517 
518 }
519 /*
520 
521 float HelixHoughFuncs_v1::calc_dzdl_error() {
522 float dzdl_err = ;
523 
524 return dzdl_err;
525 }
526 
527 float HelixHoughFuncs_v1::calc_phi_error() {
528 float phi_err =;
529 
530 return phi_err;
531 }
532 
533 */