pythiaeRHIC/pythia/pymult.f

0001
0002 C*********************************************************************
0003
0004 C...PYMULT
0005 C...Initializes treatment of multiple interactions, selects kinematics
0006 C...of hardest interaction if low-pT physics included in run, and
0007 C...generates all non-hardest interactions.
0008
0009       SUBROUTINE PYMULT(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 C...Commonblocks.
0016       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0017       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0018       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0019       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0020       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0021       COMMON/PYINT1/MINT(400),VINT(400)
0022       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0023       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
0024       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
0025       COMMON/PYINT7/SIGT(0:6,0:6,0:5)
0026       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYSUBS/,/PYPARS/,/PYINT1/,
0027      &/PYINT2/,/PYINT3/,/PYINT5/,/PYINT7/
0028 C...Local arrays and saved variables.
0029       DIMENSION NMUL(20),SIGM(20),KSTR(500,2),VINTSV(80)
0030       SAVE XT2,XT2FAC,XC2,XTS,IRBIN,RBIN,NMUL,SIGM,P83A,P83B,P83C,
0031      &CQ2I,CQ2R,PIK,BDIV,B,PLOWB,PHIGHB,PALLB,S4A,S4B,S4C,POWIP,
0032      &RPWIP,B2RPDV,B2RPMX,BAVG,VNT145,VNT146,VNT147
0033
0034 C...Initialization of multiple interaction treatment.
0035       IF(MMUL.EQ.1) THEN
0036         IF(MSTP(122).GE.1) WRITE(MSTU(11),5000) MSTP(82)
0037         ISUB=96
0038         MINT(1)=96
0039         VINT(63)=0D0
0040         VINT(64)=0D0
0041         VINT(143)=1D0
0042         VINT(144)=1D0
0043
0044 C...Loop over phase space points: xT2 choice in 20 bins.
0045   100   SIGSUM=0D0
0046         DO 120 IXT2=1,20
0047           NMUL(IXT2)=MSTP(83)
0048           SIGM(IXT2)=0D0
0049           DO 110 ITRY=1,MSTP(83)
0050             RSCA=0.05D0*((21-IXT2)-PYR(0))
0051             XT2=VINT(149)*(1D0+VINT(149))/(VINT(149)+RSCA)-VINT(149)
0052             XT2=MAX(0.01D0*VINT(149),XT2)
0053             VINT(25)=XT2
0054
0055 C...Choose tau and y*. Calculate cos(theta-hat).
0056             IF(PYR(0).LE.COEF(ISUB,1)) THEN
0057               TAUT=(2D0*(1D0+SQRT(1D0-XT2))/XT2-1D0)**PYR(0)
0058               TAU=XT2*(1D0+TAUT)**2/(4D0*TAUT)
0059             ELSE
0060               TAU=XT2*(1D0+TAN(PYR(0)*ATAN(SQRT(1D0/XT2-1D0)))**2)
0061             ENDIF
0062             VINT(21)=TAU
0063             CALL PYKLIM(2)
0064             RYST=PYR(0)
0065             MYST=1
0066             IF(RYST.GT.COEF(ISUB,8)) MYST=2
0067             IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
0068             CALL PYKMAP(2,MYST,PYR(0))
0069             VINT(23)=SQRT(MAX(0D0,1D0-XT2/TAU))*(-1)**INT(1.5D0+PYR(0))
0070
0071 C...Calculate differential cross-section.
0072             VINT(71)=0.5D0*VINT(1)*SQRT(XT2)
0073             CALL PYSIGH(NCHN,SIGS)
0074             SIGM(IXT2)=SIGM(IXT2)+SIGS
0075   110     CONTINUE
0076           SIGSUM=SIGSUM+SIGM(IXT2)
0077   120   CONTINUE
0078         SIGSUM=SIGSUM/(20D0*MSTP(83))
0079
0080 C...Reject result if sigma(parton-parton) is smaller than hadronic one.
0081         IF(SIGSUM.LT.1.1D0*SIGT(0,0,5)) THEN
0082           IF(MSTP(122).GE.1) WRITE(MSTU(11),5100)
0083      &    PARP(82)*(VINT(1)/PARP(89))**PARP(90),SIGSUM
0084           PARP(82)=0.9D0*PARP(82)
0085           VINT(149)=4D0*(PARP(82)*(VINT(1)/PARP(89))**PARP(90))**2/
0086      &    VINT(2)
0087           GOTO 100
0088         ENDIF
0089         IF(MSTP(122).GE.1) WRITE(MSTU(11),5200)
0090      &  PARP(82)*(VINT(1)/PARP(89))**PARP(90), SIGSUM
0091
0092 C...Start iteration to find k factor.
0093         YKE=SIGSUM/MAX(1D-10,SIGT(0,0,5))
0094         P83A=(1D0-PARP(83))**2
0095         P83B=2D0*PARP(83)*(1D0-PARP(83))
0096         P83C=PARP(83)**2
0097         CQ2I=1D0/PARP(84)**2
0098         CQ2R=2D0/(1D0+PARP(84)**2)
0099         SO=0.5D0
0100         XI=0D0
0101         YI=0D0
0102         XF=0D0
0103         YF=0D0
0104         XK=0.5D0
0105         IIT=0
0106   130   IF(IIT.EQ.0) THEN
0107           XK=2D0*XK
0108         ELSEIF(IIT.EQ.1) THEN
0109           XK=0.5D0*XK
0110         ELSE
0111           XK=XI+(YKE-YI)*(XF-XI)/(YF-YI)
0112         ENDIF
0113
0114 C...Evaluate overlap integrals. Find where to divide the b range.
0115         IF(MSTP(82).EQ.2) THEN
0116           SP=0.5D0*PARU(1)*(1D0-EXP(-XK))
0117           SOP=SP/PARU(1)
0118         ELSE
0119           IF(MSTP(82).EQ.3) THEN
0120             DELTAB=0.02D0
0121           ELSEIF(MSTP(82).EQ.4) THEN
0122             DELTAB=MIN(0.01D0,0.05D0*PARP(84))
0123           ELSE
0124             POWIP=MAX(0.4D0,PARP(83))
0125             RPWIP=2D0/POWIP-1D0
0126             DELTAB=MAX(0.02D0,0.02D0*(2D0/POWIP)**(1D0/POWIP))
0127             SO=0D0
0128           ENDIF
0129           SP=0D0
0130           SOP=0D0
0131           BSP=0D0
0132           SOHIGH=0D0
0133           IBDIV=0
0134           B=-0.5D0*DELTAB
0135   140     B=B+DELTAB
0136           IF(MSTP(82).EQ.3) THEN
0137             OV=EXP(-B**2)/PARU(2)
0138           ELSEIF(MSTP(82).EQ.4) THEN
0139             OV=(P83A*EXP(-MIN(50D0,B**2))+
0140      &      P83B*CQ2R*EXP(-MIN(50D0,B**2*CQ2R))+
0141      &      P83C*CQ2I*EXP(-MIN(50D0,B**2*CQ2I)))/PARU(2)
0142           ELSE
0143             OV=EXP(-B**POWIP)/PARU(2)
0144             SO=SO+PARU(2)*B*DELTAB*OV
0145           ENDIF
0146           IF(IBDIV.EQ.1) SOHIGH=SOHIGH+PARU(2)*B*DELTAB*OV
0147           PACC=1D0-EXP(-MIN(50D0,PARU(1)*XK*OV))
0148           SP=SP+PARU(2)*B*DELTAB*PACC
0149           SOP=SOP+PARU(2)*B*DELTAB*OV*PACC
0150           BSP=BSP+B*PARU(2)*B*DELTAB*PACC
0151           IF(IBDIV.EQ.0.AND.PARU(1)*XK*OV.LT.1D0) THEN
0152             IBDIV=1
0153             BDIV=B+0.5D0*DELTAB
0154           ENDIF
0155           IF(B.LT.1D0.OR.B*PACC.GT.1D-6) GOTO 140
0156         ENDIF
0157         YK=PARU(1)*XK*SO/SP
0158
0159 C...Continue iteration until convergence.
0160         IF(YK.LT.YKE) THEN
0161           XI=XK
0162           YI=YK
0163           IF(IIT.EQ.1) IIT=2
0164         ELSE
0165           XF=XK
0166           YF=YK
0167           IF(IIT.EQ.0) IIT=1
0168         ENDIF
0169         IF(ABS(YK-YKE).GE.1D-5*YKE) GOTO 130
0170
0171 C...Store some results for subsequent use.
0172         BAVG=BSP/SP
0173         VINT(145)=SIGSUM
0174         VINT(146)=SOP/SO
0175         VINT(147)=SOP/SP
0176         VNT145=VINT(145)
0177         VNT146=VINT(146)
0178         VNT147=VINT(147)
0179 C...PIK = PARU(1)*XK = (VINT(146)/VINT(147))*sigma_jet/sigma_nondiffr.
0180         PIK=(VNT146/VNT147)*YKE
0181
0182 C...Find relative weight for low and high impact parameter.
0183       PLOWB=PARU(1)*BDIV**2
0184       IF(MSTP(82).EQ.3) THEN
0185         PHIGHB=PIK*0.5*EXP(-BDIV**2)
0186       ELSEIF(MSTP(82).EQ.4) THEN
0187         S4A=P83A*EXP(-BDIV**2)
0188         S4B=P83B*EXP(-BDIV**2*CQ2R)
0189         S4C=P83C*EXP(-BDIV**2*CQ2I)
0190         PHIGHB=PIK*0.5*(S4A+S4B+S4C)
0191       ELSEIF(PARP(83).GE.1.999D0) THEN
0192         PHIGHB=PIK*SOHIGH
0193         B2RPDV=BDIV**POWIP
0194       ELSE
0195         PHIGHB=PIK*SOHIGH
0196         B2RPDV=BDIV**POWIP
0197         B2RPMX=MAX(2D0*RPWIP,B2RPDV)
0198       ENDIF
0199       PALLB=PLOWB+PHIGHB
0200
0201 C...Initialize iteration in xT2 for hardest interaction.
0202       ELSEIF(MMUL.EQ.2) THEN
0203         VINT(145)=VNT145
0204         VINT(146)=VNT146
0205         VINT(147)=VNT147
0206         IF(MSTP(82).LE.0) THEN
0207         ELSEIF(MSTP(82).EQ.1) THEN
0208           XT2=1D0
0209           SIGRAT=XSEC(96,1)/MAX(1D-10,VINT(315)*VINT(316)*SIGT(0,0,5))
0210           IF(MINT(141).NE.0.OR.MINT(142).NE.0) SIGRAT=SIGRAT*
0211      &    VINT(317)/(VINT(318)*VINT(320))
0212           XT2FAC=SIGRAT*VINT(149)/(1D0-VINT(149))
0213         ELSEIF(MSTP(82).EQ.2) THEN
0214           XT2=1D0
0215           XT2FAC=VNT146*XSEC(96,1)/MAX(1D-10,SIGT(0,0,5))*
0216      &    VINT(149)*(1D0+VINT(149))
0217         ELSE
0218           XC2=4D0*CKIN(3)**2/VINT(2)
0219           IF(CKIN(3).LE.CKIN(5).OR.MINT(82).GE.2) XC2=0D0
0220         ENDIF
0221
0222 C...Select impact parameter for hardest interaction.
0223         IF(MSTP(82).LE.2) RETURN
0224   142   IF(PYR(0)*PALLB.LT.PLOWB) THEN
0225 C...Treatment in low b region.
0226           MINT(39)=1
0227           B=BDIV*SQRT(PYR(0))
0228           IF(MSTP(82).EQ.3) THEN
0229             OV=EXP(-B**2)/PARU(2)
0230           ELSEIF(MSTP(82).EQ.4) THEN
0231             OV=(P83A*EXP(-MIN(50D0,B**2))+
0232      &      P83B*CQ2R*EXP(-MIN(50D0,B**2*CQ2R))+
0233      &      P83C*CQ2I*EXP(-MIN(50D0,B**2*CQ2I)))/PARU(2)
0234           ELSE
0235             OV=EXP(-B**POWIP)/PARU(2)
0236           ENDIF
0237           VINT(148)=OV/VNT147
0238           PACC=1D0-EXP(-MIN(50D0,PIK*OV))
0239           XT2=1D0
0240           XT2FAC=VNT146*VINT(148)*XSEC(96,1)/MAX(1D-10,SIGT(0,0,5))*
0241      &    VINT(149)*(1D0+VINT(149))
0242         ELSE
0243 C...Treatment in high b region.
0244           MINT(39)=2
0245           IF(MSTP(82).EQ.3) THEN
0246             B=SQRT(BDIV**2-LOG(PYR(0)))
0247             OV=EXP(-B**2)/PARU(2)
0248           ELSEIF(MSTP(82).EQ.4) THEN
0249             S4RNDM=PYR(0)*(S4A+S4B+S4C)
0250             IF(S4RNDM.LT.S4A) THEN
0251               B=SQRT(BDIV**2-LOG(PYR(0)))
0252             ELSEIF(S4RNDM.LT.S4A+S4B) THEN
0253               B=SQRT(BDIV**2-LOG(PYR(0))/CQ2R)
0254             ELSE
0255               B=SQRT(BDIV**2-LOG(PYR(0))/CQ2I)
0256             ENDIF
0257             OV=(P83A*EXP(-MIN(50D0,B**2))+
0258      &      P83B*CQ2R*EXP(-MIN(50D0,B**2*CQ2R))+
0259      &      P83C*CQ2I*EXP(-MIN(50D0,B**2*CQ2I)))/PARU(2)
0260           ELSEIF(PARP(83).GE.1.999D0) THEN
0261   144       B2RPW=B2RPDV-LOG(PYR(0))
0262             ACCIP=(B2RPW/B2RPDV)**RPWIP
0263             IF(ACCIP.LT.PYR(0)) GOTO 144
0264             OV=EXP(-B2RPW)/PARU(2)
0265             B=B2RPW**(1D0/POWIP)
0266           ELSE
0267   146       B2RPW=B2RPDV-2D0*LOG(PYR(0))
0268             ACCIP=(B2RPW/B2RPMX)**RPWIP*EXP(-0.5D0*(B2RPW-B2RPMX))
0269             IF(ACCIP.LT.PYR(0)) GOTO 146
0270             OV=EXP(-B2RPW)/PARU(2)
0271             B=B2RPW**(1D0/POWIP)
0272           ENDIF
0273           VINT(148)=OV/VNT147
0274           PACC=(1D0-EXP(-MIN(50D0,PIK*OV)))/(PIK*OV)
0275         ENDIF
0276         IF(PACC.LT.PYR(0)) GOTO 142
0277         VINT(139)=B/BAVG
0278
0279       ELSEIF(MMUL.EQ.3) THEN
0280 C...Low-pT or multiple interactions (first semihard interaction):
0281 C...choose xT2 according to dpT2/pT2**2*exp(-(sigma above pT2)/norm)
0282 C...or (MSTP(82)>=2) dpT2/(pT2+pT0**2)**2*exp(-....).
0283         ISUB=MINT(1)
0284         VINT(145)=VNT145
0285         VINT(146)=VNT146
0286         VINT(147)=VNT147
0287         IF(MSTP(82).LE.0) THEN
0288           XT2=0D0
0289         ELSEIF(MSTP(82).EQ.1) THEN
0290           XT2=XT2FAC*XT2/(XT2FAC-XT2*LOG(PYR(0)))
0291 C...Use with "Sudakov" for low b values when impact parameter dependence.
0292         ELSEIF(MSTP(82).EQ.2.OR.MINT(39).EQ.1) THEN
0293           IF(XT2.LT.1D0.AND.EXP(-XT2FAC*XT2/(VINT(149)*(XT2+
0294      &    VINT(149)))).GT.PYR(0)) XT2=1D0
0295           IF(XT2.GE.1D0) THEN
0296             XT2=(1D0+VINT(149))*XT2FAC/(XT2FAC-(1D0+VINT(149))*LOG(1D0-
0297      &      PYR(0)*(1D0-EXP(-XT2FAC/(VINT(149)*(1D0+VINT(149)))))))-
0298      &      VINT(149)
0299           ELSE
0300             XT2=-XT2FAC/LOG(EXP(-XT2FAC/(XT2+VINT(149)))+PYR(0)*
0301      &      (EXP(-XT2FAC/VINT(149))-EXP(-XT2FAC/(XT2+VINT(149)))))-
0302      &      VINT(149)
0303           ENDIF
0304           XT2=MAX(0.01D0*VINT(149),XT2)
0305 C...Use without "Sudakov" for high b values when impact parameter dep.
0306         ELSE
0307           XT2=(XC2+VINT(149))*(1D0+VINT(149))/(1D0+VINT(149)-
0308      &    PYR(0)*(1D0-XC2))-VINT(149)
0309           XT2=MAX(0.01D0*VINT(149),XT2)
0310         ENDIF
0311         VINT(25)=XT2
0312
0313 C...Low-pT: choose xT2, tau, y* and cos(theta-hat) fixed.
0314         IF(MSTP(82).LE.1.AND.XT2.LT.VINT(149)) THEN
0315           IF(MINT(82).EQ.1) NGEN(0,1)=NGEN(0,1)-MINT(143)
0316           IF(MINT(82).EQ.1) NGEN(ISUB,1)=NGEN(ISUB,1)-MINT(143)
0317           ISUB=95
0318           MINT(1)=ISUB
0319           VINT(21)=0.01D0*VINT(149)
0320           VINT(22)=0D0
0321           VINT(23)=0D0
0322           VINT(25)=0.01D0*VINT(149)
0323
0324         ELSE
0325 C...Multiple interactions (first semihard interaction).
0326 C...Choose tau and y*. Calculate cos(theta-hat).
0327           IF(PYR(0).LE.COEF(ISUB,1)) THEN
0328             TAUT=(2D0*(1D0+SQRT(1D0-XT2))/XT2-1D0)**PYR(0)
0329             TAU=XT2*(1D0+TAUT)**2/(4D0*TAUT)
0330           ELSE
0331             TAU=XT2*(1D0+TAN(PYR(0)*ATAN(SQRT(1D0/XT2-1D0)))**2)
0332           ENDIF
0333           VINT(21)=TAU
0334           CALL PYKLIM(2)
0335           RYST=PYR(0)
0336           MYST=1
0337           IF(RYST.GT.COEF(ISUB,8)) MYST=2
0338           IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
0339           CALL PYKMAP(2,MYST,PYR(0))
0340           VINT(23)=SQRT(MAX(0D0,1D0-XT2/TAU))*(-1)**INT(1.5D0+PYR(0))
0341         ENDIF
0342         VINT(71)=0.5D0*VINT(1)*SQRT(VINT(25))
0343
0344 C...Store results of cross-section calculation.
0345       ELSEIF(MMUL.EQ.4) THEN
0346         ISUB=MINT(1)
0347         VINT(145)=VNT145
0348         VINT(146)=VNT146
0349         VINT(147)=VNT147
0350         XTS=VINT(25)
0351         IF(ISET(ISUB).EQ.1) XTS=VINT(21)
0352         IF(ISET(ISUB).EQ.2)
0353      &  XTS=(4D0*VINT(48)+2D0*VINT(63)+2D0*VINT(64))/VINT(2)
0354         IF(ISET(ISUB).GE.3.AND.ISET(ISUB).LE.5) XTS=VINT(26)
0355         RBIN=MAX(0.000001D0,MIN(0.999999D0,XTS*(1D0+VINT(149))/
0356      &  (XTS+VINT(149))))
0357         IRBIN=INT(1D0+20D0*RBIN)
0358         IF(ISUB.EQ.96.AND.MSTP(171).EQ.0) THEN
0359           NMUL(IRBIN)=NMUL(IRBIN)+1
0360           SIGM(IRBIN)=SIGM(IRBIN)+VINT(153)
0361         ENDIF
0362
0363 C...Choose impact parameter if not already done.
0364       ELSEIF(MMUL.EQ.5) THEN
0365         ISUB=MINT(1)
0366         VINT(145)=VNT145
0367         VINT(146)=VNT146
0368         VINT(147)=VNT147
0369   150   IF(MINT(39).GT.0) THEN
0370         ELSEIF(MSTP(82).EQ.3) THEN
0371           EXPB2=PYR(0)
0372           B2=-LOG(PYR(0))
0373           VINT(148)=EXPB2/(PARU(2)*VNT147)
0374           VINT(139)=SQRT(B2)/BAVG
0375         ELSEIF(MSTP(82).EQ.4) THEN
0376           RTYPE=PYR(0)
0377           IF(RTYPE.LT.P83A) THEN
0378             B2=-LOG(PYR(0))
0379           ELSEIF(RTYPE.LT.P83A+P83B) THEN
0380             B2=-LOG(PYR(0))/CQ2R
0381           ELSE
0382             B2=-LOG(PYR(0))/CQ2I
0383           ENDIF
0384           VINT(148)=(P83A*EXP(-MIN(50D0,B2))+
0385      &    P83B*CQ2R*EXP(-MIN(50D0,B2*CQ2R))+
0386      &    P83C*CQ2I*EXP(-MIN(50D0,B2*CQ2I)))/(PARU(2)*VNT147)
0387           VINT(139)=SQRT(B2)/BAVG
0388         ELSEIF(PARP(83).GE.1.999D0) THEN
0389           POWIP=MAX(2D0,PARP(83))
0390           RPWIP=2D0/POWIP-1D0
0391           PROB1=POWIP/(2D0*EXP(-1D0)+POWIP)
0392   160     IF(PYR(0).LT.PROB1) THEN
0393             B2RPW=PYR(0)**(0.5D0*POWIP)
0394             ACCIP=EXP(-B2RPW)
0395           ELSE
0396             B2RPW=1D0-LOG(PYR(0))
0397             ACCIP=B2RPW**RPWIP
0398           ENDIF
0399           IF(ACCIP.LT.PYR(0)) GOTO 160
0400           VINT(148)=EXP(-B2RPW)/(PARU(2)*VNT147)
0401           VINT(139)=B2RPW**(1D0/POWIP)/BAVG
0402         ELSE
0403           POWIP=MAX(0.4D0,PARP(83))
0404           RPWIP=2D0/POWIP-1D0
0405           PROB1=RPWIP/(RPWIP+2D0**RPWIP*EXP(-RPWIP))
0406   170     IF(PYR(0).LT.PROB1) THEN
0407             B2RPW=2D0*RPWIP*PYR(0)
0408             ACCIP=(B2RPW/RPWIP)**RPWIP*EXP(RPWIP-B2RPW)
0409           ELSE
0410             B2RPW=2D0*(RPWIP-LOG(PYR(0)))
0411             ACCIP=(0.5D0*B2RPW/RPWIP)**RPWIP*EXP(RPWIP-0.5D0*B2RPW)
0412           ENDIF
0413           IF(ACCIP.LT .PYR(0)) GOTO 170
0414           VINT(148)=EXP(-B2RPW)/(PARU(2)*VNT147)
0415           VINT(139)=B2RPW**(1D0/POWIP)/BAVG
0416         ENDIF
0417
0418 C...Multiple interactions (variable impact parameter) : reject with
0419 C...probability exp(-overlap*cross-section above pT/normalization).
0420 C...Does not apply to low-b region, where "Sudakov" already included.
0421         VINT(150)=1D0
0422         IF(MINT(39).NE.1) THEN
0423           RNCOR=(IRBIN-20D0*RBIN)*NMUL(IRBIN)
0424           SIGCOR=(IRBIN-20D0*RBIN)*SIGM(IRBIN)
0425           DO 180 IBIN=IRBIN+1,20
0426             RNCOR=RNCOR+NMUL(IBIN)
0427             SIGCOR=SIGCOR+SIGM(IBIN)
0428   180     CONTINUE
0429           SIGABV=(SIGCOR/RNCOR)*VINT(149)*(1D0-XTS)/(XTS+VINT(149))
0430           IF(MSTP(171).EQ.1) SIGABV=SIGABV*VINT(2)/VINT(289)
0431           VINT(150)=EXP(-MIN(50D0,VNT146*VINT(148)*
0432      &    SIGABV/MAX(1D-10,SIGT(0,0,5))))
0433         ENDIF
0434         IF(MSTP(86).EQ.3.OR.(MSTP(86).EQ.2.AND.ISUB.NE.11.AND.
0435      &  ISUB.NE.12.AND.ISUB.NE.13.AND.ISUB.NE.28.AND.ISUB.NE.53
0436      &  .AND.ISUB.NE.68.AND.ISUB.NE.95.AND.ISUB.NE.96)) THEN
0437           IF(VINT(150).LT.PYR(0)) GOTO 150
0438           VINT(150)=1D0
0439         ENDIF
0440
0441 C...Generate additional multiple semihard interactions.
0442       ELSEIF(MMUL.EQ.6) THEN
0443         ISUBSV=MINT(1)
0444         VINT(145)=VNT145
0445         VINT(146)=VNT146
0446         VINT(147)=VNT147
0447         DO 190 J=11,80
0448           VINTSV(J)=VINT(J)
0449   190   CONTINUE
0450         ISUB=96
0451         MINT(1)=96
0452         VINT(151)=0D0
0453         VINT(152)=0D0
0454
0455 C...Reconstruct strings in hard scattering.
0456         NMAX=MINT(84)+4
0457         IF(ISET(ISUBSV).EQ.1) NMAX=MINT(84)+2
0458         IF(ISET(ISUBSV).EQ.11) NMAX=MINT(84)+2+MINT(3)
0459         NSTR=0
0460         DO 210 I=MINT(84)+1,NMAX
0461           KCS=KCHG(PYCOMP(K(I,2)),2)*ISIGN(1,K(I,2))
0462           IF(KCS.EQ.0) GOTO 210
0463           DO 200 J=1,4
0464             IF(KCS.EQ.1.AND.(J.EQ.2.OR.J.EQ.4)) GOTO 200
0465             IF(KCS.EQ.-1.AND.(J.EQ.1.OR.J.EQ.3)) GOTO 200
0466             IF(J.LE.2) THEN
0467               IST=MOD(K(I,J+3)/MSTU(5),MSTU(5))
0468             ELSE
0469               IST=MOD(K(I,J+1),MSTU(5))
0470             ENDIF
0471             IF(IST.LT.MINT(84).OR.IST.GT.I) GOTO 200
0472             IF(KCHG(PYCOMP(K(IST,2)),2).EQ.0) GOTO 200
0473             NSTR=NSTR+1
0474             IF(J.EQ.1.OR.J.EQ.4) THEN
0475               KSTR(NSTR,1)=I
0476               KSTR(NSTR,2)=IST
0477             ELSE
0478               KSTR(NSTR,1)=IST
0479               KSTR(NSTR,2)=I
0480             ENDIF
0481   200     CONTINUE
0482   210   CONTINUE
0483
0484 C...Set up starting values for iteration in xT2.
0485         XT2=4D0*VINT(62)/VINT(2)
0486         IF(MSTP(82).LE.1) THEN
0487           SIGRAT=XSEC(ISUB,1)/MAX(1D-10,VINT(315)*VINT(316)*SIGT(0,0,5))
0488           IF(MINT(141).NE.0.OR.MINT(142).NE.0) SIGRAT=SIGRAT*
0489      &    VINT(317)/(VINT(318)*VINT(320))
0490           XT2FAC=SIGRAT*VINT(149)/(1D0-VINT(149))
0491         ELSE
0492           XT2FAC=VNT146*VINT(148)*XSEC(ISUB,1)/
0493      &    MAX(1D-10,SIGT(0,0,5))*VINT(149)*(1D0+VINT(149))
0494         ENDIF
0495         VINT(63)=0D0
0496         VINT(64)=0D0
0497         VINT(143)=1D0-VINT(141)
0498         VINT(144)=1D0-VINT(142)
0499
0500 C...Iterate downwards in xT2.
0501   220   IF(MSTP(82).LE.1) THEN
0502           XT2=XT2FAC*XT2/(XT2FAC-XT2*LOG(PYR(0)))
0503           IF(XT2.LT.VINT(149)) GOTO 270
0504         ELSE
0505           IF(XT2.LE.0.01001D0*VINT(149)) GOTO 270
0506           XT2=XT2FAC*(XT2+VINT(149))/(XT2FAC-(XT2+VINT(149))*
0507      &    LOG(PYR(0)))-VINT(149)
0508           IF(XT2.LE.0D0) GOTO 270
0509           XT2=MAX(0.01D0*VINT(149),XT2)
0510         ENDIF
0511         VINT(25)=XT2
0512
0513 C...Choose tau and y*. Calculate cos(theta-hat).
0514         IF(PYR(0).LE.COEF(ISUB,1)) THEN
0515           TAUT=(2D0*(1D0+SQRT(1D0-XT2))/XT2-1D0)**PYR(0)
0516           TAU=XT2*(1D0+TAUT)**2/(4D0*TAUT)
0517         ELSE
0518           TAU=XT2*(1D0+TAN(PYR(0)*ATAN(SQRT(1D0/XT2-1D0)))**2)
0519         ENDIF
0520         VINT(21)=TAU
0521         CALL PYKLIM(2)
0522         RYST=PYR(0)
0523         MYST=1
0524         IF(RYST.GT.COEF(ISUB,8)) MYST=2
0525         IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
0526         CALL PYKMAP(2,MYST,PYR(0))
0527         VINT(23)=SQRT(MAX(0D0,1D0-XT2/TAU))*(-1)**INT(1.5D0+PYR(0))
0528
0529 C...Check that x not used up. Accept or reject kinematical variables.
0530         X1M=SQRT(TAU)*EXP(VINT(22))
0531         X2M=SQRT(TAU)*EXP(-VINT(22))
0532         IF(VINT(143)-X1M.LT.0.01D0.OR.VINT(144)-X2M.LT.0.01D0) GOTO 220
0533         VINT(71)=0.5D0*VINT(1)*SQRT(XT2)
0534         CALL PYSIGH(NCHN,SIGS)
0535         IF(MINT(141).NE.0.OR.MINT(142).NE.0) SIGS=SIGS*VINT(320)
0536         IF(SIGS.LT.XSEC(ISUB,1)*PYR(0)) GOTO 220
0537
0538 C...Reset K, P and V vectors. Select some variables.
0539         DO 240 I=N+1,N+2
0540           DO 230 J=1,5
0541             K(I,J)=0
0542             P(I,J)=0D0
0543             V(I,J)=0D0
0544   230     CONTINUE
0545   240   CONTINUE
0546         RFLAV=PYR(0)
0547         PT=0.5D0*VINT(1)*SQRT(XT2)
0548         PHI=PARU(2)*PYR(0)
0549         CTH=VINT(23)
0550
0551 C...Add first parton to event record.
0552         K(N+1,1)=3
0553         K(N+1,2)=21
0554         IF(RFLAV.GE.MAX(PARP(85),PARP(86))) K(N+1,2)=
0555      &  1+INT((2D0+PARJ(2))*PYR(0))
0556         P(N+1,1)=PT*COS(PHI)
0557         P(N+1,2)=PT*SIN(PHI)
0558         P(N+1,3)=0.25D0*VINT(1)*(VINT(41)*(1D0+CTH)-VINT(42)*(1D0-CTH))
0559         P(N+1,4)=0.25D0*VINT(1)*(VINT(41)*(1D0+CTH)+VINT(42)*(1D0-CTH))
0560         P(N+1,5)=0D0
0561
0562 C...Add second parton to event record.
0563         K(N+2,1)=3
0564         K(N+2,2)=21
0565         IF(K(N+1,2).NE.21) K(N+2,2)=-K(N+1,2)
0566         P(N+2,1)=-P(N+1,1)
0567         P(N+2,2)=-P(N+1,2)
0568         P(N+2,3)=0.25D0*VINT(1)*(VINT(41)*(1D0-CTH)-VINT(42)*(1D0+CTH))
0569         P(N+2,4)=0.25D0*VINT(1)*(VINT(41)*(1D0-CTH)+VINT(42)*(1D0+CTH))
0570         P(N+2,5)=0D0
0571
0572         IF(RFLAV.LT.PARP(85).AND.NSTR.GE.1) THEN
0573 C....Choose relevant string pieces to place gluons on.
0574           DO 260 I=N+1,N+2
0575             DMIN=1D8
0576             DO 250 ISTR=1,NSTR
0577               I1=KSTR(ISTR,1)
0578               I2=KSTR(ISTR,2)
0579               DIST=(P(I,4)*P(I1,4)-P(I,1)*P(I1,1)-P(I,2)*P(I1,2)-
0580      &        P(I,3)*P(I1,3))*(P(I,4)*P(I2,4)-P(I,1)*P(I2,1)-
0581      &        P(I,2)*P(I2,2)-P(I,3)*P(I2,3))/MAX(1D0,P(I1,4)*P(I2,4)-
0582      &        P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)-P(I1,3)*P(I2,3))
0583               IF(ISTR.EQ.1.OR.DIST.LT.DMIN) THEN
0584                 DMIN=DIST
0585                 IST1=I1
0586                 IST2=I2
0587                 ISTM=ISTR
0588               ENDIF
0589   250       CONTINUE
0590
0591 C....Colour flow adjustments, new string pieces.
0592             IF(K(IST1,4)/MSTU(5).EQ.IST2) K(IST1,4)=MSTU(5)*I+
0593      &      MOD(K(IST1,4),MSTU(5))
0594             IF(MOD(K(IST1,5),MSTU(5)).EQ.IST2) K(IST1,5)=
0595      &      MSTU(5)*(K(IST1,5)/MSTU(5))+I
0596             K(I,5)=MSTU(5)*IST1
0597             K(I,4)=MSTU(5)*IST2
0598             IF(K(IST2,5)/MSTU(5).EQ.IST1) K(IST2,5)=MSTU(5)*I+
0599      &      MOD(K(IST2,5),MSTU(5))
0600             IF(MOD(K(IST2,4),MSTU(5)).EQ.IST1) K(IST2,4)=
0601      &      MSTU(5)*(K(IST2,4)/MSTU(5))+I
0602             KSTR(ISTM,2)=I
0603             KSTR(NSTR+1,1)=I
0604             KSTR(NSTR+1,2)=IST2
0605             NSTR=NSTR+1
0606   260     CONTINUE
0607
0608 C...String drawing and colour flow for gluon loop.
0609         ELSEIF(K(N+1,2).EQ.21) THEN
0610           K(N+1,4)=MSTU(5)*(N+2)
0611           K(N+1,5)=MSTU(5)*(N+2)
0612           K(N+2,4)=MSTU(5)*(N+1)
0613           K(N+2,5)=MSTU(5)*(N+1)
0614           KSTR(NSTR+1,1)=N+1
0615           KSTR(NSTR+1,2)=N+2
0616           KSTR(NSTR+2,1)=N+2
0617           KSTR(NSTR+2,2)=N+1
0618           NSTR=NSTR+2
0619
0620 C...String drawing and colour flow for qqbar pair.
0621         ELSE
0622           K(N+1,4)=MSTU(5)*(N+2)
0623           K(N+2,5)=MSTU(5)*(N+1)
0624           KSTR(NSTR+1,1)=N+1
0625           KSTR(NSTR+1,2)=N+2
0626           NSTR=NSTR+1
0627         ENDIF
0628
0629 C...Global statistics.
0630         MINT(351)=MINT(351)+1
0631         VINT(351)=VINT(351)+PT
0632         IF (MINT(351).EQ.1) VINT(356)=PT
0633
0634 C...Update remaining energy; iterate.
0635         N=N+2
0636         IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
0637           CALL PYERRM(11,'(PYMULT:) no more memory left in PYJETS')
0638           MINT(51)=1
0639           RETURN
0640         ENDIF
0641         MINT(31)=MINT(31)+1
0642         VINT(151)=VINT(151)+VINT(41)
0643         VINT(152)=VINT(152)+VINT(42)
0644         VINT(143)=VINT(143)-VINT(41)
0645         VINT(144)=VINT(144)-VINT(42)
0646 C...Allow FSR for UE
0647         IF(MSTP(152).EQ.1) CALL PYSHOW(N-1,N,SQRT(PARP(71))*PT)
0648         IF(MINT(31).LT.240) GOTO 220
0649   270   CONTINUE
0650         MINT(1)=ISUBSV
0651         DO 280 J=11,80
0652           VINT(J)=VINTSV(J)
0653   280   CONTINUE
0654       ENDIF
0655
0656 C...Format statements for printout.
0657  5000 FORMAT(/1X,'****** PYMULT: initialization of multiple inter',
0658      &'actions for MSTP(82) =',I2,' ******')
0659  5100 FORMAT(8X,'pT0 =',F5.2,' GeV gives sigma(parton-parton) =',1P,
0660      &D9.2,' mb: rejected')
0661  5200 FORMAT(8X,'pT0 =',F5.2,' GeV gives sigma(parton-parton) =',1P,
0662      &D9.2,' mb: accepted')
0663
0664       RETURN
0665       END