PYTHIA/pythia/pymult.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE PYMULT(MMUL)
   5
   6 C...Initializes treatment of multiple interactions, selects kinematics
   7 C...of hardest interaction if low-pT physics included in run, and
   8 C...generates all non-hardest interactions.
   9       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
  10       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  11       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  12       COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
  13       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
  14       COMMON/PYINT1/MINT(400),VINT(400)
  15       COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
  16       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
  17       COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
  18       COMMON/PYINT7/SIGT(0:6,0:6,0:5)
  19       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
  20       SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT3/,/PYINT5/,
  21      &/PYINT7/
  22       DIMENSION NMUL(20),SIGM(20),KSTR(500,2),VINTSV(80)
  23       SAVE XT2,XT2FAC,XC2,XTS,IRBIN,RBIN,NMUL,SIGM
  24
  25 C...Initialization of multiple interaction treatment.
  26       IF(MMUL.EQ.1) THEN
  27         IF(MSTP(122).GE.1) WRITE(MSTU(11),5000) MSTP(82)
  28         ISUB=96
  29         MINT(1)=96
  30         VINT(63)=0.
  31         VINT(64)=0.
  32         VINT(143)=1.
  33         VINT(144)=1.
  34
  35 C...Loop over phase space points: xT2 choice in 20 bins.
  36   100   SIGSUM=0.
  37         DO 120 IXT2=1,20
  38         NMUL(IXT2)=MSTP(83)
  39         SIGM(IXT2)=0.
  40         DO 110 ITRY=1,MSTP(83)
  41         RSCA=0.05*((21-IXT2)-RLU(0))
  42         XT2=VINT(149)*(1.+VINT(149))/(VINT(149)+RSCA)-VINT(149)
  43         XT2=MAX(0.01*VINT(149),XT2)
  44         VINT(25)=XT2
  45
  46 C...Choose tau and y*. Calculate cos(theta-hat).
  47         IF(RLU(0).LE.COEF(ISUB,1)) THEN
  48           TAUT=(2.*(1.+SQRT(1.-XT2))/XT2-1.)**RLU(0)
  49           TAU=XT2*(1.+TAUT)**2/(4.*TAUT)
  50         ELSE
  51           TAU=XT2*(1.+TAN(RLU(0)*ATAN(SQRT(1./XT2-1.)))**2)
  52         ENDIF
  53         VINT(21)=TAU
  54         CALL PYKLIM(2)
  55         RYST=RLU(0)
  56         MYST=1
  57         IF(RYST.GT.COEF(ISUB,8)) MYST=2
  58         IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
  59         CALL PYKMAP(2,MYST,RLU(0))
  60         VINT(23)=SQRT(MAX(0.,1.-XT2/TAU))*(-1)**INT(1.5+RLU(0))
  61
  62 C...Calculate differential cross-section.
  63         VINT(71)=0.5*VINT(1)*SQRT(XT2)
  64         CALL PYSIGH(NCHN,SIGS)
  65         SIGM(IXT2)=SIGM(IXT2)+SIGS
  66   110   CONTINUE
  67         SIGSUM=SIGSUM+SIGM(IXT2)
  68   120   CONTINUE
  69         SIGSUM=SIGSUM/(20.*MSTP(83))
  70
  71 C...Reject result if sigma(parton-parton) is smaller than hadronic one.
  72         IF(SIGSUM.LT.1.1*SIGT(0,0,5)) THEN
  73           IF(MSTP(122).GE.1) WRITE(MSTU(11),5100) PARP(82),SIGSUM
  74           PARP(82)=0.9*PARP(82)
  75           VINT(149)=4.*PARP(82)**2/VINT(2)
  76           GOTO 100
  77         ENDIF
  78         IF(MSTP(122).GE.1) WRITE(MSTU(11),5200) PARP(82), SIGSUM
  79
  80 C...Start iteration to find k factor.
  81         YKE=SIGSUM/SIGT(0,0,5)
  82         SO=0.5
  83         XI=0.
  84         YI=0.
  85         XF=0.
  86         YF=0.
  87         XK=0.5
  88         IIT=0
  89   130   IF(IIT.EQ.0) THEN
  90           XK=2.*XK
  91         ELSEIF(IIT.EQ.1) THEN
  92           XK=0.5*XK
  93         ELSE
  94           XK=XI+(YKE-YI)*(XF-XI)/(YF-YI)
  95         ENDIF
  96
  97 C...Evaluate overlap integrals.
  98         IF(MSTP(82).EQ.2) THEN
  99           SP=0.5*PARU(1)*(1.-EXP(-XK))
 100           SOP=SP/PARU(1)
 101         ELSE
 102           IF(MSTP(82).EQ.3) DELTAB=0.02
 103           IF(MSTP(82).EQ.4) DELTAB=MIN(0.01,0.05*PARP(84))
 104           SP=0.
 105           SOP=0.
 106           B=-0.5*DELTAB
 107   140     B=B+DELTAB
 108           IF(MSTP(82).EQ.3) THEN
 109             OV=EXP(-B**2)/PARU(2)
 110           ELSE
 111             CQ2=PARP(84)**2
 112             OV=((1.-PARP(83))**2*EXP(-MIN(50.,B**2))+2.*PARP(83)*
 113      &      (1.-PARP(83))*2./(1.+CQ2)*EXP(-MIN(50.,B**2*2./(1.+CQ2)))+
 114      &      PARP(83)**2/CQ2*EXP(-MIN(50.,B**2/CQ2)))/PARU(2)
 115           ENDIF
 116           PACC=1.-EXP(-MIN(50.,PARU(1)*XK*OV))
 117           SP=SP+PARU(2)*B*DELTAB*PACC
 118           SOP=SOP+PARU(2)*B*DELTAB*OV*PACC
 119           IF(B.LT.1..OR.B*PACC.GT.1E-6) GOTO 140
 120         ENDIF
 121         YK=PARU(1)*XK*SO/SP
 122
 123 C...Continue iteration until convergence.
 124         IF(YK.LT.YKE) THEN
 125           XI=XK
 126           YI=YK
 127           IF(IIT.EQ.1) IIT=2
 128         ELSE
 129           XF=XK
 130           YF=YK
 131           IF(IIT.EQ.0) IIT=1
 132         ENDIF
 133         IF(ABS(YK-YKE).GE.1E-5*YKE) GOTO 130
 134
 135 C...Store some results for subsequent use.
 136         VINT(145)=SIGSUM
 137         VINT(146)=SOP/SO
 138         VINT(147)=SOP/SP
 139
 140 C...Initialize iteration in xT2 for hardest interaction.
 141       ELSEIF(MMUL.EQ.2) THEN
 142         IF(MSTP(82).LE.0) THEN
 143         ELSEIF(MSTP(82).EQ.1) THEN
 144           XT2=1.
 145           XT2FAC=XSEC(96,1)/SIGT(0,0,5)*VINT(149)/(1.-VINT(149))
 146         ELSEIF(MSTP(82).EQ.2) THEN
 147           XT2=1.
 148           XT2FAC=VINT(146)*XSEC(96,1)/SIGT(0,0,5)*VINT(149)*
 149      &    (1.+VINT(149))
 150         ELSE
 151           XC2=4.*CKIN(3)**2/VINT(2)
 152           IF(CKIN(3).LE.CKIN(5).OR.MINT(82).GE.2) XC2=0.
 153         ENDIF
 154
 155       ELSEIF(MMUL.EQ.3) THEN
 156 C...Low-pT or multiple interactions (first semihard interaction):
 157 C...choose xT2 according to dpT2/pT2**2*exp(-(sigma above pT2)/norm)
 158 C...or (MSTP(82)>=2) dpT2/(pT2+pT0**2)**2*exp(-....).
 159         ISUB=MINT(1)
 160         IF(MSTP(82).LE.0) THEN
 161           XT2=0.
 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.
 167           IF(XT2.GE.1.) THEN
 168             XT2=(1.+VINT(149))*XT2FAC/(XT2FAC-(1.+VINT(149))*LOG(1.-
 169      &      RLU(0)*(1.-EXP(-XT2FAC/(VINT(149)*(1.+VINT(149)))))))-
 170      &      VINT(149)
 171           ELSE
 172             XT2=-XT2FAC/LOG(EXP(-XT2FAC/(XT2+VINT(149)))+RLU(0)*
 173      &      (EXP(-XT2FAC/VINT(149))-EXP(-XT2FAC/(XT2+VINT(149)))))-
 174      &      VINT(149)
 175           ENDIF
 176           XT2=MAX(0.01*VINT(149),XT2)
 177         ELSE
 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)
 181         ENDIF
 182         VINT(25)=XT2
 183
 184 C...Low-pT: choose xT2, tau, y* and cos(theta-hat) fixed.
 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
 188           ISUB=95
 189           MINT(1)=ISUB
 190           VINT(21)=0.01*VINT(149)
 191           VINT(22)=0.
 192           VINT(23)=0.
 193           VINT(25)=0.01*VINT(149)
 194
 195         ELSE
 196 C...Multiple interactions (first semihard interaction).
 197 C...Choose tau and y*. Calculate cos(theta-hat).
 198           IF(RLU(0).LE.COEF(ISUB,1)) THEN
 199             TAUT=(2.*(1.+SQRT(1.-XT2))/XT2-1.)**RLU(0)
 200             TAU=XT2*(1.+TAUT)**2/(4.*TAUT)
 201           ELSE
 202             TAU=XT2*(1.+TAN(RLU(0)*ATAN(SQRT(1./XT2-1.)))**2)
 203           ENDIF
 204           VINT(21)=TAU
 205           CALL PYKLIM(2)
 206           RYST=RLU(0)
 207           MYST=1
 208           IF(RYST.GT.COEF(ISUB,8)) MYST=2
 209           IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
 210           CALL PYKMAP(2,MYST,RLU(0))
 211           VINT(23)=SQRT(MAX(0.,1.-XT2/TAU))*(-1)**INT(1.5+RLU(0))
 212         ENDIF
 213         VINT(71)=0.5*VINT(1)*SQRT(VINT(25))
 214
 215 C...Store results of cross-section calculation.
 216       ELSEIF(MMUL.EQ.4) THEN
 217         ISUB=MINT(1)
 218         XTS=VINT(25)
 219         IF(ISET(ISUB).EQ.1) XTS=VINT(21)
 220         IF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.6)
 221      &  XTS=(4.*VINT(48)+2.*VINT(63)+2.*VINT(64))/VINT(2)
 222         IF(ISET(ISUB).GE.3.AND.ISET(ISUB).LE.5) XTS=VINT(26)
 223         RBIN=MAX(0.000001,MIN(0.999999,XTS*(1.+VINT(149))/
 224      &  (XTS+VINT(149))))
 225         IRBIN=INT(1.+20.*RBIN)
 226         IF(ISUB.EQ.96.AND.MSTP(171).EQ.0) THEN
 227           NMUL(IRBIN)=NMUL(IRBIN)+1
 228           SIGM(IRBIN)=SIGM(IRBIN)+VINT(153)
 229         ENDIF
 230
 231 C...Choose impact parameter.
 232       ELSEIF(MMUL.EQ.5) THEN
 233         IF(MSTP(82).EQ.3) THEN
 234           VINT(148)=RLU(0)/(PARU(2)*VINT(147))
 235         ELSE
 236           RTYPE=RLU(0)
 237           CQ2=PARP(84)**2
 238           IF(RTYPE.LT.(1.-PARP(83))**2) THEN
 239             B2=-LOG(RLU(0))
 240           ELSEIF(RTYPE.LT.1.-PARP(83)**2) THEN
 241             B2=-0.5*(1.+CQ2)*LOG(RLU(0))
 242           ELSE
 243             B2=-CQ2*LOG(RLU(0))
 244           ENDIF
 245           VINT(148)=((1.-PARP(83))**2*EXP(-MIN(50.,B2))+2.*PARP(83)*
 246      &    (1.-PARP(83))*2./(1.+CQ2)*EXP(-MIN(50.,B2*2./(1.+CQ2)))+
 247      &    PARP(83)**2/CQ2*EXP(-MIN(50.,B2/CQ2)))/(PARU(2)*VINT(147))
 248         ENDIF
 249
 250 C...Multiple interactions (variable impact parameter) : reject with
 251 C...probability exp(-overlap*cross-section above pT/normalization).
 252         RNCOR=(IRBIN-20.*RBIN)*NMUL(IRBIN)
 253         SIGCOR=(IRBIN-20.*RBIN)*SIGM(IRBIN)
 254         DO 150 IBIN=IRBIN+1,20
 255         RNCOR=RNCOR+NMUL(IBIN)
 256         SIGCOR=SIGCOR+SIGM(IBIN)
 257   150   CONTINUE
 258         SIGABV=(SIGCOR/RNCOR)*VINT(149)*(1.-XTS)/(XTS+VINT(149))
 259         IF(MSTP(171).EQ.1) SIGABV=SIGABV*VINT(2)/VINT(289)
 260         VINT(150)=EXP(-MIN(50.,VINT(146)*VINT(148)*
 261      &  SIGABV/SIGT(0,0,5)))
 262
 263 C...Generate additional multiple semihard interactions.
 264       ELSEIF(MMUL.EQ.6) THEN
 265         ISUBSV=MINT(1)
 266         DO 160 J=11,80
 267         VINTSV(J)=VINT(J)
 268   160   CONTINUE
 269         ISUB=96
 270         MINT(1)=96
 271
 272 C...Reconstruct strings in hard scattering.
 273         NMAX=MINT(84)+4
 274         IF(ISET(ISUBSV).EQ.1) NMAX=MINT(84)+2
 275         IF(ISET(ISUBSV).EQ.11) NMAX=MINT(84)+2+MINT(3)
 276         NSTR=0
 277         DO 180 I=MINT(84)+1,NMAX
 278         KCS=KCHG(LUCOMP(K(I,2)),2)*ISIGN(1,K(I,2))
 279         IF(KCS.EQ.0) GOTO 180
 280
 281         DO 170 J=1,4
 282         IF(KCS.EQ.1.AND.(J.EQ.2.OR.J.EQ.4)) GOTO 170
 283         IF(KCS.EQ.-1.AND.(J.EQ.1.OR.J.EQ.3)) GOTO 170
 284         IF(J.LE.2) THEN
 285           IST=MOD(K(I,J+3)/MSTU(5),MSTU(5))
 286         ELSE
 287           IST=MOD(K(I,J+1),MSTU(5))
 288         ENDIF
 289         IF(IST.LT.MINT(84).OR.IST.GT.I) GOTO 170
 290         IF(KCHG(LUCOMP(K(IST,2)),2).EQ.0) GOTO 170
 291         NSTR=NSTR+1
 292         IF(J.EQ.1.OR.J.EQ.4) THEN
 293           KSTR(NSTR,1)=I
 294           KSTR(NSTR,2)=IST
 295         ELSE
 296           KSTR(NSTR,1)=IST
 297           KSTR(NSTR,2)=I
 298         ENDIF
 299   170   CONTINUE
 300   180   CONTINUE
 301
 302 C...Set up starting values for iteration in xT2.
 303         XT2=VINT(25)
 304         IF(ISET(ISUBSV).EQ.1) XT2=VINT(21)
 305         IF(ISET(ISUBSV).EQ.2.OR.ISET(ISUBSV).EQ.6)
 306      &  XT2=(4.*VINT(48)+2.*VINT(63)+2.*VINT(64))/VINT(2)
 307         IF(ISET(ISUBSV).GE.3.AND.ISET(ISUBSV).LE.5) XT2=VINT(26)
 308         IF(MSTP(82).LE.1) THEN
 309           XT2FAC=XSEC(ISUB,1)*VINT(149)/((1.-VINT(149))*SIGT(0,0,5))
 310         ELSE
 311           XT2FAC=VINT(146)*VINT(148)*XSEC(ISUB,1)/SIGT(0,0,5)*
 312      &    VINT(149)*(1.+VINT(149))
 313         ENDIF
 314         VINT(63)=0.
 315         VINT(64)=0.
 316         VINT(143)=1.-VINT(141)
 317         VINT(144)=1.-VINT(142)
 318
 319 C...Iterate downwards in xT2.
 320   190   IF(MSTP(82).LE.1) THEN
 321           XT2=XT2FAC*XT2/(XT2FAC-XT2*LOG(RLU(0)))
 322           IF(XT2.LT.VINT(149)) GOTO 240
 323         ELSE
 324           IF(XT2.LE.0.01*VINT(149)) GOTO 240
 325           XT2=XT2FAC*(XT2+VINT(149))/(XT2FAC-(XT2+VINT(149))*
 326      &    LOG(RLU(0)))-VINT(149)
 327           IF(XT2.LE.0.) GOTO 240
 328           XT2=MAX(0.01*VINT(149),XT2)
 329         ENDIF
 330         VINT(25)=XT2
 331
 332 C...Choose tau and y*. Calculate cos(theta-hat).
 333         IF(RLU(0).LE.COEF(ISUB,1)) THEN
 334           TAUT=(2.*(1.+SQRT(1.-XT2))/XT2-1.)**RLU(0)
 335           TAU=XT2*(1.+TAUT)**2/(4.*TAUT)
 336         ELSE
 337           TAU=XT2*(1.+TAN(RLU(0)*ATAN(SQRT(1./XT2-1.)))**2)
 338         ENDIF
 339         VINT(21)=TAU
 340         CALL PYKLIM(2)
 341         RYST=RLU(0)
 342         MYST=1
 343         IF(RYST.GT.COEF(ISUB,8)) MYST=2
 344         IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
 345         CALL PYKMAP(2,MYST,RLU(0))
 346         VINT(23)=SQRT(MAX(0.,1.-XT2/TAU))*(-1)**INT(1.5+RLU(0))
 347
 348 C...Check that x not used up. Accept or reject kinematical variables.
 349         X1M=SQRT(TAU)*EXP(VINT(22))
 350         X2M=SQRT(TAU)*EXP(-VINT(22))
 351         IF(VINT(143)-X1M.LT.0.01.OR.VINT(144)-X2M.LT.0.01) GOTO 190
 352         VINT(71)=0.5*VINT(1)*SQRT(XT2)
 353         CALL PYSIGH(NCHN,SIGS)
 354         IF(SIGS.LT.XSEC(ISUB,1)*RLU(0)) GOTO 190
 355
 356 C...Reset K, P and V vectors. Select some variables.
 357         DO 210 I=N+1,N+2
 358         DO 200 J=1,5
 359         K(I,J)=0
 360         P(I,J)=0.
 361         V(I,J)=0.
 362   200   CONTINUE
 363   210   CONTINUE
 364         RFLAV=RLU(0)
 365         PT=0.5*VINT(1)*SQRT(XT2)
 366         PHI=PARU(2)*RLU(0)
 367         CTH=VINT(23)
 368
 369 C...Add first parton to event record.
 370         K(N+1,1)=3
 371         K(N+1,2)=21
 372         IF(RFLAV.GE.MAX(PARP(85),PARP(86))) K(N+1,2)=
 373      &  1+INT((2.+PARJ(2))*RLU(0))
 374         P(N+1,1)=PT*COS(PHI)
 375         P(N+1,2)=PT*SIN(PHI)
 376         P(N+1,3)=0.25*VINT(1)*(VINT(41)*(1.+CTH)-VINT(42)*(1.-CTH))
 377         P(N+1,4)=0.25*VINT(1)*(VINT(41)*(1.+CTH)+VINT(42)*(1.-CTH))
 378         P(N+1,5)=0.
 379
 380 C...Add second parton to event record.
 381         K(N+2,1)=3
 382         K(N+2,2)=21
 383         IF(K(N+1,2).NE.21) K(N+2,2)=-K(N+1,2)
 384         P(N+2,1)=-P(N+1,1)
 385         P(N+2,2)=-P(N+1,2)
 386         P(N+2,3)=0.25*VINT(1)*(VINT(41)*(1.-CTH)-VINT(42)*(1.+CTH))
 387         P(N+2,4)=0.25*VINT(1)*(VINT(41)*(1.-CTH)+VINT(42)*(1.+CTH))
 388         P(N+2,5)=0.
 389
 390         IF(RFLAV.LT.PARP(85).AND.NSTR.GE.1) THEN
 391 C....Choose relevant string pieces to place gluons on.
 392           DO 230 I=N+1,N+2
 393           DMIN=1E8
 394           DO 220 ISTR=1,NSTR
 395           I1=KSTR(ISTR,1)
 396           I2=KSTR(ISTR,2)
 397           DIST=(P(I,4)*P(I1,4)-P(I,1)*P(I1,1)-P(I,2)*P(I1,2)-
 398      &    P(I,3)*P(I1,3))*(P(I,4)*P(I2,4)-P(I,1)*P(I2,1)-
 399      &    P(I,2)*P(I2,2)-P(I,3)*P(I2,3))/MAX(1.,P(I1,4)*P(I2,4)-
 400      &    P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)-P(I1,3)*P(I2,3))
 401           IF(ISTR.EQ.1.OR.DIST.LT.DMIN) THEN
 402             DMIN=DIST
 403             IST1=I1
 404             IST2=I2
 405             ISTM=ISTR
 406           ENDIF
 407   220     CONTINUE
 408
 409 C....Colour flow adjustments, new string pieces.
 410           IF(K(IST1,4)/MSTU(5).EQ.IST2) K(IST1,4)=MSTU(5)*I+
 411      &    MOD(K(IST1,4),MSTU(5))
 412           IF(MOD(K(IST1,5),MSTU(5)).EQ.IST2) K(IST1,5)=
 413      &    MSTU(5)*(K(IST1,5)/MSTU(5))+I
 414           K(I,5)=MSTU(5)*IST1
 415           K(I,4)=MSTU(5)*IST2
 416           IF(K(IST2,5)/MSTU(5).EQ.IST1) K(IST2,5)=MSTU(5)*I+
 417      &    MOD(K(IST2,5),MSTU(5))
 418           IF(MOD(K(IST2,4),MSTU(5)).EQ.IST1) K(IST2,4)=
 419      &    MSTU(5)*(K(IST2,4)/MSTU(5))+I
 420           KSTR(ISTM,2)=I
 421           KSTR(NSTR+1,1)=I
 422           KSTR(NSTR+1,2)=IST2
 423           NSTR=NSTR+1
 424   230     CONTINUE
 425
 426 C...String drawing and colour flow for gluon loop.
 427         ELSEIF(K(N+1,2).EQ.21) THEN
 428           K(N+1,4)=MSTU(5)*(N+2)
 429           K(N+1,5)=MSTU(5)*(N+2)
 430           K(N+2,4)=MSTU(5)*(N+1)
 431           K(N+2,5)=MSTU(5)*(N+1)
 432           KSTR(NSTR+1,1)=N+1
 433           KSTR(NSTR+1,2)=N+2
 434           KSTR(NSTR+2,1)=N+2
 435           KSTR(NSTR+2,2)=N+1
 436           NSTR=NSTR+2
 437
 438 C...String drawing and colour flow for qq~ pair.
 439         ELSE
 440           K(N+1,4)=MSTU(5)*(N+2)
 441           K(N+2,5)=MSTU(5)*(N+1)
 442           KSTR(NSTR+1,1)=N+1
 443           KSTR(NSTR+1,2)=N+2
 444           NSTR=NSTR+1
 445         ENDIF
 446
 447 C...Update remaining energy; iterate.
 448         N=N+2
 449         IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
 450           CALL LUERRM(11,'(PYMULT:) no more memory left in LUJETS')
 451           IF(MSTU(21).GE.1) RETURN
 452         ENDIF
 453         MINT(31)=MINT(31)+1
 454         VINT(151)=VINT(151)+VINT(41)
 455         VINT(152)=VINT(152)+VINT(42)
 456         VINT(143)=VINT(143)-VINT(41)
 457         VINT(144)=VINT(144)-VINT(42)
 458         IF(MINT(31).LT.240) GOTO 190
 459   240   CONTINUE
 460         MINT(1)=ISUBSV
 461         DO 250 J=11,80
 462         VINT(J)=VINTSV(J)
 463   250   CONTINUE
 464       ENDIF
 465
 466 C...Format statements for printout.
 467  5000 FORMAT(/1X,'****** PYMULT: initialization of multiple inter',
 468      &'actions for MSTP(82) =',I2,' ******')
 469  5100 FORMAT(8X,'pT0 =',F5.2,' GeV gives sigma(parton-parton) =',1P,
 470      &E9.2,' mb: rejected')
 471  5200 FORMAT(8X,'pT0 =',F5.2,' GeV gives sigma(parton-parton) =',1P,
 472      &E9.2,' mb: accepted')
 473
 474       RETURN
 475       END