10 double cx_det = -vertex_x;
11 double cy_det = -vertex_y;
13 double d2 = ((cx-cx_det)*(cx-cx_det) + (cy-cy_det)*(cy-cy_det));
15 if(d > (rad_det + rad_trk))
19 if(d < fabs(rad_det - rad_trk))
24 double r2 = rad_trk*rad_trk;
27 double R2 = rad_det*rad_det;
28 double a = 0.5*(R2 - r2 +
d2)*d_inv;
30 double P2x = cx_det + (cx-cx_det)*h;
31 double P2y = cy_det + (cy-cy_det)*h;;
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;
40 double P3x2 = P2x + ux*
h;
41 double P3y2 = P2y + uy*
h;
43 double d1_2 = (startx - P3x1)*(startx - P3x1) + (starty - P3y1)*(starty - P3y1);
44 double d2_2 = (startx - P3x2)*(startx - P3x2) + (starty - P3y2)*(starty - P3y2);
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)
63 for(
unsigned int i=0;i<detrad.size();++i)
67 init_ZHough(z0_nbin, theta_nbin, minzres, maxzres, zrange);
79 if(zhough!=NULL){
delete zhough;}
80 zhough =
new ZHough_Cylindrical(z0_nbin, theta_nbin, minzres, maxzres, zrange, sigma_xy, sigma_z);
82 vector<double> lxy, lz;
85 lxy.push_back(0.005/3.);
86 lz.push_back(0.0425/3.);
91 zhough->setCircleHough(
this);
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)
112 for(
unsigned int i=0;i<hits.size();i++)
114 AngleIndexPair temppair(atan2(hits[i].get_y(), hits[i].get_x()),hits[i].index);
115 angle_list[hits[i].layer].addPair( temppair );
125 void CylindricalHough::addHits(
unsigned int zlevel, vector<SimpleTrack3D>& temptracks, vector<SimpleTrack3D>& tracks, vector<float>& params,
int tracks_per_hit,
float z_cut)
127 vector<double> chi2_hit;
129 float chi2_cut = params[1];
130 unsigned int min_hits = (
unsigned int)(params[2]);
131 float max_kappa_cut = params[3];
133 for(
unsigned int trk=0;trk<temptracks.size();++trk)
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);
144 if(fabs(temptracks[trk].kappa) < 1.0
e-8){temptracks[trk].kappa = 1.0e-8;}
146 r = 1./temptracks[trk].kappa;
149 cx = (temptracks[trk].d+
r)*cos(temptracks[trk].
phi);
150 cy = (temptracks[trk].d+
r)*sin(temptracks[trk].phi);
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;
162 if(cross<0.){helicity=
false;}
165 vector<bool> layer_used;
167 for(
unsigned int h=0;
h<temptracks[trk].hits.size();++
h)
169 layer_used[temptracks[trk].hits[
h].layer]=
true;
171 vector<int> unused_layers;
172 for(
unsigned int ll=0;ll<layer_used.size();++ll)
174 if(layer_used[ll]==
false){unused_layers.push_back(ll);}
177 for(
unsigned int ll=0;ll<unused_layers.size();ll++)
179 if(((unused_layers.size() - ll) + temptracks[trk].
hits.size()) < min_hits){
break;}
181 int layer = unused_layers[ll];
182 int closest_layer = temptracks[trk].hits[0].layer;
186 for(
int cl=1;cl<(
int)(temptracks[trk].
hits.size());cl++)
188 if( fabs(layer - temptracks[trk].
hits[cl].layer) < fabs(layer - closest_layer) )
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();
197 double x_intersect=0.;
198 double y_intersect=0.;
199 double phi_error =
phicut;
205 bool try_helicity = helicity;
209 if( ( intersected==
false ) || ( x_intersect != x_intersect ) || ( y_intersect != y_intersect ) )
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;
219 double z_intersect=0.;
223 if(v >= 0.999999){v = 0.999999;}
229 double temp2 = D*0.5;
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)
247 vector<AngleIndexPair*> hit_candidates;
248 angle_list[
layer].getRangeList(intersect_phi, phi_error, hit_candidates);
250 for(
unsigned int hit_cand=0;hit_cand<hit_candidates.size();hit_cand++)
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 )
262 temptracks[trk].hits.push_back((*(hits_vec[0]))[hit_candidates[hit_cand]->index]);
264 double chi2 = zhough->fitTrack(temptracks[trk], chi2_hit);
266 if( (chi2 != chi2) || (chi2 > chi2_cut) || (temptracks[trk].kappa > max_kappa_cut))
268 temptracks[trk] = temptrack;
272 temptrack = temptracks[trk];
277 if(temptracks[trk].
hits.size() >= min_hits)
280 for(
unsigned int itrack=0; itrack<tracks.size();itrack++)
283 sort (tracks[itrack].
hits.begin(), tracks[itrack].hits.end(),SimpleHit3D_LessThan);
287 sort (tracks.begin(), tracks.end(), SimpleTrack3D_LessThan);
289 sort (temptracks[trk].
hits.begin(), temptracks[trk].hits.end(),SimpleHit3D_LessThan);
292 if(binary_search(tracks.begin(), tracks.end(), temptracks[trk], SimpleTrack3D_LessThan) ==
false)
295 tracks.push_back(temptracks[trk]);
297 for(
unsigned int h=0;
h<temptracks[trk].hits.size();
h++)
302 used_vec[temptracks[trk].hits[
h].index]++;