246 {
247 HepVector cen1(lp1.
center());
248 HepVector cen2(lp2.
center());
249 double dx = cen1(1)-cen2(1);
250 double dy = cen1(2)-cen2(2);
251 double dc = sqrt(dx*dx+dy*dy);
253 double a1 = std::sqrt(lp1.alpha()) + std::sqrt(lp1.beta());
254 double a2 = std::sqrt(lp2.alpha()) + std::sqrt(lp2.beta());
255 double a3 = lp1.alpha()*lp2.alpha() + lp1.beta()*lp2.beta();
256 double det = lp1.alpha()*lp2.beta() - lp1.beta()*lp2.alpha();
257 if(fabs(det)>1e-12) {
258 double c1 = a2 * std::sqrt(lp1.
kappa()) + a1 * std::sqrt(lp2.
kappa()) -
261 double cinv = 1.0 /
c1;
262 double c2 = std::sqrt(a3) - 0.5 * (a1 + a2) - 2.0 * a3 *
263 (lp1.gamma() * lp2.
kappa() + lp2.gamma() * lp1.
kappa());
264 double c3 = a2 * std::sqrt(lp1.gamma()) + a1 * std::sqrt(lp2.gamma()) -
265 2.0 * a3 * lp1.gamma() * lp2.gamma();
266 double root = std::sqrt(c2) - 4.0 *
c1 * c3;
269 double rad2[2];
270 rad2[0] = 0.5 * cinv * (-c2 -
root);
271 rad2[1] = 0.5 * cinv * (-c2 +
root);
272 double ab1 = -(lp2.beta() * lp1.gamma() - lp1.beta() * lp2.gamma());
273 double ab2 = (lp2.alpha() * lp1.gamma() - lp1.alpha() * lp2.gamma());
274 double ac1 = -(lp2.beta() * lp1.
kappa() - lp1.beta() * lp2.
kappa());
275 double ac2 = (lp2.alpha() * lp1.
kappa() - lp1.alpha() * lp2.
kappa());
276 double dinv = 1.0 / det;
277 v1(1) = dinv * (ab1 + ac1 * rad2[0]);
278 v1(2) = dinv * (ab2 + ac2 * rad2[0]);
279 v1(3) = 0;
280 v2(1) = dinv * (ab1 + ac1 * rad2[1]);
281 v2(2) = dinv * (ab2 + ac2 * rad2[1]);
282 v2(3) = 0;
283 double d1 = lp1.
d(v1(1),v1(2));
284 double d2 = lp2.
d(v1(1),v1(2));
285 double d3 = lp1.
d(v2(1),v2(2));
286 double d4 = lp2.
d(v2(1),v2(2));
287 double r = sqrt(rad2[0]);
288 Lpar::Cpar cp1(lp1);
289 Lpar::Cpar cp2(lp2);
290 for(int j=0;j<2;j++) {
291 double s1,s2;
292 if(j==0) {
293 s1 = lp1.
s(v1(1),v1(2));
294 s2 = lp2.
s(v1(1),v1(2));
295 } else {
296 s1 = lp1.
s(v2(1),v2(2));
297 s2 = lp2.
s(v2(1),v2(2));
298 }
299 double phi1 = cp1.fi() + 2 * cp1.cu() * s1;
300 double phi2 = cp2.fi() + 2 * cp2.cu() * s2;
301 double f = (1 + 2 * cp1.cu() * cp1.da()) *
302 (1 + 2 * cp2.cu() * cp2.da()) *
cos(cp1.fi()-cp2.fi());
303 f -= 2 * (lp1.gamma() * lp2.
kappa() + lp2.gamma() * lp1.
kappa());
304 double cosphi12 = f;
305 }
306 return 2;
307 }
308 }
309 }
310 }
311 return 0;
312}
double cos(const BesAngle a)
double d(double x, double y) const
double s(double x, double y) const