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 common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200)
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
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)
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)
90 tevb=tevb+2.*
log(alam(jt)/paru(117))
93 b0=(33.-2.*mstu(118))/6.
100 wtapq=16.*(1.-sqrt(
xb+xe))/(3.*sqrt(
xb))
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)
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
128 tevb=tevb*exp(
max(-100.,
log(
rlu(0))*b0/wtsum))
130 tevb=tevb*exp(
max(-100.,
log(
rlu(0))*b0/(5.*wtsum)))
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
151 z=
xb/(1.-
rlu(0)*(1.-sqrt(
xb+xe)))**2
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)
200 210 xfb(kfl)=xfn(kfl)
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 & 1d-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)
312 p(
n+1,3)=
p(
it,3)+
p(is(jt),3)
313 p(
n+1,4)=
p(
it,4)+
p(is(jt),4)
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)
338 CALL ludbrb(ns+1,
n,-
ulangl(
p(ir,3),
p(ir,1)),paru(2)*
rlu(0),
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)