10 IMPLICIT DOUBLE PRECISION(a-
h, o-
z)
14 parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
15 &kexcit=4000000,kdimen=5000000)
18 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
19 common/
pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
22 common/
pyssmt/zmix(4,4),umix(2,2),vmix(2,2),smz(4),smw(2),
23 &sfmix(16,4),zmixi(4,4),umixi(2,2),vmixi(2,2)
28 DOUBLE PRECISION xm(5)
29 COMPLEX*16 olpp,orpp,qll,qlr,qrr,qrl,glij,grij,propz
30 COMPLEX*16 qlls,qrrs,qlrs,qrls,qllu,qrru,qlrt,qrlt
31 COMPLEX*16 zmixc(4,4),umixc(2,2),vmixc(2,2)
32 DOUBLE PRECISION s12min,s12max,yjaco1,s23ave,s23df1,s23df2
33 DOUBLE PRECISION d1,d2,d3,p1,p2,p3,cthe1,sthe1,cthe3,sthe3
34 DOUBLE PRECISION cphi1,sphi1
35 DOUBLE PRECISION s23del,
eps
36 DOUBLE PRECISION golden,ax,bx,cx,tol,
xmin,
r,
c
37 parameter(
r=0.61803399d0,
c=1d0-
r,tol=1d-3)
38 DOUBLE PRECISION f1,
f2,
x0,x1,x2,x3
40 DATA inoid/22,23,25,35/
51 s12min=(xm(1)+xm(2))**2
52 s12max=(xm(5)-xm(3))**2
68 IF(iabs(
k(
id,2)).LT.ksusy1.OR.iabs(
k(
id,2)).GE.3000000)
THEN
72 IF(mod(
k(
n+1,2),ksusy1).EQ.inoid(
i1)) izid1=
i1
73 IF(mod(
k(
id,2),ksusy1).EQ.inoid(
i1)) izid2=
i1
75 IF(mod(
k(
n+1,2),ksusy1).EQ.24) iwid1=1
76 IF(mod(
k(
n+1,2),ksusy1).EQ.37) iwid1=2
77 IF(mod(
k(
id,2),ksusy1).EQ.24) iwid2=1
78 IF(mod(
k(
id,2),ksusy1).EQ.37) iwid2=2
85 IF(
max(abs(ia),abs(ja)).EQ.6)
THEN
87 ELSEIF(izid1*izid2.NE.0)
THEN
89 gmmz=pmas(23,1)*pmas(23,2)
91 zmixc(izid1,
i)=dcmplx(zmix(izid1,
i),zmixi(izid1,
i))
92 zmixc(izid2,
i)=dcmplx(zmix(izid2,
i),zmixi(izid2,
i))
94 olpp=(zmixc(izid1,3)*dconjg(zmixc(izid2,3))-
95 & zmixc(izid1,4)*dconjg(zmixc(izid2,4)))/2d0
97 xll2=pmas(
pycomp(ksusy1+iabs(ia)),1)**2
99 xrr2=pmas(
pycomp(ksusy2+iabs(ia)),1)**2
101 glij=(t3i*zmixc(izid1,2)-tanw*(t3i-
ei)*zmixc(izid1,1))*
102 & dconjg(t3i*zmixc(izid2,2)-tanw*(t3i-
ei)*zmixc(izid2,1))
103 grij=zmixc(izid1,1)*dconjg(zmixc(izid2,1))*(
ei*tanw)**2
104 xm1m2=smz(izid1)*smz(izid2)
105 qlls=dcmplx((t3i-
ei*xw)/xw1)*olpp
107 qlrs=-dcmplx((t3i-
ei*xw)/xw1)*orpp
109 qrls=-dcmplx((
ei*xw)/xw1)*olpp
111 qrrs=dcmplx((
ei*xw)/xw1)*orpp
113 ELSEIF(izid1*iwid2.NE.0.OR.izid2*iwid1.NE.0)
THEN
115 xm1m2=smz(izid1)*smw(iwid2)
119 xm1m2=smz(izid2)*smw(iwid1)
124 gmmz=pmas(24,1)*pmas(24,2)
126 vmixc(izid1,
i)=dcmplx(vmix(izid1,
i),vmixi(izid1,
i))
127 umixc(izid1,
i)=dcmplx(umix(izid1,
i),umixi(izid1,
i))
130 zmixc(izid2,
i)=dcmplx(zmix(izid2,
i),zmixi(izid2,
i))
132 qlls=(dconjg(zmixc(izid2,2))*vmixc(izid1,1)-
133 & dconjg(zmixc(izid2,4))*vmixc(izid1,2)*rt2i)
134 qlrs=(zmixc(izid2,2)*dconjg(umixc(izid1,1))+
135 & zmixc(izid2,3)*dconjg(umixc(izid1,2))*rt2i)
136 ej=kchg(iabs(ja),1)/3d0
137 t3j=
sign(1d0,ej+1d-6)/2d0
144 xlr2 = pmas(
pycomp(ksusy1+iabs(ja)),1)**2
145 xll2 = pmas(
pycomp(ksusy1+iabs(ia)),1)**2
146 IF(mod(ia,2).EQ.0)
THEN
147 qllu=vmixc(izid1,1)*dconjg(zmixc(izid2,1)*(
ei-t3i)*
148 & tanw+zmixc(izid2,2)*t3i)
149 qlrt=-dconjg(umixc(izid1,1))*(
150 & zmixc(izid2,1)*(ej-t3j)*tanw+zmixc(izid2,2)*t3j)
152 qllu=vmixc(izid1,1)*dconjg(zmixc(izid2,1)*(ej-t3j)*
153 & tanw+zmixc(izid2,2)*t3j)
154 qlrt=-dconjg(umixc(izid1,1))*(
155 & zmixc(izid2,1)*(
ei-t3i)*tanw+zmixc(izid2,2)*t3i)
157 ELSEIF(iwid1*iwid2.NE.0)
THEN
160 xm1m2=smw(iwid1)*smw(iwid2)
162 gmmz=pmas(23,1)*pmas(23,2)
164 vmixc(izid1,
i)=dcmplx(vmix(izid1,
i),vmixi(izid1,
i))
165 umixc(izid1,
i)=dcmplx(umix(izid1,
i),umixi(izid1,
i))
166 vmixc(izid2,
i)=dcmplx(vmix(izid2,
i),vmixi(izid2,
i))
167 umixc(izid2,
i)=dcmplx(umix(izid2,
i),umixi(izid2,
i))
169 olpp=-vmixc(izid2,1)*dconjg(vmixc(izid1,1))-
170 & vmixc(izid2,2)*dconjg(vmixc(izid1,2))/2d0
171 orpp=-umixc(izid1,1)*dconjg(umixc(izid2,1))-
172 & umixc(izid1,2)*dconjg(umixc(izid2,2))/2d0
173 qrls=-dcmplx(
ei/xw1)*orpp
174 qlls=dcmplx((t3i-xw*
ei)/xw/xw1)*orpp
175 qrrs=-dcmplx(
ei/xw1)*olpp
176 qlrs=dcmplx((t3i-xw*
ei)/xw/xw1)*olpp
177 IF(mod(ia,2).EQ.0)
THEN
178 xlr2=pmas(
pycomp(ksusy1+iabs(ia)-1),1)**2
179 qlrt=-umixc(izid2,1)*dconjg(umixc(izid1,1))*dcmplx(t3i/xw)
181 xlr2=pmas(
pycomp(ksusy1+iabs(ia)+1),1)**2
182 qlrt=-vmixc(izid2,1)*dconjg(vmixc(izid1,1))*dcmplx(t3i/xw)
184 ELSEIF(mod(
k(
n+1,2),ksusy1).EQ.21.OR.mod(
k(
id,2),ksusy1).EQ.21)
194 s12=s12min+yjaco1*(kt-1)/99
195 s23ave=xm(2)**2+xm(3)**2-(s12+xm(2)**2-xm(1)**2)
196 & *(s12+xm(3)**2-xm(5)**2)/(2d0*s12)
197 s23df1=(s12-xm(2)**2-xm(1)**2)**2
198 & -(2d0*xm(1)*xm(2))**2
199 s23df2=(s12-xm(3)**2-xm(5)**2)**2
200 & -(2d0*xm(3)*xm(5))**2
203 s23del=sqrt(
max(0d0,s23df1*s23df2))/(2d0*s12)
210 s23=s23min+yjaco2*(ks-1)/99
213 wu2 = (uh-zm12)*(uh-zm22)
214 wt2 = (th-zm12)*(th-zm22)
216 propz2 = (sh-sqmz)**2 + gmmz**2
217 propz=dcmplx(sh-sqmz,-gmmz)/dcmplx(propz2)
218 qll=qlls*propz+qllu/dcmplx(uh-xll2)
219 qlr=qlrs*propz+qlrt/dcmplx(th-xlr2)
220 qrl=qrls*propz+qrlt/dcmplx(th-xrl2)
221 qrr=qrrs*propz+qrru/dcmplx(uh-xrr2)
222 wt0=-((abs(qll)**2+abs(qrr)**2)*wu2+
223 & (abs(qrl)**2+abs(qlr)**2)*wt2+
224 & 2d0*dble(qlr*dconjg(qll)+qrl*dconjg(qrr))*ws2)
225 IF(wt0.GT.wtmax) wtmax=wt0
235 bx=s12min+0.5d0*yjaco1
238 IF(abs(cx-bx).GT.abs(bx-ax))
THEN
247 s23df1=(x1-xm(2)**2-xm(1)**2)**2
248 &-(2d0*xm(1)*xm(2))**2
249 s23df2=(x1-xm(3)**2-xm(5)**2)**2
250 &-(2d0*xm(3)*xm(5))**2
253 s23del=sqrt(
max(0d0,s23df1*s23df2))/(2d0*x1)
255 s23df1=(x2-xm(2)**2-xm(1)**2)**2
256 &-(2d0*xm(1)*xm(2))**2
257 s23df2=(x2-xm(3)**2-xm(5)**2)**2
258 &-(2d0*xm(3)*xm(5))**2
261 s23del=sqrt(
max(0d0,s23df1*s23df2))/(2d0*x2)
264 170
IF(abs(x3-
x0).GT.tol*(abs(x1)+abs(x2)))
THEN
271 s23df1=(x2-xm(2)**2-xm(1)**2)**2
272 & -(2d0*xm(1)*xm(2))**2
273 s23df2=(x2-xm(3)**2-xm(5)**2)**2
274 & -(2d0*xm(3)*xm(5))**2
277 s23del=sqrt(
max(0d0,s23df1*s23df2))/(2d0*x2)
284 s23df1=(x1-xm(2)**2-xm(1)**2)**2
285 & -(2d0*xm(1)*xm(2))**2
286 s23df2=(x1-xm(3)**2-xm(5)**2)**2
287 & -(2d0*xm(3)*xm(5))**2
290 s23del=sqrt(
max(0d0,s23df1*s23df2))/(2d0*x1)
305 180 s12=s12min+
pyr(0)*yjaco1
308 s23ave=xm(2)**2+xm(3)**2-(s12+xm(2)**2-xm(1)**2)
309 &*(s12+xm(3)**2-xm(5)**2)/(2d0*s12)
310 s23df1=(s12-xm(2)**2-xm(1)**2)**2
311 &-(2d0*xm(1)*xm(2))**2
312 s23df2=(s12-xm(3)**2-xm(5)**2)**2
313 &-(2d0*xm(3)*xm(5))**2
316 s23del=sqrt(
max(0d0,s23df1*s23df2))/(2d0*s12)
321 s23=s23min+
pyr(0)*yjaco2
325 WRITE(mstu(11),*)
' IKNT > 100 IN PYTBDY '
328 IF(yjaco2.LT.
pyr(0)*golden) goto 180
330 IF(iskip.EQ.0) goto 190
336 wu2 = (uh-zm12)*(uh-zm22)
337 wt2 = (th-zm12)*(th-zm22)
339 propz2 = (sh-sqmz)**2 + gmmz**2
340 propz=dcmplx(sh-sqmz,-gmmz)/dcmplx(propz2)
342 qll=qlls*propz+qllu/dcmplx(uh-xll2)
343 qlr=qlrs*propz+qlrt/dcmplx(th-xlr2)
344 qrl=qrls*propz+qrlt/dcmplx(th-xrl2)
345 qrr=qrrs*propz+qrru/dcmplx(uh-xrr2)
352 wt=-((abs(qll)**2+abs(qrr)**2)*wu2+
353 &(abs(qrl)**2+abs(qlr)**2)*wt2+
354 &2d0*dble(qlr*dconjg(qll)+qrl*dconjg(qrr))*ws2)
356 IF(wt.LT.
pyr(0)*wtmax) goto 180
357 IF(wt.GT.wtmax)
print*,
' WT > WTMAX ',wt,wtmax
359 190 d3=(xm(5)**2+xm(3)**2-s12)/(2d0*xm(5))
360 d1=(xm(5)**2+xm(1)**2-s23)/(2d0*xm(5))
362 p1=sqrt(
d1*
d1-xm(1)**2)
363 p2=sqrt(d2*d2-xm(2)**2)
364 p3=sqrt(d3*d3-xm(3)**2)
366 ang1=2d0*
pyr(0)*paru(1)
370 IF(arg.LT.0d0.AND.arg.GT.-1d-3) arg=0d0
372 p(
n+1,1)=p1*sthe1*cphi1
373 p(
n+1,2)=p1*sthe1*sphi1
378 ang3=2d0*
pyr(0)*paru(1)
381 cthe3=(p2**2-p1**2-p3**2)/2d0/p1/p3
383 IF(arg.LT.0d0.AND.arg.GT.-1d-3) arg=0d0
385 p(
n+3,1)=-p3*sthe3*cphi3*cthe1*cphi1
386 &+p3*sthe3*sphi3*sphi1
387 &+p3*cthe3*sthe1*cphi1
388 p(
n+3,2)=-p3*sthe3*cphi3*cthe1*sphi1
389 &-p3*sthe3*sphi3*cphi1
390 &+p3*cthe3*sthe1*sphi1
391 p(
n+3,3)=p3*sthe3*cphi3*sthe1