11 common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
13 common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
15 common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200)
17 common/pyhipars/mstp(200),parp(200),msti(200),pari(200)
19 common/pyhiint1/mint(400),vint(400)
21 common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
23 common/pyhiint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
25 common/pyhiint5/ngen(0:200,3),xsec(0:200,3)
27 dimension nmul(20),sigm(20),kstr(500,2)
28 SAVE xt2,xt2fac,xc2,xts,irbin,rbin,nmul,sigm
32 IF(mstp(122).GE.1)
WRITE(mstu(11),1000) mstp(82)
45 DO 110 itry=1,mstp(83)
46 rsca=0.05*((21-ixt2)-
rlu(0))
47 xt2=vint(149)*(1.+vint(149))/(vint(149)+rsca)-vint(149)
48 xt2=
max(0.01*vint(149),xt2)
52 IF(
rlu(0).LE.coef(isub,1))
THEN
53 taup=(2.*(1.+
sqrt(1.-xt2))/xt2-1.)**
rlu(0)
54 tau=xt2*(1.+taup)**2/(4.*taup)
62 IF(ryst.GT.coef(isub,7)) myst=2
63 IF(ryst.GT.coef(isub,7)+coef(isub,8)) myst=3
68 vint(71)=0.5*vint(1)*
sqrt(xt2)
70 110 sigm(ixt2)=sigm(ixt2)+sigs
71 120 sigsum=sigsum+sigm(ixt2)
72 sigsum=sigsum/(20.*mstp(83))
75 IF(sigsum.LT.1.1*vint(106))
THEN
76 IF(mstp(122).GE.1)
WRITE(mstu(11),1100) parp(82),sigsum
78 vint(149)=4.*parp(82)**2/vint(2)
81 IF(mstp(122).GE.1)
WRITE(mstu(11),1200) parp(82), sigsum
95 xk=xi+(yke-yi)*(xf-xi)/(yf-yi)
99 IF(mstp(82).EQ.2)
THEN
100 sp=0.5*paru(1)*(1.-
exp(-xk))
103 IF(mstp(82).EQ.3) deltab=0.02
104 IF(mstp(82).EQ.4) deltab=
min(0.01,0.05*parp(84))
109 IF(mstp(82).EQ.3)
THEN
110 ov=
exp(-
b**2)/paru(2)
113 ov=((1.-parp(83))**2*
exp(-
min(100.,
b**2))+2.*parp(83)*
114 & (1.-parp(83))*2./(1.+cq2)*
exp(-
min(100.,
b**2*2./(1.+cq2)))+
115 & parp(83)**2/cq2*
exp(-
min(100.,
b**2/cq2)))/paru(2)
117 pacc=1.-
exp(-
min(100.,paru(1)*xk*ov))
118 sp=sp+paru(2)*
b*deltab*pacc
119 sop=sop+paru(2)*
b*deltab*ov*pacc
120 IF(
b.LT.1..OR.
b*pacc.GT.1
e-6) goto 140
134 IF(abs(yk-yke).GE.1
e-5*yke) goto 130
142 ELSEIF(mmul.EQ.2)
THEN
143 IF(mstp(82).LE.0)
THEN
144 ELSEIF(mstp(82).EQ.1)
THEN
146 xt2fac=xsec(96,1)/vint(106)*vint(149)/(1.-vint(149))
147 ELSEIF(mstp(82).EQ.2)
THEN
149 xt2fac=vint(146)*xsec(96,1)/vint(106)*vint(149)*(1.+vint(149))
151 xc2=4.*ckin(3)**2/vint(2)
152 IF(ckin(3).LE.ckin(5).OR.mint(82).GE.2) xc2=0.
155 ELSEIF(mmul.EQ.3)
THEN
160 IF(mstp(82).LE.0)
THEN
162 ELSEIF(mstp(82).EQ.1)
THEN
163 xt2=xt2fac*xt2/(xt2fac-xt2*
log(
rlu(0)))
164 ELSEIF(mstp(82).EQ.2)
THEN
165 IF(xt2.LT.1..AND.
exp(-xt2fac*xt2/(vint(149)*(xt2+
166 & vint(149)))).GT.
rlu(0)) xt2=1.
168 xt2=(1.+vint(149))*xt2fac/(xt2fac-(1.+vint(149))*
log(1.-
169 &
rlu(0)*(1.-
exp(-xt2fac/(vint(149)*(1.+vint(149)))))))-
172 xt2=-xt2fac/
log(
exp(-xt2fac/(xt2+vint(149)))+
rlu(0)*
173 & (
exp(-xt2fac/vint(149))-
exp(-xt2fac/(xt2+vint(149)))))-
176 xt2=
max(0.01*vint(149),xt2)
178 xt2=(xc2+vint(149))*(1.+vint(149))/(1.+vint(149)-
179 &
rlu(0)*(1.-xc2))-vint(149)
180 xt2=
max(0.01*vint(149),xt2)
185 IF(mstp(82).LE.1.AND.xt2.LT.vint(149))
THEN
186 IF(mint(82).EQ.1) ngen(0,1)=ngen(0,1)-1
187 IF(mint(82).EQ.1) ngen(isub,1)=ngen(isub,1)-1
190 vint(21)=0.01*vint(149)
193 vint(25)=0.01*vint(149)
198 IF(
rlu(0).LE.coef(isub,1))
THEN
199 taup=(2.*(1.+
sqrt(1.-xt2))/xt2-1.)**
rlu(0)
200 tau=xt2*(1.+taup)**2/(4.*taup)
208 IF(ryst.GT.coef(isub,7)) myst=2
209 IF(ryst.GT.coef(isub,7)+coef(isub,8)) myst=3
213 vint(71)=0.5*vint(1)*
sqrt(vint(25))
216 ELSEIF(mmul.EQ.4)
THEN
219 IF(iset(isub).EQ.1) xts=vint(21)
220 IF(iset(isub).EQ.2) xts=(4.*vint(48)+2.*vint(63)+2.*vint(64))/
222 IF(iset(isub).EQ.3.OR.iset(isub).EQ.4) xts=vint(26)
223 rbin=
max(0.000001,
min(0.999999,xts*(1.+vint(149))/
225 irbin=
int(1.+20.*rbin)
226 IF(isub.EQ.96) nmul(irbin)=nmul(irbin)+1
227 IF(isub.EQ.96) sigm(irbin)=sigm(irbin)+vint(153)
230 ELSEIF(mmul.EQ.5)
THEN
231 IF(mstp(82).EQ.3)
THEN
232 vint(148)=
rlu(0)/(paru(2)*vint(147))
236 IF(rtype.LT.(1.-parp(83))**2)
THEN
238 ELSEIF(rtype.LT.1.-parp(83)**2)
THEN
239 b2=-0.5*(1.+cq2)*
log(
rlu(0))
243 vint(148)=((1.-parp(83))**2*
exp(-
min(100.,b2))+2.*parp(83)*
244 & (1.-parp(83))*2./(1.+cq2)*
exp(-
min(100.,b2*2./(1.+cq2)))+
245 & parp(83)**2/cq2*
exp(-
min(100.,b2/cq2)))/(paru(2)*vint(147))
250 rncor=(irbin-20.*rbin)*nmul(irbin)
251 sigcor=(irbin-20.*rbin)*sigm(irbin)
252 DO 150 ibin=irbin+1,20
253 rncor=rncor+nmul(ibin)
254 150 sigcor=sigcor+sigm(ibin)
255 sigabv=(sigcor/rncor)*vint(149)*(1.-xts)/(xts+vint(149))
256 vint(150)=
exp(-
min(100.,vint(146)*vint(148)*sigabv/vint(106)))
259 ELSEIF(mmul.EQ.6)
THEN
264 IF(iset(isub).EQ.1)
nmax=mint(84)+2
266 DO 170 i=mint(84)+1,
nmax
267 kcs=kchg(
lucomp(
k(i,2)),2)*isign(1,
k(i,2))
268 IF(kcs.EQ.0) goto 170
270 IF(kcs.EQ.1.AND.(j.EQ.2.OR.j.EQ.4)) goto 160
271 IF(kcs.EQ.-1.AND.(j.EQ.1.OR.j.EQ.3)) goto 160
273 ist=
mod(
k(i,j+3)/mstu(5),mstu(5))
275 ist=
mod(
k(i,j+1),mstu(5))
277 IF(ist.LT.mint(84).OR.ist.GT.i) goto 160
278 IF(kchg(
lucomp(
k(ist,2)),2).EQ.0) goto 160
280 IF(j.EQ.1.OR.j.EQ.4)
THEN
292 IF(iset(isub).EQ.1) xt2=vint(21)
293 IF(iset(isub).EQ.2) xt2=(4.*vint(48)+2.*vint(63)+2.*vint(64))/
295 IF(iset(isub).EQ.3.OR.iset(isub).EQ.4) xt2=vint(26)
298 IF(mstp(82).LE.1)
THEN
299 xt2fac=xsec(isub,1)*vint(149)/((1.-vint(149))*vint(106))
301 xt2fac=vint(146)*vint(148)*xsec(isub,1)/vint(106)*
302 & vint(149)*(1.+vint(149))
308 vint(143)=1.-vint(141)
309 vint(144)=1.-vint(142)
312 180
IF(mstp(82).LE.1)
THEN
313 xt2=xt2fac*xt2/(xt2fac-xt2*
log(
rlu(0)))
314 IF(xt2.LT.vint(149)) goto 220
316 IF(xt2.LE.0.01*vint(149)) goto 220
317 xt2=xt2fac*(xt2+vint(149))/(xt2fac-(xt2+vint(149))*
319 IF(xt2.LE.0.) goto 220
320 xt2=
max(0.01*vint(149),xt2)
325 IF(
rlu(0).LE.coef(isub,1))
THEN
326 taup=(2.*(1.+
sqrt(1.-xt2))/xt2-1.)**
rlu(0)
327 tau=xt2*(1.+taup)**2/(4.*taup)
335 IF(ryst.GT.coef(isub,7)) myst=2
336 IF(ryst.GT.coef(isub,7)+coef(isub,8)) myst=3
343 IF(vint(143)-x1m.LT.0.01.OR.vint(144)-x2m.LT.0.01) goto 180
344 vint(71)=0.5*vint(1)*
sqrt(xt2)
346 IF(sigs.LT.xsec(isub,1)*
rlu(0)) goto 180
362 IF(rflav.GE.
max(parp(85),parp(86)))
k(
n+1,2)=
363 & 1+
int((2.+parj(2))*
rlu(0))
366 p(
n+1,3)=0.25*vint(1)*(vint(41)*(1.+cth)-vint(42)*(1.-cth))
367 p(
n+1,4)=0.25*vint(1)*(vint(41)*(1.+cth)+vint(42)*(1.-cth))
373 IF(
k(
n+1,2).NE.21)
k(
n+2,2)=-
k(
n+1,2)
376 p(
n+2,3)=0.25*vint(1)*(vint(41)*(1.-cth)-vint(42)*(1.+cth))
377 p(
n+2,4)=0.25*vint(1)*(vint(41)*(1.-cth)+vint(42)*(1.+cth))
380 IF(rflav.LT.parp(85).AND.nstr.GE.1)
THEN
387 dist=(
p(i,4)*
p(i1,4)-
p(i,1)*
p(i1,1)-
p(i,2)*
p(i1,2)-
388 &
p(i,3)*
p(i1,3))*(
p(i,4)*
p(i2,4)-
p(i,1)*
p(i2,1)-
389 &
p(i,2)*
p(i2,2)-
p(i,3)*
p(i2,3))/
max(1.,
p(i1,4)*
p(i2,4)-
390 &
p(i1,1)*
p(i2,1)-
p(i1,2)*
p(i2,2)-
p(i1,3)*
p(i2,3))
391 IF(istr.EQ.1.OR.dist.LT.dmin)
THEN
400 IF(
k(ist1,4)/mstu(5).EQ.ist2)
k(ist1,4)=mstu(5)*i+
401 &
mod(
k(ist1,4),mstu(5))
402 IF(
mod(
k(ist1,5),mstu(5)).EQ.ist2)
k(ist1,5)=
403 & mstu(5)*(
k(ist1,5)/mstu(5))+i
406 IF(
k(ist2,5)/mstu(5).EQ.ist1)
k(ist2,5)=mstu(5)*i+
407 &
mod(
k(ist2,5),mstu(5))
408 IF(
mod(
k(ist2,4),mstu(5)).EQ.ist1)
k(ist2,4)=
409 & mstu(5)*(
k(ist2,4)/mstu(5))+i
416 ELSEIF(
k(
n+1,2).EQ.21)
THEN
417 k(
n+1,4)=mstu(5)*(
n+2)
418 k(
n+1,5)=mstu(5)*(
n+2)
419 k(
n+2,4)=mstu(5)*(
n+1)
420 k(
n+2,5)=mstu(5)*(
n+1)
429 k(
n+1,4)=mstu(5)*(
n+2)
430 k(
n+2,5)=mstu(5)*(
n+1)
438 IF(
n.GT.mstu(4)-mstu(32)-10)
THEN
439 CALL
luerrm(11,
'(PYHIMULT:) no more memory left in LUJETS')
440 IF(mstu(21).GE.1)
RETURN
443 vint(151)=vint(151)+vint(41)
444 vint(152)=vint(152)+vint(42)
445 vint(143)=vint(143)-vint(41)
446 vint(144)=vint(144)-vint(42)
447 IF(mint(31).LT.240) goto 180
452 1000
FORMAT(/1
x,
'****** PYHIMULT: initialization of multiple inter',
453 &
'actions for MSTP(82) =',i2,
' ******')
454 1100
FORMAT(8
x,
'pT0 =',f5.2,
' GeV gives sigma(parton-parton) =',1
p,
455 &e9.2,
' mb: rejected')
456 1200
FORMAT(8
x,
'pT0 =',f5.2,
' GeV gives sigma(parton-parton) =',1
p,
457 &e9.2,
' mb: accepted')