7 IMPLICIT DOUBLE PRECISION(
d)
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 common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200)
16 common/pyhipars/mstp(200),parp(200),msti(200),pari(200)
18 common/pyhiint1/mint(400),vint(400)
20 common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
22 common/pyhiint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
24 dimension kfls(4),
is(2),xs(2),zs(2),q2s(2),tevs(2),robo(5),
25 &xfs(2,-6:6),xfa(-6:6),xfb(-6:6),xfn(-6:6),wtap(-6:6),wtsf(-6:6),
26 &the2(2),alam(2),dq2(3),dpc(3),dpd(4),dpb(4)
33 IF(iset(isub).EQ.1)
THEN
35 ELSEIF(iset(isub).EQ.3.OR.iset(isub).EQ.4)
THEN
37 IF(isub.EQ.8.OR.isub.EQ.76.OR.isub.EQ.77) q2e=pmas(24,1)**2
39 tmax=
log(parp(67)*parp(63)*q2e/parp(61)**2)
40 IF(parp(67)*q2e.LT.
max(parp(62)**2,2.*parp(61)**2).OR.
44 xe0=2.*parp(65)/vint(1)
59 110 xfs(jt,kfl)=xsfx(jt,kfl)
61 IF(iset(isub).EQ.3.OR.iset(isub).EQ.4) dsh=vint(26)*vint(2)
66 IF(
n.GT.
ns+1.AND.q2s(2).GT.q2s(1)) jt=2
70 130 xfb(kfl)=xfs(jt,kfl)
73 xe=
max(xe0,
xb*(1./(1.-parp(66))-1.))
74 IF(
xb+xe.GE.0.999)
THEN
80 IF(mstp(62).LE.1)
THEN
81 q2b=0.5*(1./zs(jt)+1.)*q2s(jt)+0.5*(1./zs(jt)-1.)*(q2s(3-jt)-
82 & sngl(dsh)+
sqrt((sngl(dsh)+q2s(1)+q2s(2))**2+8.*q2s(1)*q2s(2)*
83 & zs(jt)/(1.-zs(jt))))
84 tevb=
log(parp(63)*q2b/alam(jt)**2)
89 alsdum=
ulalps(parp(63)*q2b)
90 tevb=tevb+2.*
log(alam(jt)/paru(117))
93 b0=(33.-2.*mstu(118))/6.
101 DO 150 kfl=-mstp(54),mstp(54)
102 IF(kfl.EQ.0) wtap(kfl)=6.*
log((1.-
xb)/xe)
103 150
IF(kfl.NE.0) wtap(kfl)=wtapq
105 wtap(0)=0.5*
xb*(1./(
xb+xe)-1.)
106 wtap(kflb)=8.*
log((1.-
xb)*(
xb+xe)/xe)/3.
109 IF(kflb.NE.21) xfbo=xfb(kflb)
110 IF(kflb.EQ.21) xfbo=xfb(0)
115 WRITE(mstu(11),1001) kflb,xfb(kflb)
119 DO 170 kfl=-mstp(54),mstp(54)
120 wtsf(kfl)=xfb(kfl)/xfbo
121 170 wtsum=wtsum+wtap(kfl)*wtsf(kfl)
122 wtsum=
max(0.0001,wtsum)
125 180
IF(mstp(64).LE.0)
THEN
126 tevb=tevb+
log(
rlu(0))*paru(2)/(paru(111)*wtsum)
127 ELSEIF(mstp(64).EQ.1)
THEN
132 190 q2ref=alam(jt)**2*
exp(tevb)
136 IF(q2b.LT.parp(62)**2)
THEN
142 wtran=wtran-wtap(kfla)*wtsf(kfla)
143 IF(kfla.LT.mstp(54).AND.wtran.GT.0.) goto 200
144 IF(kfla.EQ.0) kfla=21
147 IF(kflb.EQ.21.AND.kfla.EQ.21)
THEN
150 ELSEIF(kflb.EQ.21)
THEN
152 wtz=0.5*(1.+(1.-
z)**2)*
sqrt(
z)
153 ELSEIF(kfla.EQ.21)
THEN
162 IF(mstp(65).GE.1)
THEN
164 IF(kflb.NE.21) rsoft=8./3.
165 z=
z*(tevb/tevs(jt))**(rsoft*xe/((
xb+xe)*b0))
170 IF(mstp(64).GE.2)
THEN
171 IF((1.-
z)*q2b.LT.parp(62)**2) goto 180
172 alprat=tevb/(tevb+
log(1.-
z))
173 IF(alprat.LT.5.*
rlu(0)) goto 180
174 IF(alprat.GT.5.) wtz=wtz*alprat/5.
178 IF(mstp(62).GE.3)
THEN
179 the2t=(4.*
z**2*q2b)/(vint(2)*(1.-
z)*
xb**2)
180 IF(the2t.GT.the2(jt)) goto 180
185 IF(kflb.NE.21) xfbn=xfn(kflb)
186 IF(kflb.EQ.21) xfbn=xfn(0)
187 IF(xfbn.LT.1
e-20)
THEN
188 IF(kfla.EQ.kflb)
THEN
192 ELSEIF(tevbsv-tevb.GT.0.2)
THEN
193 tevb=0.5*(tevbsv+tevb)
199 DO 210 kfl=-mstp(54),mstp(54)
200 210 xfb(kfl)=xfn(kfl)
202 CALL
pyhistfu(mint(10+jt),xa,q2ref,xfa,jt)
203 IF(kfla.NE.21) xfan=xfa(kfla)
204 IF(kfla.EQ.21) xfan=xfa(0)
205 IF(xfan.LT.1
e-20) goto 160
206 IF(kfla.NE.21) wtsfa=wtsf(kfla)
207 IF(kfla.EQ.21) wtsfa=wtsf(0)
208 IF(wtz*xfan/xfbn.LT.
rlu(0)*wtsfa) goto 160
212 220
IF(
n.EQ.
ns+2)
THEN
214 dplcm=
sqrt((dsh+dq2(1)+dq2(2))**2-4d0*dq2(1)*dq2(2))/dshr
217 IF(jr.EQ.1) ipo=ipus1
218 IF(jr.EQ.2) ipo=ipus2
227 p(i,3)=dplcm*(-1)**(jr+1)
228 p(i,4)=(dsh+dq2(3-jr)-dq2(jr))/dshr
229 p(i,5)=-
sqrt(sngl(dq2(jr)))
232 k(ipo,4)=
mod(
k(ipo,4),mstu(5))+mstu(5)*i
233 240
k(ipo,5)=
mod(
k(ipo,5),mstu(5))+mstu(5)*i
236 ELSEIF(
n.GT.
ns+2)
THEN
241 dpc(3)=0.5*(abs(
p(
is(1),3))+abs(
p(
is(2),3)))
242 dpd(1)=dsh+dq2(jr)+dq2(jt)
243 dpd(2)=dshz+dq2(jr)+dq2(3)
244 dpd(3)=
sqrt(dpd(1)**2-4d0*dq2(jr)*dq2(jt))
245 dpd(4)=
sqrt(dpd(2)**2-4d0*dq2(jr)*dq2(3))
247 IF(q2s(jr).GE.(0.5*parp(62))**2.AND.dpd(1)-dpd(3).GE.
248 & 1
d-10*dpd(1)) ikin=1
249 IF(ikin.EQ.0) dmsma=(dq2(jt)/dble(zs(jt))-dq2(3))*(dsh/
250 & (dsh+dq2(jt))-dsh/(dshz+dq2(3)))
251 IF(ikin.EQ.1) dmsma=(dpd(1)*dpd(2)-dpd(3)*dpd(4))/(2.*
252 & dq2(jr))-dq2(jt)-dq2(3)
262 IF(kflb.EQ.21.AND.kfls(jt+2).NE.21)
k(
it,2)=-kfls(jt+2)
263 IF(kflb.NE.21.AND.kfls(jt+2).EQ.21)
k(
it,2)=kflb
265 IF(sngl(dmsma).LE.
p(
it,5)**2) goto 100
266 IF(mstp(63).GE.1)
THEN
267 p(
it,4)=(dshz-dsh-
p(
it,5)**2)/dshr
269 IF(mstp(63).EQ.1)
THEN
271 ELSEIF(mstp(63).EQ.2)
THEN
272 q2tim=
min(sngl(dmsma),parp(71)*q2s(jt))
283 IF(ikin.EQ.0) dpt2=(dmsma-dms)*(dshz+dq2(3))/(dsh+dq2(jt))
284 IF(ikin.EQ.1) dpt2=(dmsma-dms)*(0.5*dpd(1)*dpd(2)+0.5*dpd(3)*
285 & dpd(4)-dq2(jr)*(dq2(jt)+dq2(3)+dms))/(4.*dsh*dpc(3)**2)
286 IF(dpt2.LT.0.) goto 100
287 dpb(1)=(0.5*dpd(2)-dpc(jr)*(dshz+dq2(jr)-dq2(jt)-dms)/
288 & dshr)/dpc(3)-dpc(3)
290 p(
it,3)=dpb(1)*(-1)**(jt+1)
291 p(
it,4)=(dshz-dsh-dms)/dshr
293 dpb(1)=
sqrt(dpb(1)**2+dpt2)
294 dpb(2)=
sqrt(dpb(1)**2+dms)
296 dpb(4)=
sqrt(dpb(3)**2+dms)
297 dbez=(dpb(4)*dpb(1)-dpb(3)*dpb(2))/(dpb(4)*dpb(2)-dpb(3)*
299 CALL ludbrb(
it+1,
n,0.,0.,0d0,0d0,dbez)
301 CALL ludbrb(
it+1,
n,the,0.,0d0,0d0,0d0)
314 p(
n+1,5)=-
sqrt(sngl(dq2(3)))
320 IF((
k(
n+1,2).GT.0.AND.
k(
n+1,2).NE.21.AND.
k(id1,2).GT.0.AND.
321 &
k(id1,2).NE.21).OR.(
k(
n+1,2).LT.0.AND.
k(id1,2).EQ.21).OR.
322 & (
k(
n+1,2).EQ.21.AND.
k(id1,2).EQ.21.AND.
rlu(0).GT.0.5).OR.
323 & (
k(
n+1,2).EQ.21.AND.
k(id1,2).LT.0)) id1=
is(jt)
325 k(
n+1,4)=
k(
n+1,4)+id1
326 k(
n+1,5)=
k(
n+1,5)+id2
327 k(id1,4)=
k(id1,4)+mstu(5)*(
n+1)
328 k(id1,5)=
k(id1,5)+mstu(5)*id2
329 k(id2,4)=
k(id2,4)+mstu(5)*id1
330 k(id2,5)=
k(id2,5)+mstu(5)*(
n+1)
334 CALL ludbrb(
ns+1,
n,0.,0.,-dble((
p(
n,1)+
p(
is(jr),1))/(
p(
n,4)+
335 &
p(
is(jr),4))),0d0,-dble((
p(
n,3)+
p(
is(jr),3))/(
p(
n,4)+
337 ir=
n+(jt-1)*(
is(1)-
n)
346 IF(mstp(62).GE.3) the2(jt)=the2t
348 IF(q2b.GE.(0.5*parp(62))**2)
THEN
354 270 xfs(jt,kfl)=xfa(kfl)
360 IF(
n.GT.mstu(4)-mstu(32)-10)
THEN
361 CALL
luerrm(11,
'(PYHISSPA:) no more memory left in LUJETS')
362 IF(mstu(21).GE.1)
n=
ns
363 IF(mstu(21).GE.1)
RETURN
365 IF(
max(q2s(1),q2s(2)).GE.(0.5*parp(62))**2.OR.
n.LE.
ns+1) goto 120
369 280 robo(j+2)=(
p(
ns+1,j)+
p(
ns+2,j))/(
p(
ns+1,4)+
p(
ns+2,4))
371 290
p(
n+2,j)=
p(
ns+1,j)
372 robot=robo(3)**2+robo(4)**2+robo(5)**2
373 IF(robot.GE.0.999999)
THEN
374 robot=1.00001*
sqrt(robot)
375 robo(3)=robo(3)/robot
376 robo(4)=robo(4)/robot
377 robo(5)=robo(5)/robot
379 CALL ludbrb(
n+2,
n+2,0.,0.,-dble(robo(3)),-dble(robo(4)),
383 CALL ludbrb(mint(83)+5,
ns,robo(1),robo(2),dble(robo(3)),
384 &dble(robo(4)),dble(robo(5)))
391 300 vint(140+jt)=xs(jt)
393 1000
FORMAT(5
x,
'structure function has a zero point here')
394 1001
FORMAT(5
x,
'xf(x,i=',i5,
')=',f10.5)