11 IMPLICIT DOUBLE PRECISION(a-
h, o-
z)
17 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
21 common/pyctag/nct,mct(4000,2)
23 common/pyintm/kfival(2,3),nmi(2),imi(2,800,2),nvc(2,-6:6),
24 & xassoc(2,-6:6,240),xpsvc(-6:6,-1:240),pvctot(2,-1:1),
25 & xmi(2,240),pt2mi(240),imisep(0:240)
28 dimension
w(0:2,0:2),vb(3),nnxt(2),ivalq(2),icomq(2)
36 sigpt(
q)=
max(parj(21),2.1d0*
q/(7d0+
q))
41 flam(a,b,
c)=a**2+b**2+
c**2-2d0*(a*b+b*
c+
c*a)
59 IF(
mstp(91).LE.0) goto 180
68 IF (
mstp(91).GE.11.AND.
mstp(91).LE.20)
THEN
79 111
IF (
mstp(91).EQ.1.OR.
mstp(91).EQ.11)
THEN
83 ELSEIF (
mstp(91).EQ.2)
THEN
84 CALL
pyerrm(11,
'(PYMIRM:) Sorry, MSTP(91)=2 not '//
85 &
'available, using MSTP(91)=1.')
88 ELSEIF(
mstp(91).EQ.3.OR.
mstp(91).EQ.13)
THEN
94 112 p12=1d0/(1d0+27d0/40d0*
sigma**6/
eps**6)
95 IF (
pyr(0).LT.p12)
THEN
99 IF (
pyr(0).GT.wt) goto 112
104 IF (
pyr(0).GT.wt) goto 112
107 IF (
pyr(0).GT.0.5d0) msign=-1
108 IF (ixy.EQ.1) ptx=msign*
pt
109 IF (ixy.EQ.2) pty=msign*
pt
111 ELSEIF (
mstp(91).EQ.4.OR.
mstp(91).EQ.14)
THEN
112 ptx=
sigma*(sqrt(6d0)*
pyr(0)-sqrt(3d0/2d0))
113 pty=
sigma*(sqrt(6d0)*
pyr(0)-sqrt(3d0/2d0))
116 pt=sqrt(ptx**2+pty**2)
119 IF(isub.EQ.95.AND.im.EQ.1) wt=0d0
122 pt=sqrt(ptx**2+pty**2)
130 IF (mcorr.EQ.0.OR.isub.EQ.95)
THEN
131 ptcx=-ptx/(nmi(js)-1)
132 ptcy=-pty/(nmi(js)-1)
134 ptcx=-ptx/(nmi(js)-2)
135 ptcy=-pty/(nmi(js)-2)
138 IF (imc.EQ.im) goto 120
139 IF(isub.EQ.95.AND.imc.EQ.1) goto 120
140 p(imi(js,imc,1),1)=
p(imi(js,imc,1),1)+ptcx
141 p(imi(js,imc,1),2)=
p(imi(js,imc,1),2)+ptcy
143 ELSEIF (mcorr.GE.1)
THEN
148 130 imo=
k(imo,msid)/mstu(5)
149 IF (imo.EQ.0) goto 140
150 nnxt(msid-3)=nnxt(msid-3)+1
152 IF (mcorr.EQ.1.AND.
k(imo,2).EQ.21) goto 130
160 IF (nnxt(msid-3).EQ.0) goto 160
161 ptcx=-(nnxt(msid-3)*ptx)/nsum
162 ptcy=-(nnxt(msid-3)*pty)/nsum
167 fac=(1d0-rs)/(rs*(1-rs**nnxt(msid-3)))
170 imo=
k(imo,msid)/mstu(5)
171 IF (imo.EQ.0) goto 160
173 IF (
k(imo,2).NE.88)
THEN
174 p(imo,1)=
p(imo,1)+fac*ptcx
175 p(imo,2)=
p(imo,2)+fac*ptcy
176 IF (mcorr.EQ.1.AND.
k(imo,2).EQ.21) goto 150
180 l1=mod(
k(imo,4),mstu(5))
182 l3=mod(
k(imo,5),mstu(5))
183 p(l1,1)=
p(l1,1)+0.5d0*fac*ptcx
184 p(l1,2)=
p(l1,2)+0.5d0*fac*ptcy
185 p(l2,1)=
p(l2,1)+0.5d0*fac*ptcx
186 p(l2,2)=
p(l2,2)+0.5d0*fac*ptcy
187 p(l3,1)=
p(l3,1)+0.5d0*fac*ptcx
188 p(l3,2)=
p(l3,2)+0.5d0*fac*ptcy
189 p(ida,1)=
p(ida,1)-0.5d0*fac*ptcx
190 p(ida,2)=
p(ida,2)-0.5d0*fac*ptcy
201 shat=xmi(1,im)*xmi(2,im)*
vint(2)
202 pt1=sqrt(
p(imi(1,im,1),1)**2+
p(imi(1,im,1),2)**2)
203 pt2=sqrt(
p(imi(2,im,1),1)**2+
p(imi(2,im,1),2)**2)
204 pt1pt2=
p(imi(1,im,1),1)*
p(imi(2,im,1),1)
205 & +
p(imi(1,im,1),2)*
p(imi(2,im,1),2)
206 IF (shat.LT.2d0*(pt1*
pt2-pt1pt2).AND.
ntry.LE.100)
THEN
210 &
'(PYMIRM:) No consistent (x,kT) sets found')
222 pt2 = (
p(imi(1,im,1),1)+
p(imi(2,im,1),1))**2
223 & +(
p(imi(1,im,1),2)+
p(imi(2,im,1),2))**2
224 st=xmi(1,im)*xmi(2,im)*
vint(2)+
pt2
225 w(0,1)=
w(0,1)-sqrt(xmi(1,im)/xmi(2,im)*st)
226 w(0,2)=
w(0,2)-sqrt(xmi(2,im)/xmi(1,im)*st)
231 IF (
w(0,0).LT.0d0.AND.
ntry.LE.100)
THEN
235 &
'(PYMIRM:) Negative beam remnant mass squared unavoidable')
249 DO 270 im=
mint(31)+1,nmi(js)
256 IF (kfa.EQ.21.AND.
k(
i,1).EQ.14) goto 270
257 IF (kfa.EQ.88) goto 270
268 ELSEIF(kfa.LE.6.AND.icomp.EQ.0)
THEN
271 ELSEIF(kfa.LE.6)
THEN
274 ELSEIF(kfa.GT.1000.AND.mod(kfa/10,10).EQ.0.AND.
276 ivalq(1)=isign(kfa/1000,kf)
277 ivalq(2)=isign(mod(kfa/100,10),kf)
279 ELSEIF(kfa.GT.1000.AND.mod(kfa/10,10).EQ.0.AND.
280 & icomp.LT.mstu(5))
THEN
281 IF(kfa/1000.EQ.iabs(
k(icomp,2)))
THEN
282 ivalq(1)=isign(mod(kfa/100,10),kf)
284 ivalq(1)=isign(kfa/1000,kf)
289 DO 220 im1=1,
mint(31)
290 IF(imi(js,im1,2).EQ.
i.AND.imi(js,im1,1).NE.icomp)
THEN
291 icomq(2)=imi(js,im1,1)
296 ELSEIF(kfa.GT.1000.AND.mod(kfa/10,10).EQ.0)
THEN
297 icomq(1)=mod(icomp,mstu(5))
298 icomq(2)=icomp/mstu(5)
303 IF(mod(kfa/1000,10).EQ.0)
THEN
306 kfl1=mod(kfa,10000)-10*kfl3-1
307 IF(mod(kfa/1000,10).EQ.mod(kfa/100,10).AND.
308 & mod(kfa,10).EQ.2) kfl1=kfl1+2
310 pr=
p(
i,5)**2+
p(
i,1)**2+
p(
i,2)**2
318 IF(ivalq(iq).NE.0)
THEN
320 IF(kfival(js,1).EQ.ivalq(iq)) nval=nval+1
321 IF(kfival(js,2).EQ.ivalq(iq)) nval=nval+1
322 IF(kfival(js,3).EQ.ivalq(iq)) nval=nval+1
324 IF(kfival(js,3).EQ.0)
THEN
327 ELSEIF(nval.EQ.3)
THEN
328 mdu=int(
pyr(0)+5d0/3d0)
330 ELSEIF(nval.EQ.2)
THEN
333 ELSEIF(kfival(js,1).EQ.kfival(js,2).OR.kfival(js,2).EQ.
334 & kfival(js,3).OR.kfival(js,1).EQ.kfival(js,3))
THEN
338 mdu=int(
pyr(0)+5d0/3d0)
341 IF(mdu.EQ.1) xpow=3.5d0
342 IF(mdu.EQ.2) xpow=2d0
344 IF((1d0-xx)**xpow.LT.
pyr(0)) goto 230
349 IF(icomq(iq).NE.0)
THEN
351 DO 240 im1=1,
mint(31)
352 IF(imi(js,im1,1).EQ.icomq(iq)) xcomp=xmi(js,im1)
355 250 xx=xcomp*(1d0/(1d0-
pyr(0)*(1d0-xcomp))-1d0)
356 corr=((1d0-xcomp-xx)/(1d0-xcomp))**npow*
357 & (xcomp**2+xx**2)/(xcomp+xx)**2
364 IF (kfa.GT.100)
x=
parp(79)*
x
369 w(js,3-js)=
w(js,3-js)+(
p(
i,5)**2+
p(
i,1)**2+
p(
i,2)**2)/
x
371 w(js,js)=
w(js,js)*
w(0,js)
372 w(js,3-js)=
w(js,3-js)/
w(0,js)
373 w(js,0)=
w(js,1)*
w(js,2)
378 IF (sqrt(
w(1,0))+sqrt(
w(2,0)).GT.sqrt(
w(0,0)).AND.ntryx.LE.100)
381 ELSEIF (ntryx.GT.100.AND.
ntry.LE.100)
THEN
383 ELSEIF (ntryx.GT.100)
THEN
384 CALL
pyerrm(11,
'(PYMIRM:) No consistent (x,kT) sets found')
391 comtrm=
w(0,0)+sqrt(flam(
w(0,0),
w(1,0),
w(2,0)))
392 r1=(comtrm+
w(1,0)-
w(2,0))/(2d0*
w(1,1)*
w(0,2))
393 r2=(comtrm+
w(2,0)-
w(1,0))/(2d0*
w(2,2)*
w(0,1))
395 IF (
r1.LT.0.OR.
r2.LT.0)
THEN
396 CALL
pyerrm(19,
'(PYMIRM:) negative rescaling factors !')
408 DO 290 im=
mint(31)+1,
max(nmi(1),nmi(2))
409 xmi(1,im)=xmi(1,im)*
r1
410 xmi(2,im)=xmi(2,im)*
r2
418 st=xmi(1,im)*xmi(2,im)*
vint(2)+(
p(
i1,1)+
p(
i2,1))**2+
420 pt12=
p(
i1,1)**2+
p(
i1,2)**2
421 pt22=
p(
i2,1)**2+
p(
i2,2)**2
423 p(
i1,3)=sqrt(flam(st,pt12,pt22)/(4d0*st))
426 p(
i1,4)=sqrt(pt12+
p(
i1,3)**2)
427 p(
i2,4)=sqrt(pt22+
p(
i2,3)**2)
430 vb(1)=(
p(
i1,1)+
p(
i2,1))/sqrt(st)
431 vb(2)=(
p(
i1,2)+
p(
i2,2))/sqrt(st)
445 vb(3)=(xmi(1,im)-xmi(2,im))/(xmi(1,im)+xmi(2,im))
447 IF (im.EQ.1) imin=
mint(83)+5
449 CALL
pyrobo(imin,imax,the,
phi,vb(1),vb(2),0d0)
450 CALL
pyrobo(imin,imax,0d0,0d0,0d0,0d0,vb(3))
457 DO 310 im=
mint(31)+1,nmi(js)
460 IF (
k(
i,2).EQ.21.AND.
k(
i,1).EQ.14) goto 310
461 IF (kfa.EQ.88) goto 310
462 rmt2=
p(
i,5)**2+
p(
i,1)**2+
p(
i,2)**2
463 p(
i,4)=0.5d0*(xmi(js,im)*
w(0,js)+rmt2/(xmi(js,im)*
w(0,js)))
464 p(
i,3)=0.5d0*(xmi(js,im)*
w(0,js)-rmt2/(xmi(js,im)*
w(0,js)))
465 IF (js.EQ.2)
p(
i,3)=-
p(
i,3)
488 IF (
mstp(95).NE.1.OR.
mint(31).LE.1) goto 380
500 IF (iiter.LE.
parp(78)*ntot)
THEN
505 kct=mod(int(jct+
pyr(0)*nct),nct)+1
512 IF (
k(
i,1).EQ.3)
THEN
513 IF (mct(
i,1).EQ.jct) ij1=
i
514 IF (mct(
i,2).EQ.jct) ij2=
i
515 IF (mct(
i,1).EQ.kct) ik1=
i
516 IF (mct(
i,2).EQ.kct) ik2=
i
520 IF (ij1.EQ.0.OR.ij2.EQ.0.OR.ik1.EQ.0.OR.ik2.EQ.0) goto 360
522 rlo=2d0*four(ij1,ij2)*2d0*four(ik1,ik2)
523 rln=2d0*four(ij1,ik2)*2d0*four(ik1,ij2)
524 IF (rln.LT.rlo.AND.mct(ij2,1).NE.kct.AND.mct(ik2,1).NE.jct)
THEN
530 IF (
mint(34).LE.1000)
THEN
533 CALL
pyerrm(4,
'(PYMIRM:) caught in infinite loop')
537 IF (nrecp.LT.
mint(34)) goto 350