1 /*
2 * Find the circle passing through two points and a fixed impact parameter.
3 * Typically, there are two circles that satisfy these criteria. The exception
4 * is when a straight line connecting the two points has the specified impact
5 * parameter, where there is only one solution.
6 *
7 * Optionally, you can provide a minimum circle radius. If neither circle
8 * satisfies the minimum circle radius cut, the fit fails. If one circle
9 * satisfies the minimum radius cut, but the other solution does not, then a
10 * second solution is found by setting the radius to the minimum radius and
11 * allowing the impact parameter to be smaller than the specified value.
12 *
13 * This class is used by the seedtracker track finding algorithm to check if
14 * a two hit trial seed satisfies the specified pT and impact parameter cuts.
15 * The two solutions found will generally correspond to the minimum and maximum
16 * radius circles given the two hits and the impact parameter cut. However,
17 * if the two hits are consistent with a straight-line with an impact parameter
18 * less than or equal to the specified value, there is no upper limit on the
19 * circle radius. A method is provided to check for this special case and
20 * return the straight-line parameters.
21 */
22
23 package org.lcsim.fit.twopointcircle;
24
25 import java.util.ArrayList;
26 import java.util.List;
27
28 import org.lcsim.event.TrackerHit;
29
30 /**
31 * @author Richard Partridge
32 */
33 public class TwoPointCircleFitter {
34 private double _rmin;
35 private List<TwoPointCircleFit> _circlefits;
36 private TwoPointLineFit _linefit;
37 private boolean _debug = false;
38 private static double _eps = 1e-6;
39 private static double twopi = 2. * Math.PI;
40
41 /**
42 * Constructor specifying a minimum radius. If the minimum radius is >0,
43 * the algorithm will enforce this constraint even if it means having a
44 * smaller impact parameter or fewer / no successful fits.
45 *
46 * @param rmin minimum circle radius
47 */
48 public TwoPointCircleFitter(double rmin) {
49 _rmin = rmin;
50 _circlefits = new ArrayList<TwoPointCircleFit>();
51 }
52
53 /**
54 * Constructor with no minimum radius.
55 */
56 public TwoPointCircleFitter() {
57 this(0.);
58 }
59
60 /**
61 * Fit a circle given two TrackerHits and the impact parameter.
62 *
63 * @param hit1 hit #1
64 * @param hit2 hit #2
65 * @param dmax impact parameter
66 * @return fit status - true if at least one circle fit is found
67 */
68 public boolean FitCircle(TrackerHit hit1, TrackerHit hit2, double dmax) {
69 double[] pos1 = hit1.getPosition();
70 double[] pos2 = hit2.getPosition();
71 return FitCircle(pos1[0], pos1[1], pos2[0], pos2[1], dmax);
72 }
73
74 /**
75 * Fit a circle given coordinates for two points and the impact parameter.
76 *
77 * @param x1 x coordinate of first point
78 * @param y1 y coordinate of first point
79 * @param x2 x coordinate of second point
80 * @param y2 y coordinate of second point
81 * @param dmax impact parameter
82 * @return fit status - true if at least one circle fit is found
83 */
84 public boolean FitCircle(double x1, double y1, double x2, double y2, double dmax) {
85
86 // Clear the list of fits from the previous time the method was called
87 _circlefits.clear();
88 _linefit = null;
89
90 // Calculate some useful quantities
91 double r1sq = x1*x1 + y1*y1;
92 double r2sq = x2*x2 + y2*y2;
93 double dmaxsq = dmax*dmax;
94
95 // If either point is inside maximum IP cut, throw an exception since
96 // we are under-constrained.
97 if (r1sq < dmaxsq || r2sq < dmaxsq)
98 throw new RuntimeException("Cannot handle case where hit is inside maximum impact parameter cut");
99
100 // Get the u = (r2 - r1)/|r2 - r1| unit vector
101 double ux = x2 - x1;
102 double uy = y2 - y1;
103 double u = Math.sqrt(ux*ux + uy*uy);
104 ux = ux / u;
105 uy = uy / u;
106 if (u<_eps)return false;
107
108 // get the midpoint vector
109 double mx = 0.5 * (x1 + x2);
110 double my = 0.5 * (y1 + y2);
111 double msq = mx*mx + my*my;
112
113 // Get the unit vector normal to u
114 double vx = uy;
115 double vy = -ux;
116
117 // Get the impact parameter for an infinite momentum track with hits 1 & 2
118 double ipinf = (x1*y2 - y1*x2) / u;
119
120 // Create an array to hold alpha, the distance of the circle center from the midpoint between hits 1 and 2
121 int nalpha = 2;
122 double alpha[] = new double[nalpha];
123
124 // First calculate the denominator used in the calculation of alpha
125 double denom = 2. * (ipinf*ipinf - dmaxsq);
126
127 // Check if we are consistent with a straight-line track
128 boolean sltrack = denom < _eps;
129
130 // Check for the singular case that occurs when a straight-line track going through the two hits
131 // is tangent to a circle of radius _dMax
132 if (Math.abs(denom) < _eps) {
133
134 // Singular case - only one finite momentum solution
135 double beta = msq - 0.25 * u*u - dmaxsq;
136 nalpha = 1;
137 alpha[0] = (u*u * dmaxsq - beta*beta) / (4. * beta * ipinf);
138
139 } else {
140
141 // Non-singular case - find the two solutions for alpha
142 double r1dotr2 = x1*x2 + y1*y2;
143 double term1 = -ipinf * (r1dotr2 - dmaxsq) / denom;
144 double term2 = Math.abs(dmax * Math.sqrt((r1sq - dmaxsq) * (r2sq - dmaxsq)) / denom);
145
146 // Make sure alpha[0] is the solution with the largest circle radius
147 if (term1 < 0.) term2 = -term2;
148 alpha[0] = term1 + term2;
149 alpha[1] = term1 - term2;
150 }
151
152 // Loop over the two solutions
153 for (int i = 0; i<nalpha; i++) {
154
155 // Find the circle radius and check if it exceeds the minimum value
156 double rcurv = Math.sqrt(alpha[i]*alpha[i] + 0.25 * u*u);
157 if (rcurv < _rmin) {
158
159 // The first iteration has the largest circle radius - if it doesn't pass the cut and the track
160 // isn't consistent with a straight-line track, we can't form a circle that passes the pT cut
161 if (i == 0 && !sltrack) return false;
162
163 // If we get here, either the first solution passed the pT cut or we are consistent with a
164 // straight-line track. Find the circle with the minimum radius and the new alpha
165 rcurv = _rmin;
166 double newalpha = Math.sqrt(_rmin*_rmin - 0.25 * u*u);
167
168 // Figure out the sign of alpha. If we are consistent with a
169 // straight-line track, then the sign should be opposite for the
170 // second (low radius) solution. If we are not consistent with
171 // a straight-line track, the second solution should have the same
172 // sign as the first solution.
173 int isign = 1;
174 if (i == 1 && sltrack && alpha[0] >= 0.) isign = -1;
175 if (i == 1 && !sltrack && alpha[0] < 0.) isign = -1;
176 alpha[i] = isign * newalpha;
177 }
178
179 // Find the center of the circle
180 double xc = mx + alpha[i] * vx;
181 double yc = my + alpha[i] * vy;
182 double rc = Math.sqrt(xc*xc + yc*yc);
183
184 // Find the point of closest approach
185 double x0 = xc * (1. - rcurv / rc);
186 double y0 = yc * (1. - rcurv / rc);
187
188 // Make some checks if we have debugging turned on
189 if (_debug) {
190 double c1 = Math.sqrt((x1-xc)*(x1-xc)+(y1-yc)*(y1-yc));
191 if (Math.abs(c1-rcurv) > _eps * c1) throw new RuntimeException("Error in circle finding - c1 = "+c1+" rcurv = "+rcurv);
192 double c2 = Math.sqrt((x2-xc)*(x2-xc)+(y2-yc)*(y2-yc));
193 if (Math.abs(c2-rcurv) > _eps * c2) throw new RuntimeException("Error in circle finding - c2 = "+c2+" rcurv = "+rcurv);
194 double c3 = Math.sqrt((x0-xc)*(x0-xc)+(y0-yc)*(y0-yc));
195 if (Math.abs(c3-rcurv) > _eps * c3) throw new RuntimeException("Error in circle finding - c3 = "+c3+" rcurv = "+rcurv);
196 double r0 = Math.sqrt(x0*x0 + y0*y0);
197 if (r0 > dmax+100*_eps) throw new RuntimeException("Invalid DCA point for solution "+i+" r0: "+r0+" dmax: "+dmax);
198 if (rcurv < _rmin) throw new RuntimeException("Invalid circle radius for solution "+i+" rcurv: "+rcurv+" rmin: "+_rmin);
199 }
200
201 // Find the x-y arc lengths to hits 1 and 2
202 // First find azimuthal angles for the dca and hit positions relative to the circle center
203 double phi0 = Math.atan2(y0-yc, x0-xc);
204 double phi1 = Math.atan2(y1-yc, x1-xc);
205 double phi2 = Math.atan2(y2-yc, x2-xc);
206
207 // Find the angle between the hits and the DCA under the assumption that |dphi| < pi
208 double dphi1 = phi1 - phi0;
209 if (dphi1 > Math.PI) dphi1 -= twopi;
210 if (dphi1 < -Math.PI) dphi1 += twopi;
211 double dphi2 = phi2 - phi0;
212 if (dphi2 > Math.PI) dphi2 -= twopi;
213 if (dphi2 < -Math.PI) dphi2 += twopi;
214
215 // Find the hit closest to the DCA and use this hit to determine the circle "direction"
216 boolean cw;
217 if (Math.abs(dphi1) < Math.abs(dphi2)) cw = dphi1 < 0.;
218 else cw = dphi2 < 0.;
219
220 // Find the arc lengths to points 1 and 2
221 double s1 = -dphi1 * rcurv;
222 double s2 = -dphi2 * rcurv;
223 if (!cw) {
224 s1 = -s1;
225 s2 = -s2;
226 }
227
228 // Fix the case when dphi1 & dphi2 have opposite signs, making an arc length negative
229 if (s1 < 0.) s1 += twopi * rcurv;
230 if (s2 < 0.) s2 += twopi * rcurv;
231
232 _circlefits.add(new TwoPointCircleFit(xc, yc, rcurv, cw, s1, s2));
233 }
234
235 // If we are consistent with a straight-line track, calculate the straight-line
236 // distances to the hits
237 if (sltrack) {
238
239 // Find the distances from the straight-line DCA to the hits
240 double s1 = (x1*ux + y1*uy);
241 double s2 = (x2*ux + y2*uy);
242
243 // If these distances have opposite sign, we violate causality for a hit originating at the DCA
244 if (s1*s2 < 0.) {
245
246 // Calculate the parameters of the line fit
247 double x0 = x1 - s1 * ux;
248 double y0 = y1 - s1 * uy;
249 double phi = Math.atan2(uy, ux);
250
251 // Flip the signs of the path lengths and direction to make path lengths positive
252 if (s1 < 0.) {
253 s1 = -s1;
254 s2 = -s2;
255 phi = phi + Math.PI;
256 if (phi > Math.PI) phi += twopi;
257 }
258
259 // Save a new TwoPointLineFit
260 _linefit = new TwoPointLineFit(x0, y0, phi, s1, s2);
261 }
262
263 // We should always find a valid circle above for this case - check to make sure
264 if (_debug & _circlefits.size() == 0) throw new RuntimeException("No circle found for hits consistent with infinite momentum");
265 }
266
267 return _circlefits.size() > 0;
268 }
269
270 /**
271 * Get the list of TwoPointCircleFits that are found.
272 *
273 * @return list of circle fits
274 */
275 public List<TwoPointCircleFit> getCircleFits() {
276 return _circlefits;
277 }
278
279 /**
280 * Get the TwoPointLineFit if the two hits are consistent with a straight-line
281 * track that passes within the circle defined by the impact parameter cut.
282 * If the line connecting the two hits has a larger impact parameter, then
283 * a null pointer is returned.
284 *
285 * @return pointer to line fit (null if no line fit is found)
286 */
287 public TwoPointLineFit getLineFit() {
288 return _linefit;
289 }
290
291 /**
292 * Set the minimum circle radius.
293 *
294 * @param rmin minimum radius
295 */
296 public void setRMin(double rmin) {
297 _rmin = rmin;
298 }
299
300 /**
301 * Turn on/off the debugging checks.
302 *
303 * @param debug state of debug flag
304 */
305 public void setDebug(boolean debug) {
306 _debug = debug;
307 }
308 }