7 IMPLICIT DOUBLE PRECISION(d)
8 common/lujets/
n,
k(9000,5),
p(9000,5),
v(9000,5)
10 common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
12 common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
14 dimension dps(5),kfl(3),pmq(3),
px(3),
py(3),gam(3),ie(2),pr(2),
15 &
in(9),dhm(4),dhg(4),dp(5,5),irank(2),mju(4),iju(3),pju(5,5),
16 &tju(5),kfjh(2),njs(2),kfjs(2),pjs(4,5)
20 dfour(
i,
j)=dp(
i,4)*dp(
j,4)-dp(
i,1)*dp(
j,1)-dp(
i,2)*dp(
j,2)-
34 IF(
i.GT.min(
n,mstu(4)-mstu(32)))
THEN
35 CALL
luerrm(12,
'(LUSTRF:) failed to reconstruct jet system')
36 IF(mstu(21).GE.1)
RETURN
38 IF(
k(
i,1).NE.1.AND.
k(
i,1).NE.2.AND.
k(
i,1).NE.41) goto 110
41 kq=kchg(kc,2)*isign(1,
k(
i,2))
43 IF(
n+5*
np+11.GT.mstu(4)-mstu(32)-5)
THEN
44 CALL
luerrm(11,
'(LUSTRF:) no more memory left in LUJETS')
45 IF(mstu(21).GE.1)
RETURN
53 120 dps(
j)=dps(
j)+
p(
i,
j)
60 IF(kq.NE.2) kqsum=kqsum+kq
63 IF(kqsum.EQ.kq) mju(1)=
n+
np
64 IF(kqsum.NE.kq) mju(2)=
n+
np
66 IF(
k(
i,1).EQ.2.OR.
k(
i,1).EQ.41) goto 110
68 CALL
luerrm(12,
'(LUSTRF:) unphysical flavour combination')
69 IF(mstu(21).GE.1)
RETURN
73 CALL ludbrb(
n+1,
n+
np,0.,0.,-dps(1)/dps(4),-dps(2)/dps(4),
86 IF(
i.EQ.
n+nr.AND.iabs(
k(
n+1,2)).NE.21) goto 140
89 IF(
k(
i,1).EQ.41.OR.
k(
i1,1).EQ.41) goto 140
90 IF(mju(1).NE.0.AND.
i1.LT.mju(1).AND.iabs(
k(
i1,2)).NE.21)
92 IF(mju(2).NE.0.AND.
i.GT.mju(2).AND.iabs(
k(
i,2)).NE.21) goto 140
93 pap=sqrt((
p(
i,1)**2+
p(
i,2)**2+
p(
i,3)**2)*(
p(
i1,1)**2+
96 pdr=4.*(pap-pvp)**2/(paru13**2*pap+2.*(pap-pvp))
97 IF(
pdr.LT.pdrmin)
THEN
104 IF(pdrmin.LT.paru12.AND.ir.EQ.
n+nr)
THEN
107 p(
n+1,5)=sqrt(
max(0.,
p(
n+1,4)**2-
p(
n+1,1)**2-
p(
n+1,2)**2-
111 ELSEIF(pdrmin.LT.paru12)
THEN
113 160
p(ir,
j)=
p(ir,
j)+
p(ir+1,
j)
114 p(ir,5)=sqrt(
max(0.,
p(ir,4)**2-
p(ir,1)**2-
p(ir,2)**2-
120 IF(ir.EQ.
n+nr-1)
k(ir,2)=
k(
n+nr,2)
122 IF(mju(1).GT.ir) mju(1)=mju(1)-1
123 IF(mju(2).GT.ir) mju(2)=mju(2)-1
134 IF(
ntry.GT.100.AND.ntryr.LE.4)
THEN
138 ELSEIF(
ntry.GT.100)
THEN
139 CALL
luerrm(14,
'(LUSTRF:) caught in infinite loop')
140 IF(mstu(21).GE.1)
RETURN
143 IF(mju(1).EQ.0.AND.mju(2).EQ.0) goto 500
146 IF(mju(jt).EQ.0) goto 490
155 DO 200
i1=
n+1+(jt-1)*(nr-1),
n+nr+(jt-1)*(1-nr),js
156 IF(
k(
i1,2).NE.21.AND.iu.LE.2)
THEN
161 200 pju(iu,
j)=pju(iu,
j)+
p(
i1,
j)
163 210 pju(iu,5)=sqrt(pju(iu,1)**2+pju(iu,2)**2+pju(iu,3)**2)
164 IF(
k(iju(3),2)/100.NE.10*
k(iju(1),2)+
k(iju(2),2).AND.
165 &
k(iju(3),2)/100.NE.10*
k(iju(2),2)+
k(iju(1),2))
THEN
166 CALL
luerrm(12,
'(LUSTRF:) unphysical flavour combination')
167 IF(mstu(21).GE.1)
RETURN
171 t12=(pju(1,1)*pju(2,1)+pju(1,2)*pju(2,2)+pju(1,3)*pju(2,3))/
173 t13=(pju(1,1)*pju(3,1)+pju(1,2)*pju(3,2)+pju(1,3)*pju(3,3))/
175 t23=(pju(2,1)*pju(3,1)+pju(2,2)*pju(3,2)+pju(2,3)*pju(3,3))/
177 t11=sqrt((2./3.)*(1.-t12)*(1.-t13)/(1.-t23))
178 t22=sqrt((2./3.)*(1.-t12)*(1.-t23)/(1.-t13))
179 tsq=sqrt((2.*t11*t22+t12-1.)*(1.+t12))
180 t1f=(tsq-t22*(1.+t12))/(1.-t12**2)
181 t2f=(tsq-t11*(1.+t12))/(1.-t12**2)
183 220 tju(
j)=-(t1f*pju(1,
j)/pju(1,5)+t2f*pju(2,
j)/pju(2,5))
184 tju(4)=sqrt(1.+tju(1)**2+tju(2)**2+tju(3)**2)
186 230 pju(iu,5)=tju(4)*pju(iu,4)-tju(1)*pju(iu,1)-tju(2)*pju(iu,2)-
190 IF(pju(1,5)+pju(2,5).GT.pju(1,4)+pju(2,4))
THEN
210 IF(is.EQ.1) dp(1,
j)=
p(is1,
j)
212 250
IF(is.EQ.ns) dp(2,
j)=-pju(iu,
j)
213 IF(is.EQ.ns) dp(2,4)=sqrt(pju(iu,1)**2+pju(iu,2)**2+pju(iu,3)**2)
214 IF(is.EQ.ns) dp(2,5)=0.
218 IF(dp(3,5)+2.*dhkc+dp(4,5).LE.0.)
THEN
219 dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
220 dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
225 dhks=sqrt(dhkc**2-dp(3,5)*dp(4,5))
226 dhk1=0.5*((dp(4,5)+dhkc)/dhks-1.)
227 dhk2=0.5*((dp(3,5)+dhkc)/dhks-1.)
229 p(in1,5)=sqrt(dp(3,5)+2.*dhkc+dp(4,5))
231 p(in1,
j)=(1.+dhk1)*dp(1,
j)-dhk2*dp(2,
j)
232 260
p(in1+1,
j)=(1.+dhk2)*dp(2,
j)-dhk1*dp(1,
j)
237 IF(
ntry.GT.100.AND.ntryr.LE.4)
THEN
241 ELSEIF(
ntry.GT.100)
THEN
242 CALL
luerrm(14,
'(LUSTRF:) caught in infinite loop')
243 IF(mstu(21).GE.1)
RETURN
247 ie(1)=
k(
n+1+(jt/2)*(
np-1),3)
252 DO 280 in1=
n+nr+2+jq,
n+nr+4*ns-2+jq,4
269 dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
270 dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
271 dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
272 dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
273 dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
274 IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1.
275 IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1.
276 IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1.
277 IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1.
279 dhcx1=dfour(3,1)/dhc12
280 dhcx2=dfour(3,2)/dhc12
281 dhcxx=1d0/sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
282 dhcy1=dfour(4,1)/dhc12
283 dhcy2=dfour(4,2)/dhc12
284 dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
285 dhcyy=1d0/sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
287 dp(3,
j)=dhcxx*(dp(3,
j)-dhcx2*dp(1,
j)-dhcx1*dp(2,
j))
289 310
p(
in(6)+1,
j)=dhcyy*(dp(4,
j)-dhcy2*dp(1,
j)-dhcy1*dp(2,
j)-
294 IF(2*
i-nsav.GE.mstu(4)-mstu(32)-5)
THEN
295 CALL
luerrm(11,
'(LUSTRF:) no more memory left in LUJETS')
296 IF(mstu(21).GE.1)
RETURN
305 330 CALL
lukfdi(kfl(1),0,kfl(3),
k(
i,2))
306 IF(
k(
i,2).EQ.0) goto 270
307 IF(mstj(12).GE.3.AND.irankj.EQ.1.AND.iabs(kfl(1)).LE.10.AND.
308 &iabs(kfl(3)).GT.10)
THEN
309 IF(
rlu(0).GT.parj(19)) goto 330
313 pr(1)=
p(
i,5)**2+(
px(1)+
px(3))**2+(
py(1)+
py(3))**2
314 CALL
luzdis(kfl(1),kfl(3),pr(1),
z)
315 gam(3)=(1.-
z)*(gam(1)+pr(1)/
z)
320 IF(
in(1)+1.EQ.
in(2).AND.
z*
p(
in(1)+2,3)*
p(
in(2)+2,3)*
321 &
p(
in(1),5)**2.GE.pr(1))
THEN
323 p(
in(2)+2,4)=pr(1)/(
p(
in(1)+2,4)*
p(
in(1),5)**2)
327 ELSEIF(
in(1)+1.EQ.
in(2))
THEN
331 IF(
in(2).GT.
n+nr+4*ns) goto 270
332 IF(four(
in(1),
in(2)).LE.1
e-2)
THEN
340 360
IF(
in(1).GT.
n+nr+4*ns.OR.
in(2).GT.
n+nr+4*ns.OR.
341 &
in(1).GT.
in(2)) goto 270
342 IF(
in(1).NE.
in(4).OR.
in(2).NE.
in(5))
THEN
348 dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
349 dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
351 IF(dhc12.LE.1
e-2)
THEN
358 dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
359 dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
360 dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
361 IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1.
362 IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1.
363 IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1.
364 IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1.
365 dhcx1=dfour(3,1)/dhc12
366 dhcx2=dfour(3,2)/dhc12
367 dhcxx=1d0/sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
368 dhcy1=dfour(4,1)/dhc12
369 dhcy2=dfour(4,2)/dhc12
370 dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
371 dhcyy=1d0/sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
373 dp(3,
j)=dhcxx*(dp(3,
j)-dhcx2*dp(1,
j)-dhcx1*dp(2,
j))
375 380
p(
in(3)+1,
j)=dhcyy*(dp(4,
j)-dhcy2*dp(1,
j)-dhcy1*dp(2,
j)-
378 pxp=-(
px(3)*four(
in(6),
in(3))+
py(3)*four(
in(6)+1,
in(3)))
380 IF(abs(pxp**2+
pyp**2-
px(3)**2-
py(3)**2).LT.0.01)
THEN
391 DO 390 in1=
in(4),
in(1)-4,4
393 DO 400 in2=
in(5),
in(2)-4,4
396 dhm(2)=2.*four(
i,
in(1))
397 dhm(3)=2.*four(
i,
in(2))
398 dhm(4)=2.*four(
in(1),
in(2))
401 DO 410 in2=
in(1)+1,
in(2),4
402 DO 410 in1=
in(1),in2-1,4
404 dhg(1)=dhg(1)+
p(in1+2,1)*
p(in2+2,1)*dhc
405 IF(in1.EQ.
in(1)) dhg(2)=dhg(2)-
p(in2+2,1)*dhc
406 IF(in2.EQ.
in(2)) dhg(3)=dhg(3)+
p(in1+2,1)*dhc
407 410
IF(in1.EQ.
in(1).AND.in2.EQ.
in(2)) dhg(4)=dhg(4)-dhc
410 dhs1=dhm(3)*dhg(4)-dhm(4)*dhg(3)
411 IF(abs(dhs1).LT.1
e-4) goto 270
412 dhs2=dhm(4)*(gam(3)-dhg(1))-dhm(2)*dhg(3)-dhg(4)*
413 &(
p(
i,5)**2-dhm(1))+dhg(2)*dhm(3)
414 dhs3=dhm(2)*(gam(3)-dhg(1))-dhg(2)*(
p(
i,5)**2-dhm(1))
415 p(
in(2)+2,4)=0.5*(sqrt(
max(0d0,dhs2**2-4.*dhs1*dhs3))/abs(dhs1)-
417 IF(dhm(2)+dhm(4)*
p(
in(2)+2,4).LE.0.) goto 270
418 p(
in(1)+2,4)=(
p(
i,5)**2-dhm(1)-dhm(3)*
p(
in(2)+2,4))/
419 &(dhm(2)+dhm(4)*
p(
in(2)+2,4))
422 IF(
p(
in(2)+2,4).GT.
p(
in(2)+2,3))
THEN
426 IF(
in(2).GT.
n+nr+4*ns) goto 270
427 IF(four(
in(1),
in(2)).LE.1
e-2)
THEN
433 ELSEIF(
p(
in(1)+2,4).GT.
p(
in(1)+2,3))
THEN
443 430 pju(iu+3,
j)=pju(iu+3,
j)+
p(
i,
j)
444 IF(
p(
i,4).LE.0.) goto 270
445 pju(iu+3,5)=tju(4)*pju(iu+3,4)-tju(1)*pju(iu+3,1)-
446 &tju(2)*pju(iu+3,2)-tju(3)*pju(iu+3,3)
447 IF(pju(iu+3,5).LT.pju(iu,5))
THEN
452 IF(
in(3).NE.
in(6))
THEN
459 p(
in(jq)+2,3)=
p(
in(jq)+2,3)-
p(
in(jq)+2,4)
460 450
p(
in(jq)+2,1)=
p(
in(jq)+2,1)-(3-2*jq)*
p(
in(jq)+2,4)
465 IF(iabs(kfl(1)).GT.10) goto 270
469 460 pju(iu+3,
j)=pju(iu+3,
j)-
p(
i+1,
j)
474 kfjs(jt)=
k(
k(mju(jt+2),3),2)
475 kfls=2*int(
rlu(0)+3.*parj(4)/(1.+3.*parj(4)))+1
476 IF(kfjh(1).EQ.kfjh(2)) kfls=3
477 IF(ista.NE.
i) kfjs(jt)=isign(1000*
max(iabs(kfjh(1)),
478 &iabs(kfjh(2)))+100*min(iabs(kfjh(1)),iabs(kfjh(2)))+
481 pjs(jt,
j)=pju(1,
j)+pju(2,
j)+
p(mju(jt),
j)
482 480 pjs(jt+2,
j)=pju(4,
j)+pju(5,
j)
483 pjs(jt,5)=sqrt(
max(0.,pjs(jt,4)**2-pjs(jt,1)**2-pjs(jt,2)**2-
488 500
IF(mju(1).NE.0.AND.mju(2).NE.0)
THEN
491 ELSEIF(mju(1).NE.0)
THEN
494 ELSEIF(mju(2).NE.0)
THEN
497 ELSEIF(iabs(
k(
n+1,2)).NE.21)
THEN
504 p(
n+nr+is,1)=0.5*four(
n+is,
n+is+1-nr*(is/nr))
505 510 w2sum=w2sum+
p(
n+nr+is,1)
509 w2sum=w2sum-
p(
n+nr+nb,1)
510 IF(w2sum.GT.w2ran.AND.nb.LT.nr) goto 520
515 is1=
n+is+nb-1-nr*((is+nb-2)/nr)
516 is2=
n+is+nb-nr*((is+nb-1)/nr)
519 IF(iabs(
k(is1,2)).EQ.21) dp(1,
j)=0.5*dp(1,
j)
520 IF(is1.EQ.mju(1)) dp(1,
j)=pjs(1,
j)-pjs(3,
j)
522 IF(iabs(
k(is2,2)).EQ.21) dp(2,
j)=0.5*dp(2,
j)
523 530
IF(is2.EQ.mju(2)) dp(2,
j)=pjs(2,
j)-pjs(4,
j)
527 IF(dp(3,5)+2.*dhkc+dp(4,5).LE.0.)
THEN
530 dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2+dp(1,5)**2)
531 dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2+dp(2,5)**2)
534 dhks=sqrt(dhkc**2-dp(3,5)*dp(4,5))
535 dhk1=0.5*((dp(4,5)+dhkc)/dhks-1.)
536 dhk2=0.5*((dp(3,5)+dhkc)/dhks-1.)
538 p(in1,5)=sqrt(dp(3,5)+2.*dhkc+dp(4,5))
540 p(in1,
j)=(1.+dhk1)*dp(1,
j)-dhk2*dp(2,
j)
541 540
p(in1+1,
j)=(1.+dhk2)*dp(2,
j)-dhk1*dp(1,
j)
546 IF(
ntry.GT.100.AND.ntryr.LE.4)
THEN
550 ELSEIF(
ntry.GT.100)
THEN
551 CALL
luerrm(14,
'(LUSTRF:) caught in infinite loop')
552 IF(mstu(21).GE.1)
RETURN
561 IF(mju(jt).NE.0) irank(jt)=njs(jt)
562 IF(ns.GT.nr) irank(jt)=1
563 ie(jt)=
k(
n+1+(jt/2)*(
np-1),3)
564 in(3*jt+1)=
n+nr+1+4*(jt/2)*(ns-1)
565 in(3*jt+2)=
in(3*jt+1)+1
566 in(3*jt+3)=
n+nr+4*ns+2*jt-1
567 DO 570 in1=
n+nr+2+jt,
n+nr+4*ns-2+jt,4
576 IF(ns.EQ.1.AND.mju(1)+mju(2).EQ.0) CALL
luptdi(0,
px(1),
py(1))
581 IF(mju(jt).NE.0) kfl(jt)=kfjs(jt)
588 kfl(3)=int(1.+(2.+parj(2))*
rlu(0))*(-1)**int(
rlu(0)+0.5)
589 CALL
lukfdi(kfl(3),0,kfl(1),kdump)
591 IF(iabs(kfl(1)).GT.10.AND.
rlu(0).GT.0.5)
THEN
592 kfl(2)=-(kfl(1)+isign(10000,kfl(1)))
593 ELSEIF(iabs(kfl(1)).GT.10)
THEN
594 kfl(1)=-(kfl(2)+isign(10000,kfl(2)))
599 pr3=min(25.,0.1*
p(
n+nr+1,5)**2)
600 590 CALL
luzdis(kfl(1),kfl(2),pr3,
z)
601 zr=pr3/(
z*
p(
n+nr+1,5)**2)
602 IF(zr.GE.1.) goto 590
607 in1=
n+nr+3+4*(jt/2)*(ns-1)
610 p(in1,3)=(2-jt)*(1.-
z)+(jt-1)*
z
613 600
p(in1+1,3)=(2-jt)*(1.-zr)+(jt-1)*zr
618 IF(jt.EQ.1.OR.ns.EQ.nr-1)
THEN
626 dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
627 dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
628 dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
629 dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
630 dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
631 IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1.
632 IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1.
633 IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1.
634 IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1.
636 dhcx1=dfour(3,1)/dhc12
637 dhcx2=dfour(3,2)/dhc12
638 dhcxx=1d0/sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
639 dhcy1=dfour(4,1)/dhc12
640 dhcy2=dfour(4,2)/dhc12
641 dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
642 dhcyy=1d0/sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
644 dp(3,
j)=dhcxx*(dp(3,
j)-dhcx2*dp(1,
j)-dhcx1*dp(2,
j))
646 620
p(in3+1,
j)=dhcyy*(dp(4,
j)-dhcy2*dp(1,
j)-dhcy1*dp(2,
j)-
651 630
p(in3+3,
j)=
p(in3+1,
j)
656 IF(mju(1)+mju(2).GT.0)
THEN
658 IF(njs(jt).EQ.0) goto 660
660 650
p(
n+nrs,
j)=
p(
n+nrs,
j)-pjs(jt+2,
j)
666 IF(2*
i-nsav.GE.mstu(4)-mstu(32)-5)
THEN
667 CALL
luerrm(11,
'(LUSTRF:) no more memory left in LUJETS')
668 IF(mstu(21).GE.1)
RETURN
671 IF(iabs(kfl(3-jt)).GT.10) jt=3-jt
674 irank(jt)=irank(jt)+1
681 680 CALL
lukfdi(kfl(jt),0,kfl(3),
k(
i,2))
682 IF(
k(
i,2).EQ.0) goto 550
683 IF(mstj(12).GE.3.AND.irank(jt).EQ.1.AND.iabs(kfl(jt)).LE.10.AND.
684 &iabs(kfl(3)).GT.10)
THEN
685 IF(
rlu(0).GT.parj(19)) goto 680
689 pr(jt)=
p(
i,5)**2+(
px(jt)+
px(3))**2+(
py(jt)+
py(3))**2
694 wmin=parj(32+mstj(11))+pmq(1)+pmq(2)+parj(36)*pmq(3)
695 IF(iabs(kfl(jt)).GT.10.AND.iabs(kfl(3)).GT.10) wmin=
696 &wmin-0.5*parj(36)*pmq(3)
697 wrem2=four(
n+nrs,
n+nrs)
698 IF(wrem2.LT.0.10) goto 550
699 IF(wrem2.LT.
max(wmin*(1.+(2.*
rlu(0)-1.)*parj(37)),
700 &parj(32)+pmq(1)+pmq(2))**2) goto 810
703 CALL
luzdis(kfl(jt),kfl(3),pr(jt),
z)
706 IF(
max(mod(kfl1a,10),mod(kfl1a/1000,10),mod(kfl2a,10),
707 &mod(kfl2a/1000,10)).GE.4)
THEN
708 pr(jr)=(pmq(jr)+pmq(3))**2+(
px(jr)-
px(3))**2+(
py(jr)-
py(3))**2
709 pw12=sqrt(
max(0.,(wrem2-pr(1)-pr(2))**2-4.*pr(1)*pr(2)))
710 z=(wrem2+pr(jt)-pr(jr)+pw12*(2.*
z-1.))/(2.*wrem2)
711 pr(jr)=(pmq(jr)+parj(32+mstj(11)))**2+(
px(jr)-
px(3))**2+
713 IF((1.-
z)*(wrem2-pr(jt)/
z).LT.pr(jr)) goto 810
715 gam(3)=(1.-
z)*(gam(jt)+pr(jt)/
z)
720 IF(
in(1)+1.EQ.
in(2).AND.
z*
p(
in(1)+2,3)*
p(
in(2)+2,3)*
721 &
p(
in(1),5)**2.GE.pr(jt))
THEN
723 p(
in(jr)+2,4)=pr(jt)/(
p(
in(jt)+2,4)*
p(
in(1),5)**2)
727 ELSEIF(
in(1)+1.EQ.
in(2))
THEN
728 p(
in(jr)+2,4)=
p(
in(jr)+2,3)
731 IF(js*
in(jr).GT.js*
in(4*jr)) goto 550
732 IF(four(
in(1),
in(2)).LE.1
e-2)
THEN
733 p(
in(jt)+2,4)=
p(
in(jt)+2,3)
740 710
IF(js*
in(1).GT.js*
in(3*jr+1).OR.js*
in(2).GT.js*
in(3*jr+2).OR.
741 &
in(1).GT.
in(2)) goto 550
742 IF(
in(1).NE.
in(3*jt+1).OR.
in(2).NE.
in(3*jt+2))
THEN
748 dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
749 dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
751 IF(dhc12.LE.1
e-2)
THEN
752 p(
in(jt)+2,4)=
p(
in(jt)+2,3)
758 dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
759 dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
760 dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
761 IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1.
762 IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1.
763 IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1.
764 IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1.
765 dhcx1=dfour(3,1)/dhc12
766 dhcx2=dfour(3,2)/dhc12
767 dhcxx=1d0/sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
768 dhcy1=dfour(4,1)/dhc12
769 dhcy2=dfour(4,2)/dhc12
770 dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
771 dhcyy=1d0/sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
773 dp(3,
j)=dhcxx*(dp(3,
j)-dhcx2*dp(1,
j)-dhcx1*dp(2,
j))
775 730
p(
in(3)+1,
j)=dhcyy*(dp(4,
j)-dhcy2*dp(1,
j)-dhcy1*dp(2,
j)-
778 pxp=-(
px(3)*four(
in(3*jt+3),
in(3))+
py(3)*
779 & four(
in(3*jt+3)+1,
in(3)))
781 & four(
in(3*jt+3)+1,
in(3)+1))
782 IF(abs(pxp**2+
pyp**2-
px(3)**2-
py(3)**2).LT.0.01)
THEN
793 DO 740 in1=
in(3*jt+1),
in(1)-4*js,4*js
795 DO 750 in2=
in(3*jt+2),
in(2)-4*js,4*js
798 dhm(2)=2.*four(
i,
in(1))
799 dhm(3)=2.*four(
i,
in(2))
800 dhm(4)=2.*four(
in(1),
in(2))
803 DO 760 in2=
in(1)+1,
in(2),4
804 DO 760 in1=
in(1),in2-1,4
806 dhg(1)=dhg(1)+
p(in1+2,jt)*
p(in2+2,jt)*dhc
807 IF(in1.EQ.
in(1)) dhg(2)=dhg(2)-js*
p(in2+2,jt)*dhc
808 IF(in2.EQ.
in(2)) dhg(3)=dhg(3)+js*
p(in1+2,jt)*dhc
809 760
IF(in1.EQ.
in(1).AND.in2.EQ.
in(2)) dhg(4)=dhg(4)-dhc
812 dhs1=dhm(jr+1)*dhg(4)-dhm(4)*dhg(jr+1)
813 IF(abs(dhs1).LT.1
e-4) goto 550
814 dhs2=dhm(4)*(gam(3)-dhg(1))-dhm(jt+1)*dhg(jr+1)-dhg(4)*
815 &(
p(
i,5)**2-dhm(1))+dhg(jt+1)*dhm(jr+1)
816 dhs3=dhm(jt+1)*(gam(3)-dhg(1))-dhg(jt+1)*(
p(
i,5)**2-dhm(1))
817 p(
in(jr)+2,4)=0.5*(sqrt(
max(0d0,dhs2**2-4.*dhs1*dhs3))/abs(dhs1)-
819 IF(dhm(jt+1)+dhm(4)*
p(
in(jr)+2,4).LE.0.) goto 550
820 p(
in(jt)+2,4)=(
p(
i,5)**2-dhm(1)-dhm(jr+1)*
p(
in(jr)+2,4))/
821 &(dhm(jt+1)+dhm(4)*
p(
in(jr)+2,4))
824 IF(
p(
in(jr)+2,4).GT.
p(
in(jr)+2,3))
THEN
825 p(
in(jr)+2,4)=
p(
in(jr)+2,3)
828 IF(js*
in(jr).GT.js*
in(4*jr)) goto 550
829 IF(four(
in(1),
in(2)).LE.1
e-2)
THEN
830 p(
in(jt)+2,4)=
p(
in(jt)+2,3)
835 ELSEIF(
p(
in(jt)+2,4).GT.
p(
in(jt)+2,3))
THEN
836 p(
in(jt)+2,4)=
p(
in(jt)+2,3)
846 IF(
p(
i,4).LE.0.) goto 550
852 IF(
in(3).NE.
in(3*jt+3))
THEN
859 p(
in(jq)+2,3)=
p(
in(jq)+2,3)-
p(
in(jq)+2,4)
860 800
p(
in(jq)+2,jt)=
p(
in(jq)+2,jt)-js*(3-2*jq)*
p(
in(jq)+2,4)
869 CALL
lukfdi(kfl(jr),-kfl(3),kfldmp,
k(
i,2))
870 IF(
k(
i,2).EQ.0) goto 550
872 pr(jr)=
p(
i,5)**2+(
px(jr)-
px(3))**2+(
py(jr)-
py(3))**2
876 IF(
p(
in(4)+2,3)*
p(
in(5)+2,3)*four(
in(4),
in(5)).LT.
p(
in(7),3)*
877 &
p(
in(8),3)*four(
in(7),
in(8))) jq=2
878 dhc12=four(
in(3*jq+1),
in(3*jq+2))
879 dhr1=four(
n+nrs,
in(3*jq+2))/dhc12
880 dhr2=four(
n+nrs,
in(3*jq+1))/dhc12
881 IF(
in(4).NE.
in(7).OR.
in(5).NE.
in(8))
THEN
882 px(3-jq)=-four(
n+nrs,
in(3*jq+3))-
px(jq)
883 py(3-jq)=-four(
n+nrs,
in(3*jq+3)+1)-
py(jq)
884 pr(3-jq)=
p(
i+(jt+jq-3)**2-1,5)**2+(
px(3-jq)+(2*jq-3)*js*
885 &
px(3))**2+(
py(3-jq)+(2*jq-3)*js*
py(3))**2
889 wrem2=wrem2+(
px(1)+
px(2))**2+(
py(1)+
py(2))**2
890 fd=(sqrt(pr(1))+sqrt(pr(2)))/sqrt(wrem2)
891 IF(mju(1)+mju(2).NE.0.AND.
i.EQ.isav+2.AND.
fd.GE.1.) goto 180
892 IF(
fd.GE.1.) goto 550
893 fa=wrem2+pr(jt)-pr(jr)
894 IF(mstj(11).EQ.2)
prev=0.5*
fd**parj(37+mstj(11))
896 &parj(37+mstj(11))*(pr(1)+pr(2))**2))
900 IF(
max(mod(kfl1a,10),mod(kfl1a/1000,10),mod(kfl2a,10),
901 &mod(kfl2a/1000,10)).GE.6) fb=
sign(sqrt(
max(0.,fa**2-
902 &4.*wrem2*pr(jt))),float(js))
905 &
p(
in(3*jq+3)+1,
j)+0.5*(dhr1*(fa+fb)*
p(
in(3*jq+1),
j)+
906 &dhr2*(fa-fb)*
p(
in(3*jq+2),
j))/wrem2
911 DO 830
i=nsav+1,nsav+
np
914 IF(mstu(16).NE.2)
THEN
932 840
v(nsav,
j)=
v(ip,
j)
933 p(nsav,5)=sqrt(
max(0d0,dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2))
947 DO 880
i=
n+1,2*
n-nsav
948 IF(
k(
i,3).NE.ie(1)) goto 880
953 IF(mstu(16).NE.2)
k(
i1,3)=nsav
955 DO 900
i=2*
n-nsav,
n+1,-1
956 IF(
k(
i,3).EQ.ie(1)) goto 900
961 IF(mstu(16).NE.2)
k(
i1,3)=nsav
965 CALL ludbrb(nsav+1,
n,0.,0.,dps(1)/dps(4),dps(2)/dps(4),