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