pythiaeRHIC/pythia/pymign.f

0001
0002 C*********************************************************************
0003
0004 C...PYMIGN
0005 C...Initializes treatment of new multiple interactions scenario,
0006 C...selects kinematics of hardest interaction if low-pT physics
0007 C...included in run, and generates all non-hardest interactions.
0008
0009       SUBROUTINE PYMIGN(MMUL)
0010
0011 C...Double precision and integer declarations.
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013       IMPLICIT INTEGER(I-N)
0014       INTEGER PYK,PYCHGE,PYCOMP
0015       EXTERNAL PYALPS
0016       DOUBLE PRECISION PYALPS
0017 C...Commonblocks.
0018       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0019       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0020       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0021       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0022       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0023       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0024       COMMON/PYINT1/MINT(400),VINT(400)
0025       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0026       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
0027       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
0028       COMMON/PYINT7/SIGT(0:6,0:6,0:5)
0029       COMMON/PYINTM/KFIVAL(2,3),NMI(2),IMI(2,800,2),NVC(2,-6:6),
0030      &     XASSOC(2,-6:6,240),XPSVC(-6:6,-1:240),PVCTOT(2,-1:1),
0031      &     XMI(2,240),PT2MI(240),IMISEP(0:240)
0032       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYSUBS/,/PYPARS/,
0033      &/PYINT1/,/PYINT2/,/PYINT3/,/PYINT5/,/PYINT7/,/PYINTM/
0034 C...Local arrays and saved variables.
0035       DIMENSION NMUL(20),SIGM(20),KSTR(500,2),VINTSV(80),
0036      &WDTP(0:400),WDTE(0:400,0:5),XPQ(-25:25),KSAV(4,5),PSAV(4,5)
0037       SAVE XT2,XT2FAC,XC2,XTS,IRBIN,RBIN,NMUL,SIGM,P83A,P83B,P83C,
0038      &CQ2I,CQ2R,PIK,BDIV,B,PLOWB,PHIGHB,PALLB,S4A,S4B,S4C,POWIP,
0039      &RPWIP,B2RPDV,B2RPMX,BAVG,VNT145,VNT146,VNT147
0040
0041 C...Initialization of multiple interaction treatment.
0042       IF(MMUL.EQ.1) THEN
0043         IF(MSTP(122).GE.1) WRITE(MSTU(11),5000) MSTP(82)
0044         ISUB=96
0045         MINT(1)=96
0046         VINT(63)=0D0
0047         VINT(64)=0D0
0048         VINT(143)=1D0
0049         VINT(144)=1D0
0050
0051 C...Loop over phase space points: xT2 choice in 20 bins.
0052   100   SIGSUM=0D0
0053         DO 120 IXT2=1,20
0054           NMUL(IXT2)=MSTP(83)
0055           SIGM(IXT2)=0D0
0056           DO 110 ITRY=1,MSTP(83)
0057             RSCA=0.05D0*((21-IXT2)-PYR(0))
0058             XT2=VINT(149)*(1D0+VINT(149))/(VINT(149)+RSCA)-VINT(149)
0059             XT2=MAX(0.01D0*VINT(149),XT2)
0060             VINT(25)=XT2
0061
0062 C...Choose tau and y*. Calculate cos(theta-hat).
0063             IF(PYR(0).LE.COEF(ISUB,1)) THEN
0064               TAUT=(2D0*(1D0+SQRT(1D0-XT2))/XT2-1D0)**PYR(0)
0065               TAU=XT2*(1D0+TAUT)**2/(4D0*TAUT)
0066             ELSE
0067               TAU=XT2*(1D0+TAN(PYR(0)*ATAN(SQRT(1D0/XT2-1D0)))**2)
0068             ENDIF
0069             VINT(21)=TAU
0070             CALL PYKLIM(2)
0071             RYST=PYR(0)
0072             MYST=1
0073             IF(RYST.GT.COEF(ISUB,8)) MYST=2
0074             IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
0075             CALL PYKMAP(2,MYST,PYR(0))
0076             VINT(23)=SQRT(MAX(0D0,1D0-XT2/TAU))*(-1)**INT(1.5D0+PYR(0))
0077
0078 C...Calculate differential cross-section.
0079             VINT(71)=0.5D0*VINT(1)*SQRT(XT2)
0080             CALL PYSIGH(NCHN,SIGS)
0081             SIGM(IXT2)=SIGM(IXT2)+SIGS
0082   110     CONTINUE
0083           SIGSUM=SIGSUM+SIGM(IXT2)
0084   120   CONTINUE
0085         SIGSUM=SIGSUM/(20D0*MSTP(83))
0086
0087 C...Reject result if sigma(parton-parton) is smaller than hadronic one.
0088         IF(SIGSUM.LT.1.1D0*SIGT(0,0,5)) THEN
0089           IF(MSTP(122).GE.1) WRITE(MSTU(11),5100)
0090      &    PARP(82)*(VINT(1)/PARP(89))**PARP(90),SIGSUM
0091           PARP(82)=0.9D0*PARP(82)
0092           VINT(149)=4D0*(PARP(82)*(VINT(1)/PARP(89))**PARP(90))**2/
0093      &    VINT(2)
0094           GOTO 100
0095         ENDIF
0096         IF(MSTP(122).GE.1) WRITE(MSTU(11),5200)
0097      &  PARP(82)*(VINT(1)/PARP(89))**PARP(90), SIGSUM
0098
0099 C...Start iteration to find k factor.
0100         YKE=SIGSUM/MAX(1D-10,SIGT(0,0,5))
0101         P83A=(1D0-PARP(83))**2
0102         P83B=2D0*PARP(83)*(1D0-PARP(83))
0103         P83C=PARP(83)**2
0104         CQ2I=1D0/PARP(84)**2
0105         CQ2R=2D0/(1D0+PARP(84)**2)
0106         SO=0.5D0
0107         XI=0D0
0108         YI=0D0
0109         XF=0D0
0110         YF=0D0
0111         XK=0.5D0
0112         IIT=0
0113   130   IF(IIT.EQ.0) THEN
0114           XK=2D0*XK
0115         ELSEIF(IIT.EQ.1) THEN
0116           XK=0.5D0*XK
0117         ELSE
0118           XK=XI+(YKE-YI)*(XF-XI)/(YF-YI)
0119         ENDIF
0120
0121 C...Evaluate overlap integrals. Find where to divide the b range.
0122         IF(MSTP(82).EQ.2) THEN
0123           SP=0.5D0*PARU(1)*(1D0-EXP(-XK))
0124           SOP=SP/PARU(1)
0125         ELSE
0126           IF(MSTP(82).EQ.3) THEN
0127             DELTAB=0.02D0
0128           ELSEIF(MSTP(82).EQ.4) THEN
0129             DELTAB=MIN(0.01D0,0.05D0*PARP(84))
0130           ELSE
0131             POWIP=MAX(0.4D0,PARP(83))
0132             RPWIP=2D0/POWIP-1D0
0133             DELTAB=MAX(0.02D0,0.02D0*(2D0/POWIP)**(1D0/POWIP))
0134             SO=0D0
0135           ENDIF
0136           SP=0D0
0137           SOP=0D0
0138           BSP=0D0
0139           SOHIGH=0D0
0140           IBDIV=0
0141           B=-0.5D0*DELTAB
0142   140     B=B+DELTAB
0143           IF(MSTP(82).EQ.3) THEN
0144             OV=EXP(-B**2)/PARU(2)
0145           ELSEIF(MSTP(82).EQ.4) THEN
0146             OV=(P83A*EXP(-MIN(50D0,B**2))+
0147      &      P83B*CQ2R*EXP(-MIN(50D0,B**2*CQ2R))+
0148      &      P83C*CQ2I*EXP(-MIN(50D0,B**2*CQ2I)))/PARU(2)
0149           ELSE
0150             OV=EXP(-B**POWIP)/PARU(2)
0151             SO=SO+PARU(2)*B*DELTAB*OV
0152           ENDIF
0153           IF(IBDIV.EQ.1) SOHIGH=SOHIGH+PARU(2)*B*DELTAB*OV
0154           PACC=1D0-EXP(-MIN(50D0,PARU(1)*XK*OV))
0155           SP=SP+PARU(2)*B*DELTAB*PACC
0156           SOP=SOP+PARU(2)*B*DELTAB*OV*PACC
0157           BSP=BSP+B*PARU(2)*B*DELTAB*PACC
0158           IF(IBDIV.EQ.0.AND.PARU(1)*XK*OV.LT.1D0) THEN
0159             IBDIV=1
0160             BDIV=B+0.5D0*DELTAB
0161           ENDIF
0162           IF(B.LT.1D0.OR.B*PACC.GT.1D-6) GOTO 140
0163         ENDIF
0164         YK=PARU(1)*XK*SO/SP
0165
0166 C...Continue iteration until convergence.
0167         IF(YK.LT.YKE) THEN
0168           XI=XK
0169           YI=YK
0170           IF(IIT.EQ.1) IIT=2
0171         ELSE
0172           XF=XK
0173           YF=YK
0174           IF(IIT.EQ.0) IIT=1
0175         ENDIF
0176         IF(ABS(YK-YKE).GE.1D-5*YKE) GOTO 130
0177
0178 C...Store some results for subsequent use.
0179         BAVG=BSP/SP
0180         VINT(145)=SIGSUM
0181         VINT(146)=SOP/SO
0182         VINT(147)=SOP/SP
0183         VNT145=VINT(145)
0184         VNT146=VINT(146)
0185         VNT147=VINT(147)
0186 C...PIK = PARU(1)*XK = (VINT(146)/VINT(147))*sigma_jet/sigma_nondiffr.
0187         PIK=(VNT146/VNT147)*YKE
0188
0189 C...Find relative weight for low and high impact parameter..
0190       PLOWB=PARU(1)*BDIV**2
0191       IF(MSTP(82).EQ.3) THEN
0192         PHIGHB=PIK*0.5*EXP(-BDIV**2)
0193       ELSEIF(MSTP(82).EQ.4) THEN
0194         S4A=P83A*EXP(-BDIV**2)
0195         S4B=P83B*EXP(-BDIV**2*CQ2R)
0196         S4C=P83C*EXP(-BDIV**2*CQ2I)
0197         PHIGHB=PIK*0.5*(S4A+S4B+S4C)
0198       ELSEIF(PARP(83).GE.1.999D0) THEN
0199         PHIGHB=PIK*SOHIGH
0200         B2RPDV=BDIV**POWIP
0201       ELSE
0202         PHIGHB=PIK*SOHIGH
0203         B2RPDV=BDIV**POWIP
0204         B2RPMX=MAX(2D0*RPWIP,B2RPDV)
0205       ENDIF
0206       PALLB=PLOWB+PHIGHB
0207
0208 C...Initialize iteration in xT2 for hardest interaction.
0209       ELSEIF(MMUL.EQ.2) THEN
0210         VINT(145)=VNT145
0211         VINT(146)=VNT146
0212         VINT(147)=VNT147
0213         IF(MSTP(82).LE.0) THEN
0214         ELSEIF(MSTP(82).EQ.1) THEN
0215           XT2=1D0
0216           SIGRAT=XSEC(96,1)/MAX(1D-10,VINT(315)*VINT(316)*SIGT(0,0,5))
0217           IF(MINT(141).NE.0.OR.MINT(142).NE.0) SIGRAT=SIGRAT*
0218      &    VINT(317)/(VINT(318)*VINT(320))
0219           XT2FAC=SIGRAT*VINT(149)/(1D0-VINT(149))
0220         ELSEIF(MSTP(82).EQ.2) THEN
0221           XT2=1D0
0222           XT2FAC=VNT146*XSEC(96,1)/MAX(1D-10,SIGT(0,0,5))*
0223      &    VINT(149)*(1D0+VINT(149))
0224         ELSE
0225           XC2=4D0*CKIN(3)**2/VINT(2)
0226           IF(CKIN(3).LE.CKIN(5).OR.MINT(82).GE.2) XC2=0D0
0227         ENDIF
0228
0229 C...Select impact parameter for hardest interaction.
0230         IF(MSTP(82).LE.2) RETURN
0231   142   IF(PYR(0)*PALLB.LT.PLOWB) THEN
0232 C...Treatment in low b region.
0233           MINT(39)=1
0234           B=BDIV*SQRT(PYR(0))
0235           IF(MSTP(82).EQ.3) THEN
0236             OV=EXP(-B**2)/PARU(2)
0237           ELSEIF(MSTP(82).EQ.4) THEN
0238             OV=(P83A*EXP(-MIN(50D0,B**2))+
0239      &      P83B*CQ2R*EXP(-MIN(50D0,B**2*CQ2R))+
0240      &      P83C*CQ2I*EXP(-MIN(50D0,B**2*CQ2I)))/PARU(2)
0241           ELSE
0242             OV=EXP(-B**POWIP)/PARU(2)
0243           ENDIF
0244           VINT(148)=OV/VNT147
0245           PACC=1D0-EXP(-MIN(50D0,PIK*OV))
0246           XT2=1D0
0247           XT2FAC=VNT146*VINT(148)*XSEC(96,1)/MAX(1D-10,SIGT(0,0,5))*
0248      &    VINT(149)*(1D0+VINT(149))
0249         ELSE
0250 C...Treatment in high b region.
0251           MINT(39)=2
0252           IF(MSTP(82).EQ.3) THEN
0253             B=SQRT(BDIV**2-LOG(PYR(0)))
0254             OV=EXP(-B**2)/PARU(2)
0255           ELSEIF(MSTP(82).EQ.4) THEN
0256             S4RNDM=PYR(0)*(S4A+S4B+S4C)
0257             IF(S4RNDM.LT.S4A) THEN
0258               B=SQRT(BDIV**2-LOG(PYR(0)))
0259             ELSEIF(S4RNDM.LT.S4A+S4B) THEN
0260               B=SQRT(BDIV**2-LOG(PYR(0))/CQ2R)
0261             ELSE
0262               B=SQRT(BDIV**2-LOG(PYR(0))/CQ2I)
0263             ENDIF
0264             OV=(P83A*EXP(-MIN(50D0,B**2))+
0265      &      P83B*CQ2R*EXP(-MIN(50D0,B**2*CQ2R))+
0266      &      P83C*CQ2I*EXP(-MIN(50D0,B**2*CQ2I)))/PARU(2)
0267           ELSEIF(PARP(83).GE.1.999D0) THEN
0268   144       B2RPW=B2RPDV-LOG(PYR(0))
0269             ACCIP=(B2RPW/B2RPDV)**RPWIP
0270             IF(ACCIP.LT.PYR(0)) GOTO 144
0271             OV=EXP(-B2RPW)/PARU(2)
0272             B=B2RPW**(1D0/POWIP)
0273           ELSE
0274   146       B2RPW=B2RPDV-2D0*LOG(PYR(0))
0275             ACCIP=(B2RPW/B2RPMX)**RPWIP*EXP(-0.5D0*(B2RPW-B2RPMX))
0276             IF(ACCIP.LT.PYR(0)) GOTO 146
0277             OV=EXP(-B2RPW)/PARU(2)
0278             B=B2RPW**(1D0/POWIP)
0279           ENDIF
0280           VINT(148)=OV/VNT147
0281           PACC=(1D0-EXP(-MIN(50D0,PIK*OV)))/(PIK*OV)
0282         ENDIF
0283         IF(PACC.LT.PYR(0)) GOTO 142
0284         VINT(139)=B/BAVG
0285
0286       ELSEIF(MMUL.EQ.3) THEN
0287 C...Low-pT or multiple interactions (first semihard interaction):
0288 C...choose xT2 according to dpT2/pT2**2*exp(-(sigma above pT2)/norm)
0289 C...or (MSTP(82)>=2) dpT2/(pT2+pT0**2)**2*exp(-....).
0290         ISUB=MINT(1)
0291         VINT(145)=VNT145
0292         VINT(146)=VNT146
0293         VINT(147)=VNT147
0294         IF(MSTP(82).LE.0) THEN
0295           XT2=0D0
0296         ELSEIF(MSTP(82).EQ.1) THEN
0297           XT2=XT2FAC*XT2/(XT2FAC-XT2*LOG(PYR(0)))
0298 C...Use with "Sudakov" for low b values when impact parameter dependence.
0299         ELSEIF(MSTP(82).EQ.2.OR.MINT(39).EQ.1) THEN
0300           IF(XT2.LT.1D0.AND.EXP(-XT2FAC*XT2/(VINT(149)*(XT2+
0301      &    VINT(149)))).GT.PYR(0)) XT2=1D0
0302           IF(XT2.GE.1D0) THEN
0303             XT2=(1D0+VINT(149))*XT2FAC/(XT2FAC-(1D0+VINT(149))*LOG(1D0-
0304      &      PYR(0)*(1D0-EXP(-XT2FAC/(VINT(149)*(1D0+VINT(149)))))))-
0305      &      VINT(149)
0306           ELSE
0307             XT2=-XT2FAC/LOG(EXP(-XT2FAC/(XT2+VINT(149)))+PYR(0)*
0308      &      (EXP(-XT2FAC/VINT(149))-EXP(-XT2FAC/(XT2+VINT(149)))))-
0309      &      VINT(149)
0310           ENDIF
0311           XT2=MAX(0.01D0*VINT(149),XT2)
0312 C...Use without "Sudakov" for high b values when impact parameter dep.
0313         ELSE
0314           XT2=(XC2+VINT(149))*(1D0+VINT(149))/(1D0+VINT(149)-
0315      &    PYR(0)*(1D0-XC2))-VINT(149)
0316           XT2=MAX(0.01D0*VINT(149),XT2)
0317         ENDIF
0318         VINT(25)=XT2
0319
0320 C...Low-pT: choose xT2, tau, y* and cos(theta-hat) fixed.
0321         IF(MSTP(82).LE.1.AND.XT2.LT.VINT(149)) THEN
0322           IF(MINT(82).EQ.1) NGEN(0,1)=NGEN(0,1)-MINT(143)
0323           IF(MINT(82).EQ.1) NGEN(ISUB,1)=NGEN(ISUB,1)-MINT(143)
0324           ISUB=95
0325           MINT(1)=ISUB
0326           VINT(21)=1D-12*VINT(149)
0327           VINT(22)=0D0
0328           VINT(23)=0D0
0329           VINT(25)=1D-12*VINT(149)
0330
0331         ELSE
0332 C...Multiple interactions (first semihard interaction).
0333 C...Choose tau and y*. Calculate cos(theta-hat).
0334           IF(PYR(0).LE.COEF(ISUB,1)) THEN
0335             TAUT=(2D0*(1D0+SQRT(1D0-XT2))/XT2-1D0)**PYR(0)
0336             TAU=XT2*(1D0+TAUT)**2/(4D0*TAUT)
0337           ELSE
0338             TAU=XT2*(1D0+TAN(PYR(0)*ATAN(SQRT(1D0/XT2-1D0)))**2)
0339           ENDIF
0340           VINT(21)=TAU
0341           CALL PYKLIM(2)
0342           RYST=PYR(0)
0343           MYST=1
0344           IF(RYST.GT.COEF(ISUB,8)) MYST=2
0345           IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
0346           CALL PYKMAP(2,MYST,PYR(0))
0347           VINT(23)=SQRT(MAX(0D0,1D0-XT2/TAU))*(-1)**INT(1.5D0+PYR(0))
0348         ENDIF
0349         VINT(71)=0.5D0*VINT(1)*SQRT(VINT(25))
0350
0351 C...Store results of cross-section calculation.
0352       ELSEIF(MMUL.EQ.4) THEN
0353         ISUB=MINT(1)
0354         VINT(145)=VNT145
0355         VINT(146)=VNT146
0356         VINT(147)=VNT147
0357         XTS=VINT(25)
0358         IF(ISET(ISUB).EQ.1) XTS=VINT(21)
0359         IF(ISET(ISUB).EQ.2)
0360      &  XTS=(4D0*VINT(48)+2D0*VINT(63)+2D0*VINT(64))/VINT(2)
0361         IF(ISET(ISUB).GE.3.AND.ISET(ISUB).LE.5) XTS=VINT(26)
0362         RBIN=MAX(0.000001D0,MIN(0.999999D0,XTS*(1D0+VINT(149))/
0363      &  (XTS+VINT(149))))
0364         IRBIN=INT(1D0+20D0*RBIN)
0365         IF(ISUB.EQ.96.AND.MSTP(171).EQ.0) THEN
0366           NMUL(IRBIN)=NMUL(IRBIN)+1
0367           SIGM(IRBIN)=SIGM(IRBIN)+VINT(153)
0368         ENDIF
0369
0370 C...Choose impact parameter if not already done.
0371       ELSEIF(MMUL.EQ.5) THEN
0372         ISUB=MINT(1)
0373         VINT(145)=VNT145
0374         VINT(146)=VNT146
0375         VINT(147)=VNT147
0376   150   IF(MINT(39).GT.0) THEN
0377         ELSEIF(MSTP(82).EQ.3) THEN
0378           EXPB2=PYR(0)
0379           B2=-LOG(PYR(0))
0380           VINT(148)=EXPB2/(PARU(2)*VNT147)
0381           VINT(139)=SQRT(B2)/BAVG
0382         ELSEIF(MSTP(82).EQ.4) THEN
0383           RTYPE=PYR(0)
0384           IF(RTYPE.LT.P83A) THEN
0385             B2=-LOG(PYR(0))
0386           ELSEIF(RTYPE.LT.P83A+P83B) THEN
0387             B2=-LOG(PYR(0))/CQ2R
0388           ELSE
0389             B2=-LOG(PYR(0))/CQ2I
0390           ENDIF
0391           VINT(148)=(P83A*EXP(-MIN(50D0,B2))+
0392      &    P83B*CQ2R*EXP(-MIN(50D0,B2*CQ2R))+
0393      &    P83C*CQ2I*EXP(-MIN(50D0,B2*CQ2I)))/(PARU(2)*VNT147)
0394           VINT(139)=SQRT(B2)/BAVG
0395         ELSEIF(PARP(83).GE.1.999D0) THEN
0396           POWIP=MAX(2D0,PARP(83))
0397           RPWIP=2D0/POWIP-1D0
0398           PROB1=POWIP/(2D0*EXP(-1D0)+POWIP)
0399   160     IF(PYR(0).LT.PROB1) THEN
0400             B2RPW=PYR(0)**(0.5D0*POWIP)
0401             ACCIP=EXP(-B2RPW)
0402           ELSE
0403             B2RPW=1D0-LOG(PYR(0))
0404             ACCIP=B2RPW**RPWIP
0405           ENDIF
0406           IF(ACCIP.LT.PYR(0)) GOTO 160
0407           VINT(148)=EXP(-B2RPW)/(PARU(2)*VNT147)
0408           VINT(139)=B2RPW**(1D0/POWIP)/BAVG
0409         ELSE
0410           POWIP=MAX(0.4D0,PARP(83))
0411           RPWIP=2D0/POWIP-1D0
0412           PROB1=RPWIP/(RPWIP+2D0**RPWIP*EXP(-RPWIP))
0413   170     IF(PYR(0).LT.PROB1) THEN
0414             B2RPW=2D0*RPWIP*PYR(0)
0415             ACCIP=(B2RPW/RPWIP)**RPWIP*EXP(RPWIP-B2RPW)
0416           ELSE
0417             B2RPW=2D0*(RPWIP-LOG(PYR(0)))
0418             ACCIP=(0.5D0*B2RPW/RPWIP)**RPWIP*EXP(RPWIP-0.5D0*B2RPW)
0419           ENDIF
0420           IF(ACCIP.LT .PYR(0)) GOTO 170
0421           VINT(148)=EXP(-B2RPW)/(PARU(2)*VNT147)
0422           VINT(139)=B2RPW**(1D0/POWIP)/BAVG
0423         ENDIF
0424
0425 C...Multiple interactions (variable impact parameter) : reject with
0426 C...probability exp(-overlap*cross-section above pT/normalization).
0427 C...Does not apply to low-b region, where "Sudakov" already included.
0428         VINT(150)=1D0
0429         IF(MINT(39).NE.1) THEN
0430           RNCOR=(IRBIN-20D0*RBIN)*NMUL(IRBIN)
0431           SIGCOR=(IRBIN-20D0*RBIN)*SIGM(IRBIN)
0432           DO 180 IBIN=IRBIN+1,20
0433             RNCOR=RNCOR+NMUL(IBIN)
0434             SIGCOR=SIGCOR+SIGM(IBIN)
0435   180     CONTINUE
0436           SIGABV=(SIGCOR/RNCOR)*VINT(149)*(1D0-XTS)/(XTS+VINT(149))
0437           IF(MSTP(171).EQ.1) SIGABV=SIGABV*VINT(2)/VINT(289)
0438           VINT(150)=EXP(-MIN(50D0,VNT146*VINT(148)*
0439      &    SIGABV/MAX(1D-10,SIGT(0,0,5))))
0440         ENDIF
0441         IF(MSTP(86).EQ.3.OR.(MSTP(86).EQ.2.AND.ISUB.NE.11.AND.
0442      &  ISUB.NE.12.AND.ISUB.NE.13.AND.ISUB.NE.28.AND.ISUB.NE.53
0443      &  .AND.ISUB.NE.68.AND.ISUB.NE.95.AND.ISUB.NE.96)) THEN
0444           IF(VINT(150).LT.PYR(0)) GOTO 150
0445           VINT(150)=1D0
0446         ENDIF
0447
0448 C...Generate additional multiple semihard interactions.
0449       ELSEIF(MMUL.EQ.6) THEN
0450
0451 C...Save data for hardest initeraction, to be restored.
0452         ISUBSV=MINT(1)
0453         VINT(145)=VNT145
0454         VINT(146)=VNT146
0455         VINT(147)=VNT147
0456         M13SV=MINT(13)
0457         M14SV=MINT(14)
0458         M15SV=MINT(15)
0459         M16SV=MINT(16)
0460         M21SV=MINT(21)
0461         M22SV=MINT(22)
0462         DO 190 J=11,80
0463           VINTSV(J)=VINT(J)
0464   190   CONTINUE
0465         V141SV=VINT(141)
0466         V142SV=VINT(142)
0467
0468 C...Store data on hardest interaction.
0469         XMI(1,1)=VINT(141)
0470         XMI(2,1)=VINT(142)
0471         PT2MI(1)=VINT(54)
0472         IMISEP(0)=MINT(84)
0473         IMISEP(1)=N
0474
0475 C...Change process to generate; sum of x values so far.
0476         ISUB=96
0477         MINT(1)=96
0478         VINT(143)=1D0-VINT(141)
0479         VINT(144)=1D0-VINT(142)
0480         VINT(151)=0D0
0481         VINT(152)=0D0
0482
0483 C...Initialize factors for PDF reshaping.
0484         DO 230 JS=1,2
0485           KFBEAM=MINT(10+JS)
0486           KFABM=IABS(KFBEAM)
0487           KFSBM=ISIGN(1,KFBEAM)
0488
0489 C...Zero flavour content of incoming beam particle.
0490           KFIVAL(JS,1)=0
0491           KFIVAL(JS,2)=0
0492           KFIVAL(JS,3)=0
0493 C...Flavour content of baryon.
0494           IF(KFABM.GT.1000) THEN
0495             KFIVAL(JS,1)=KFSBM*MOD(KFABM/1000,10)
0496             KFIVAL(JS,2)=KFSBM*MOD(KFABM/100,10)
0497             KFIVAL(JS,3)=KFSBM*MOD(KFABM/10,10)
0498 C...Flavour content of pi+-, K+-.
0499           ELSEIF(KFABM.EQ.211) THEN
0500             KFIVAL(JS,1)=KFSBM*2
0501             KFIVAL(JS,2)=-KFSBM
0502           ELSEIF(KFABM.EQ.321) THEN
0503             KFIVAL(JS,1)=-KFSBM*3
0504             KFIVAL(JS,2)=KFSBM*2
0505 C...Flavour content of pi0, gamma, K0S, K0L not defined yet.
0506           ENDIF
0507
0508 C...Zero initial valence and companion content.
0509           DO 200 IFL=-6,6
0510             NVC(JS,IFL)=0
0511   200     CONTINUE
0512
0513 C...Initiate listing of all incoming partons from two sides.
0514           NMI(JS)=0
0515           DO 210 I=MINT(84)+1,N
0516             IF(K(I,3).EQ.MINT(83)+2+JS) THEN
0517               IMI(JS,1,1)=I
0518               IMI(JS,1,2)=0
0519             ENDIF
0520   210     CONTINUE
0521
0522 C...Decide whether quarks in hard scattering were valence or sea.
0523           IFL=K(IMI(JS,1,1),2)
0524           IF (IABS(IFL).GT.6) GOTO 230
0525
0526 C...Get PDFs at X and Q2 of the parton shower initiator for the
0527 C...hard scattering.
0528           X=VINT(140+JS)
0529           IF(MSTP(61).GE.1) THEN
0530             Q2=PARP(62)**2
0531           ELSE
0532             Q2=VINT(54)
0533           ENDIF
0534 C...Note: XPSVC = x*pdf.
0535           MINT(30)=JS
0536           CALL PYPDFU(KFBEAM,X,Q2,XPQ)
0537           SEA=XPSVC(IFL,-1)
0538           VAL=XPSVC(IFL,0)
0539
0540 C...Decide (Extra factor x cancels in the division).
0541           RVCS=PYR(0)*(SEA+VAL)
0542           IVNOW=1
0543   220     IF (RVCS.LE.VAL.AND.IVNOW.GE.1) THEN
0544 C...Safety check that valence present; pi0/gamma/K0S/K0L special cases.
0545             IVNOW=0
0546             IF(KFIVAL(JS,1).EQ.IFL) IVNOW=IVNOW+1
0547             IF(KFIVAL(JS,2).EQ.IFL) IVNOW=IVNOW+1
0548             IF(KFIVAL(JS,3).EQ.IFL) IVNOW=IVNOW+1
0549             IF(KFIVAL(JS,1).EQ.0) THEN
0550               IF(KFBEAM.EQ.111.AND.IABS(IFL).LE.2) IVNOW=1
0551               IF(KFBEAM.EQ.22.AND.IABS(IFL).LE.5) IVNOW=1
0552               IF((KFBEAM.EQ.130.OR.KFBEAM.EQ.310).AND.
0553      &        (IABS(IFL).EQ.1.OR.IABS(IFL).EQ.3)) IVNOW=1
0554             ENDIF
0555             IF(IVNOW.EQ.0) GOTO 220
0556 C...Mark valence.
0557             IMI(JS,1,2)=0
0558 C...Sets valence content of gamma, pi0, K0S, K0L if not done.
0559             IF(KFIVAL(JS,1).EQ.0) THEN
0560               IF(KFBEAM.EQ.111.OR.KFBEAM.EQ.22) THEN
0561                 KFIVAL(JS,1)=IFL
0562                 KFIVAL(JS,2)=-IFL
0563               ELSEIF(KFBEAM.EQ.130.OR.KFBEAM.EQ.310) THEN
0564                 KFIVAL(JS,1)=IFL
0565                 IF(IABS(IFL).EQ.1) KFIVAL(JS,2)=ISIGN(3,-IFL)
0566                 IF(IABS(IFL).NE.1) KFIVAL(JS,2)=ISIGN(1,-IFL)
0567               ENDIF
0568             ENDIF
0569
0570 C...If sea, add opposite sign companion parton. Store X and I.
0571           ELSE
0572             NVC(JS,-IFL)=NVC(JS,-IFL)+1
0573             XASSOC(JS,-IFL,NVC(JS,-IFL))=X
0574 C...Set pointer to companion
0575             IMI(JS,1,2)=-NVC(JS,-IFL)
0576           ENDIF
0577   230   CONTINUE
0578
0579 C...Update counter number of multiple interactions.
0580         NMI(1)=1
0581         NMI(2)=1
0582
0583 C...Set up starting values for iteration in xT2.
0584         IF(MSTP(86).EQ.3.OR.(MSTP(86).EQ.2.AND.ISUBSV.NE.11.AND.
0585      &  ISUBSV.NE.12.AND.ISUBSV.NE.13.AND.ISUBSV.NE.28.AND.
0586      &  ISUBSV.NE.53.AND.ISUBSV.NE.68.AND.ISUBSV.NE.95.AND.
0587      &  ISUBSV.NE.96)) THEN
0588           XT2=(1D0-VINT(141))*(1D0-VINT(142))
0589         ELSE
0590           XT2=VINT(25)
0591           IF(ISET(ISUBSV).EQ.1) XT2=VINT(21)
0592           IF(ISET(ISUBSV).EQ.2)
0593      &    XT2=(4D0*VINT(48)+2D0*VINT(63)+2D0*VINT(64))/VINT(2)
0594           IF(ISET(ISUBSV).GE.3.AND.ISET(ISUBSV).LE.5) XT2=VINT(26)
0595         ENDIF
0596         IF(MSTP(82).LE.1) THEN
0597           SIGRAT=XSEC(ISUB,1)/MAX(1D-10,VINT(315)*VINT(316)*SIGT(0,0,5))
0598           IF(MINT(141).NE.0.OR.MINT(142).NE.0) SIGRAT=SIGRAT*
0599      &    VINT(317)/(VINT(318)*VINT(320))
0600           XT2FAC=SIGRAT*VINT(149)/(1D0-VINT(149))
0601         ELSE
0602           XT2FAC=VNT146*VINT(148)*XSEC(ISUB,1)/
0603      &    MAX(1D-10,SIGT(0,0,5))*VINT(149)*(1D0+VINT(149))
0604         ENDIF
0605         VINT(63)=0D0
0606         VINT(64)=0D0
0607
0608 C...Iterate downwards in xT2.
0609   240   IF((MINT(35).EQ.2.AND.MSTP(81).EQ.10).OR.ISUBSV.EQ.95) THEN
0610           XT2=0D0
0611           GOTO 440
0612         ELSEIF(MSTP(82).LE.1) THEN
0613           XT2=XT2FAC*XT2/(XT2FAC-XT2*LOG(PYR(0)))
0614           IF(XT2.LT.VINT(149)) GOTO 440
0615         ELSE
0616           IF(XT2.LE.0.01001D0*VINT(149)) GOTO 440
0617           XT2=XT2FAC*(XT2+VINT(149))/(XT2FAC-(XT2+VINT(149))*
0618      &    LOG(PYR(0)))-VINT(149)
0619           IF(XT2.LE.0D0) GOTO 440
0620           XT2=MAX(0.01D0*VINT(149),XT2)
0621         ENDIF
0622         VINT(25)=XT2
0623
0624 C...Choose tau and y*. Calculate cos(theta-hat).
0625         IF(PYR(0).LE.COEF(ISUB,1)) THEN
0626           TAUT=(2D0*(1D0+SQRT(1D0-XT2))/XT2-1D0)**PYR(0)
0627           TAU=XT2*(1D0+TAUT)**2/(4D0*TAUT)
0628         ELSE
0629           TAU=XT2*(1D0+TAN(PYR(0)*ATAN(SQRT(1D0/XT2-1D0)))**2)
0630         ENDIF
0631         VINT(21)=TAU
0632 C...New: require shat > 1.
0633         IF(TAU*VINT(2).LT.1D0) GOTO 240
0634         CALL PYKLIM(2)
0635         RYST=PYR(0)
0636         MYST=1
0637         IF(RYST.GT.COEF(ISUB,8)) MYST=2
0638         IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
0639         CALL PYKMAP(2,MYST,PYR(0))
0640         VINT(23)=SQRT(MAX(0D0,1D0-XT2/TAU))*(-1)**INT(1.5D0+PYR(0))
0641
0642 C...Check that x not used up. Accept or reject kinematical variables.
0643         X1M=SQRT(TAU)*EXP(VINT(22))
0644         X2M=SQRT(TAU)*EXP(-VINT(22))
0645         IF(VINT(143)-X1M.LT.0.01D0.OR.VINT(144)-X2M.LT.0.01D0) GOTO 240
0646         VINT(71)=0.5D0*VINT(1)*SQRT(XT2)
0647         CALL PYSIGH(NCHN,SIGS)
0648         IF(MINT(141).NE.0.OR.MINT(142).NE.0) SIGS=SIGS*VINT(320)
0649         IF(SIGS.LT.XSEC(ISUB,1)*PYR(0)) GOTO 240
0650         IF(MINT(141).NE.0.OR.MINT(142).NE.0) SIGS=SIGS/VINT(320)
0651
0652 C...Reset K, P and V vectors.
0653         DO 260 I=N+1,N+4
0654           DO 250 J=1,5
0655             K(I,J)=0
0656             P(I,J)=0D0
0657             V(I,J)=0D0
0658   250     CONTINUE
0659   260   CONTINUE
0660         PT=0.5D0*VINT(1)*SQRT(XT2)
0661
0662 C...Choose flavour of reacting partons (and subprocess).
0663         RSIGS=SIGS*PYR(0)
0664         DO 270 ICHN=1,NCHN
0665           KFL1=ISIG(ICHN,1)
0666           KFL2=ISIG(ICHN,2)
0667           ICONMI=ISIG(ICHN,3)
0668           RSIGS=RSIGS-SIGH(ICHN)
0669           IF(RSIGS.LE.0D0) GOTO 280
0670   270   CONTINUE
0671
0672 C...Reassign to appropriate process codes.
0673   280   ISUBMI=ICONMI/10
0674         ICONMI=MOD(ICONMI,10)
0675
0676 C...Choose new quark flavour for annihilation graphs
0677         IF(ISUBMI.EQ.12.OR.ISUBMI.EQ.53) THEN
0678           SH=TAU*VINT(2)
0679           CALL PYWIDT(21,SH,WDTP,WDTE)
0680   290     RKFL=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))*PYR(0)
0681           DO 300 I=1,MDCY(21,3)
0682             KFLF=KFDP(I+MDCY(21,2)-1,1)
0683             RKFL=RKFL-(WDTE(I,1)+WDTE(I,2)+WDTE(I,4))
0684             IF(RKFL.LE.0D0) GOTO 310
0685   300     CONTINUE
0686   310     IF(ISUBMI.EQ.53.AND.ICONMI.LE.2) THEN
0687             IF(KFLF.GE.4) GOTO 290
0688           ELSEIF(ISUBMI.EQ.53.AND.ICONMI.LE.4) THEN
0689             KFLF=4
0690             ICONMI=ICONMI-2
0691           ELSEIF(ISUBMI.EQ.53) THEN
0692             KFLF=5
0693             ICONMI=ICONMI-4
0694           ENDIF
0695         ENDIF
0696
0697 C...Final state flavours and colour flow: default values
0698         JS=1
0699         KFL3=KFL1
0700         KFL4=KFL2
0701         KCC=20
0702         KCS=ISIGN(1,KFL1)
0703
0704         IF(ISUBMI.EQ.11) THEN
0705 C...f + f' -> f + f' (g exchange); th = (p(f)-p(f))**2
0706           KCC=ICONMI
0707           IF(KFL1*KFL2.LT.0) KCC=KCC+2
0708
0709         ELSEIF(ISUBMI.EQ.12) THEN
0710 C...f + fbar -> f' + fbar'; th = (p(f)-p(f'))**2
0711           KFL3=ISIGN(KFLF,KFL1)
0712           KFL4=-KFL3
0713           KCC=4
0714
0715         ELSEIF(ISUBMI.EQ.13) THEN
0716 C...f + fbar -> g + g; th arbitrary
0717           KFL3=21
0718           KFL4=21
0719           KCC=ICONMI+4
0720
0721         ELSEIF(ISUBMI.EQ.28) THEN
0722 C...f + g -> f + g; th = (p(f)-p(f))**2
0723           IF(KFL1.EQ.21) JS=2
0724           KCC=ICONMI+6
0725           IF(KFL1.EQ.21) KCC=KCC+2
0726           IF(KFL1.NE.21) KCS=ISIGN(1,KFL1)
0727           IF(KFL2.NE.21) KCS=ISIGN(1,KFL2)
0728
0729         ELSEIF(ISUBMI.EQ.53) THEN
0730 C...g + g -> f + fbar; th arbitrary
0731           KCS=(-1)**INT(1.5D0+PYR(0))
0732           KFL3=ISIGN(KFLF,KCS)
0733           KFL4=-KFL3
0734           KCC=ICONMI+10
0735
0736         ELSEIF(ISUBMI.EQ.68) THEN
0737 C...g + g -> g + g; th arbitrary
0738           KCC=ICONMI+12
0739           KCS=(-1)**INT(1.5D0+PYR(0))
0740         ENDIF
0741
0742 C...Store flavours of scattering.
0743         MINT(13)=KFL1
0744         MINT(14)=KFL2
0745         MINT(15)=KFL1
0746         MINT(16)=KFL2
0747         MINT(21)=KFL3
0748         MINT(22)=KFL4
0749
0750 C...Set flavours and mothers of scattering partons.
0751         K(N+1,1)=14
0752         K(N+2,1)=14
0753         K(N+3,1)=3
0754         K(N+4,1)=3
0755         K(N+1,2)=KFL1
0756         K(N+2,2)=KFL2
0757         K(N+3,2)=KFL3
0758         K(N+4,2)=KFL4
0759         K(N+1,3)=MINT(83)+1
0760         K(N+2,3)=MINT(83)+2
0761         K(N+3,3)=N+1
0762         K(N+4,3)=N+2
0763
0764 C...Store colour connection indices.
0765         DO 320 J=1,2
0766           JC=J
0767           IF(KCS.EQ.-1) JC=3-J
0768           IF(ICOL(KCC,1,JC).NE.0) K(N+1,J+3)=N+ICOL(KCC,1,JC)
0769           IF(ICOL(KCC,2,JC).NE.0) K(N+2,J+3)=N+ICOL(KCC,2,JC)
0770           IF(ICOL(KCC,3,JC).NE.0) K(N+3,J+3)=MSTU(5)*(N+ICOL(KCC,3,JC))
0771           IF(ICOL(KCC,4,JC).NE.0) K(N+4,J+3)=MSTU(5)*(N+ICOL(KCC,4,JC))
0772   320   CONTINUE
0773
0774 C...Store incoming and outgoing partons in their CM-frame.
0775         SHR=SQRT(TAU)*VINT(1)
0776         P(N+1,3)=0.5D0*SHR
0777         P(N+1,4)=0.5D0*SHR
0778         P(N+2,3)=-0.5D0*SHR
0779         P(N+2,4)=0.5D0*SHR
0780         P(N+3,5)=PYMASS(K(N+3,2))
0781         P(N+4,5)=PYMASS(K(N+4,2))
0782         IF(P(N+3,5)+P(N+4,5).GE.SHR) GOTO 240
0783         P(N+3,4)=0.5D0*(SHR+(P(N+3,5)**2-P(N+4,5)**2)/SHR)
0784         P(N+3,3)=SQRT(MAX(0D0,P(N+3,4)**2-P(N+3,5)**2))
0785         P(N+4,4)=SHR-P(N+3,4)
0786         P(N+4,3)=-P(N+3,3)
0787
0788 C...Rotate outgoing partons using cos(theta)=(th-uh)/lam(sh,sqm3,sqm4)
0789         PHI=PARU(2)*PYR(0)
0790         CALL PYROBO(N+3,N+4,ACOS(VINT(23)),PHI,0D0,0D0,0D0)
0791
0792 C...Set up default values before showers.
0793         MINT(31)=MINT(31)+1
0794         IPU1=N+1
0795         IPU2=N+2
0796         IPU3=N+3
0797         IPU4=N+4
0798         VINT(141)=VINT(41)
0799         VINT(142)=VINT(42)
0800         N=N+4
0801
0802 C...Showering of initial state partons (optional).
0803 C...Note: no showering of final state partons here; it comes later.
0804         IF(MSTP(84).GE.1.AND.MSTP(61).GE.1) THEN
0805           MINT(51)=0
0806           ALAMSV=PARJ(81)
0807           PARJ(81)=PARP(72)
0808           NSAV=N
0809           DO 340 I=1,4
0810             DO 330 J=1,5
0811               KSAV(I,J)=K(N-4+I,J)
0812               PSAV(I,J)=P(N-4+I,J)
0813   330       CONTINUE
0814   340     CONTINUE
0815           CALL PYSSPA(IPU1,IPU2)
0816           PARJ(81)=ALAMSV
0817 C...If shower failed then restore to situation before shower.
0818           IF(MINT(51).GE.1) THEN
0819             N=NSAV
0820             DO 360 I=1,4
0821               DO 350 J=1,5
0822                 K(N-4+I,J)=KSAV(I,J)
0823                 P(N-4+I,J)=PSAV(I,J)
0824   350         CONTINUE
0825   360       CONTINUE
0826             IPU1=N-3
0827             IPU2=N-2
0828             VINT(141)=VINT(41)
0829             VINT(142)=VINT(42)
0830           ENDIF
0831         ENDIF
0832
0833 C...Keep track of loose colour ends and information on scattering.
0834   370   IMI(1,MINT(31),1)=IPU1
0835         IMI(2,MINT(31),1)=IPU2
0836         IMI(1,MINT(31),2)=0
0837         IMI(2,MINT(31),2)=0
0838         XMI(1,MINT(31))=VINT(141)
0839         XMI(2,MINT(31))=VINT(142)
0840         PT2MI(MINT(31))=VINT(54)
0841         IMISEP(MINT(31))=N
0842
0843 C...Decide whether quarks in last scattering were valence, companion or
0844 C...sea.
0845         DO 430 JS=1,2
0846           KFBEAM=MINT(10+JS)
0847           KFSBM=ISIGN(1,MINT(10+JS))
0848           IFL=K(IMI(JS,MINT(31),1),2)
0849           IMI(JS,MINT(31),2)=0
0850           IF (IABS(IFL).GT.6) GOTO 430
0851
0852 C...Get PDFs at X and Q2 of the parton shower initiator for the
0853 C...last scattering. At this point VINT(143:144) do not yet
0854 C...include the scattered x values VINT(141:142).
0855           X=VINT(140+JS)/VINT(142+JS)
0856           IF(MSTP(84).GE.1.AND.MSTP(61).GE.1) THEN
0857             Q2=PARP(62)**2
0858           ELSE
0859             Q2=VINT(54)
0860           ENDIF
0861 C...Note: XPSVC = x*pdf.
0862           MINT(30)=JS
0863           CALL PYPDFU(KFBEAM,X,Q2,XPQ)
0864           SEA=XPSVC(IFL,-1)
0865           VAL=XPSVC(IFL,0)
0866           CMP=0D0
0867           DO 380 IVC=1,NVC(JS,IFL)
0868             CMP=CMP+XPSVC(IFL,IVC)
0869   380     CONTINUE
0870
0871 C...Decide (Extra factor x cancels in the dvision).
0872           RVCS=PYR(0)*(SEA+VAL+CMP)
0873           IVNOW=1
0874   390     IF (RVCS.LE.VAL.AND.IVNOW.GE.1) THEN
0875 C...Safety check that valence present; pi0/gamma/K0S/K0L special cases.
0876             IVNOW=0
0877             IF(KFIVAL(JS,1).EQ.IFL) IVNOW=IVNOW+1
0878             IF(KFIVAL(JS,2).EQ.IFL) IVNOW=IVNOW+1
0879             IF(KFIVAL(JS,3).EQ.IFL) IVNOW=IVNOW+1
0880             IF(KFIVAL(JS,1).EQ.0) THEN
0881               IF(KFBEAM.EQ.111.AND.IABS(IFL).LE.2) IVNOW=1
0882               IF(KFBEAM.EQ.22.AND.IABS(IFL).LE.5) IVNOW=1
0883               IF((KFBEAM.EQ.130.OR.KFBEAM.EQ.310).AND.
0884      &        (IABS(IFL).EQ.1.OR.IABS(IFL).EQ.3)) IVNOW=1
0885             ELSE
0886               DO 400 I1=1,NMI(JS)
0887                 IF (K(IMI(JS,I1,1),2).EQ.IFL.AND.IMI(JS,I1,2).EQ.0)
0888      &            IVNOW=IVNOW-1
0889   400         CONTINUE
0890             ENDIF
0891             IF(IVNOW.EQ.0) GOTO 390
0892 C...Mark valence.
0893             IMI(JS,MINT(31),2)=0
0894 C...Sets valence content of gamma, pi0, K0S, K0L if not done.
0895             IF(KFIVAL(JS,1).EQ.0) THEN
0896               IF(KFBEAM.EQ.111.OR.KFBEAM.EQ.22) THEN
0897                 KFIVAL(JS,1)=IFL
0898                 KFIVAL(JS,2)=-IFL
0899               ELSEIF(KFBEAM.EQ.130.OR.KFBEAM.EQ.310) THEN
0900                 KFIVAL(JS,1)=IFL
0901                 IF(IABS(IFL).EQ.1) KFIVAL(JS,2)=ISIGN(3,-IFL)
0902                 IF(IABS(IFL).NE.1) KFIVAL(JS,2)=ISIGN(1,-IFL)
0903               ENDIF
0904             ENDIF
0905
0906           ELSEIF (RVCS.LE.VAL+SEA.OR.NVC(JS,IFL).EQ.0) THEN
0907 C...If sea, add opposite sign companion parton. Store X and I.
0908             NVC(JS,-IFL)=NVC(JS,-IFL)+1
0909             XASSOC(JS,-IFL,NVC(JS,-IFL))=X
0910 C...Set pointer to companion
0911             IMI(JS,MINT(31),2)=-NVC(JS,-IFL)
0912           ELSE
0913 C...If companion, decide which one.
0914             CMPSUM=VAL+SEA
0915             ISEL=0
0916   410       ISEL=ISEL+1
0917             CMPSUM=CMPSUM+XPSVC(IFL,ISEL)
0918             IF (RVCS.GT.CMPSUM.AND.ISEL.LT.NVC(JS,IFL)) GOTO 410
0919 C...Find original sea (anti-)quark:
0920             IASSOC=0
0921             DO 420 I1=1,NMI(JS)
0922               IF (K(IMI(JS,I1,1),2).NE.-IFL) GOTO 420
0923               IF (-IMI(JS,I1,2).EQ.ISEL) THEN
0924                 IMI(JS,MINT(31),2)=IMI(JS,I1,1)
0925                 IMI(JS,I1,2)=IMI(JS,MINT(31),1)
0926               ENDIF
0927   420       CONTINUE
0928 C...Change X to what associated companion had, so that the correct
0929 C...amount of momentum can be subtracted from the companion sum below.
0930             X=XASSOC(JS,IFL,ISEL)
0931 C...Mark companion read.
0932             XASSOC(JS,IFL,ISEL)=0D0
0933           ENDIF
0934  430    CONTINUE
0935
0936 C...Global statistics.
0937         MINT(351)=MINT(351)+1
0938         VINT(351)=VINT(351)+PT
0939         IF (MINT(351).EQ.1) VINT(356)=PT
0940
0941 C...Update remaining energy and other counters.
0942         IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
0943           CALL PYERRM(11,'(PYMIGN:) no more memory left in PYJETS')
0944           MINT(51)=1
0945           RETURN
0946         ENDIF
0947         NMI(1)=NMI(1)+1
0948         NMI(2)=NMI(2)+1
0949         VINT(151)=VINT(151)+VINT(41)
0950         VINT(152)=VINT(152)+VINT(42)
0951         VINT(143)=VINT(143)-VINT(141)
0952         VINT(144)=VINT(144)-VINT(142)
0953
0954 C...Iterate, with more interactions allowed.
0955         IF(MINT(31).LT.240) GOTO 240
0956  440    CONTINUE
0957
0958 C...Restore saved quantities for hardest interaction.
0959         MINT(1)=ISUBSV
0960         MINT(13)=M13SV
0961         MINT(14)=M14SV
0962         MINT(15)=M15SV
0963         MINT(16)=M16SV
0964         MINT(21)=M21SV
0965         MINT(22)=M22SV
0966         DO 450 J=11,80
0967           VINT(J)=VINTSV(J)
0968   450   CONTINUE
0969         VINT(141)=V141SV
0970         VINT(142)=V142SV
0971
0972       ENDIF
0973
0974 C...Format statements for printout.
0975  5000 FORMAT(/1X,'****** PYMIGN: initialization of multiple inter',
0976      &'actions for MSTP(82) =',I2,' ******')
0977  5100 FORMAT(8X,'pT0 =',F5.2,' GeV gives sigma(parton-parton) =',1P,
0978      &D9.2,' mb: rejected')
0979  5200 FORMAT(8X,'pT0 =',F5.2,' GeV gives sigma(parton-parton) =',1P,
0980      &D9.2,' mb: accepted')
0981
0982       RETURN
0983       END