pythiaeRHIC/pythia/pyptis.f

0001 C*********************************************************************
0002
0003 C...PYPTIS
0004 C...Generates pT-ordered spacelike initial-state parton showers and
0005 C...trial joinings.
0006 C...MODE=-1: Initialize ISR from scratch, starting from the hardest
0007 C...         interaction initiators at PT2NOW.
0008 C...MODE= 0: Generate a trial branching on interaction MINT(36), side
0009 C...         MINT(30). Start evolution at PT2NOW, solve Sudakov for PT2.
0010 C...         Store in /PYISMX/ if PT2 is largest so far. Abort if PT2
0011 C...         is below PT2CUT.
0012 C...         (Also generate test joinings if MSTP(96)=1.)
0013 C...MODE= 1: Accept stored shower branching. Update event record etc.
0014 C...PT2NOW : Starting (max) PT2 scale for evolution.
0015 C...PT2CUT : Lower limit for evolution.
0016 C...PT2    : Result of evolution. Generated PT2 for trial emission.
0017 C...IFAIL  : Status return code. IFAIL=0 when all is well.
0018
0019       SUBROUTINE PYPTIS(MODE,PT2NOW,PT2CUT,PT2,IFAIL)
0020
0021 C...Double precision and integer declarations.
0022       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0023       IMPLICIT INTEGER(I-N)
0024       INTEGER PYK,PYCHGE,PYCOMP
0025 C...Parameter statement for maximum size of showers.
0026       PARAMETER (MAXNUR=1000)
0027 C...Commonblocks.
0028       COMMON/PYPART/NPART,NPARTD,IPART(MAXNUR),PTPART(MAXNUR)
0029       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0030       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0031       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0032       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0033       COMMON/PYINT1/MINT(400),VINT(400)
0034       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0035       COMMON/PYINTM/KFIVAL(2,3),NMI(2),IMI(2,800,2),NVC(2,-6:6),
0036      &     XASSOC(2,-6:6,240),XPSVC(-6:6,-1:240),PVCTOT(2,-1:1),
0037      &     XMI(2,240),PT2MI(240),IMISEP(0:240)
0038       COMMON/PYISMX/MIMX,JSMX,KFLAMX,KFLCMX,KFBEAM(2),NISGEN(2,240),
0039      &     PT2MX,PT2AMX,ZMX,RM2CMX,Q2BMX,PHIMX
0040       COMMON/PYCTAG/NCT,MCT(4000,2)
0041       COMMON/PYISJN/MJN1MX,MJN2MX,MJOIND(2,240)
0042       SAVE /PYPART/,/PYJETS/,/PYDAT1/,/PYDAT2/,/PYPARS/,/PYINT1/,
0043      &     /PYINT2/,/PYINTM/,/PYISMX/,/PYCTAG/,/PYISJN/
0044 C...Local variables
0045       DIMENSION ZSAV(2,240),PT2SAV(2,240),
0046      &     XFB(-25:25),XFA(-25:25),XFN(-25:25),XFJ(-25:25),
0047      &     WTAP(-25:25),WTPDF(-25:25),SHTNOW(240),
0048      &     WTAPJ(240),WTPDFJ(240),X1(240),Y(240)
0049       SAVE ZSAV,PT2SAV,XFB,XFA,XFN,WTAP,WTPDF,XMXC,SHTNOW,
0050      &     RMB2,RMC2,ALAM3,ALAM4,ALAM5,TMIN,PTEMAX,WTEMAX,AEM2PI
0051 C...For check on excessive weights.
0052       CHARACTER CHWT*12
0053
0054       DATA PTEMAX /0D0/
0055       DATA WTEMAX /0D0/
0056
0057       IFAIL=-1
0058
0059 C----------------------------------------------------------------------
0060 C...MODE=-1: Initialize initial state showers from scratch, i.e.
0061 C...starting from the hardest interaction initiators.
0062       IF (MODE.EQ.-1) THEN
0063 C...Set hard scattering SHAT.
0064         SHTNOW(1)=VINT(44)
0065 C...Mass thresholds and Lambda for QCD evolution.
0066         AEM2PI=PARU(101)/PARU(2)
0067         RMB=PMAS(5,1)
0068         RMC=PMAS(4,1)
0069         ALAM4=PARP(61)
0070         IF(MSTU(112).LT.4) ALAM4=PARP(61)*(PARP(61)/RMC)**(2D0/25D0)
0071         IF(MSTU(112).GT.4) ALAM4=PARP(61)*(RMB/PARP(61))**(2D0/25D0)
0072         ALAM5=ALAM4*(ALAM4/RMB)**(2D0/23D0)
0073         ALAM3=ALAM4*(RMC/ALAM4)**(2D0/27D0)
0074         RMB2=RMB**2
0075         RMC2=RMC**2
0076 C...Massive quark forced creation threshold (in M**2).
0077         TMIN=1.01D0
0078 C...Set upper limit for X (ensures some X left for beam remnant).
0079         XMXC=1D0-2D0*PARP(111)/VINT(1)
0080
0081         IF (MSTP(61).GE.1) THEN
0082 C...Initial values: flavours, momenta, virtualities.
0083           DO 100 JS=1,2
0084             NISGEN(JS,1)=0
0085
0086 C...Special kinematics check for c/b quarks (that g -> c cbar or
0087 C...b bbar kinematically possible).
0088             KFLB=K(IMI(JS,1,1),2)
0089             KFLCB=IABS(KFLB)
0090             IF(KFBEAM(JS).NE.22.AND.(KFLCB.EQ.4.OR.KFLCB.EQ.5)) THEN
0091 C...Check PT2MAX > mQ^2
0092               IF (VINT(56).LT.1.05D0*PMAS(PYCOMP(KFLCB),1)**2) THEN
0093                 CALL PYERRM(9,'(PYPTIS:) PT2MAX < 1.05 * MQ**2. '//
0094      &               'No Q creation possible.')
0095                 MINT(51)=1
0096                 RETURN
0097               ELSE
0098 C...Check for physical z values (m == MQ / sqrt(s))
0099 C...For creation diagram, x < z < (1-m)/(1+m(1-m))
0100                 FMQ=PMAS(KFLCB,1)/SQRT(SHTNOW(1))
0101                 ZMXCR=(1D0-FMQ)/(1D0+FMQ*(1D0-FMQ))
0102                 IF (XMI(JS,1).GT.0.9D0*ZMXCR) THEN
0103                   CALL PYERRM(9,'(PYPTIS:) No physical z value for '//
0104      &                 'Q creation.')
0105                   MINT(51)=1
0106                   RETURN
0107                 ENDIF
0108               ENDIF
0109             ENDIF
0110   100     CONTINUE
0111         ENDIF
0112
0113         MINT(354)=0
0114 C...Zero joining array
0115         DO 110 MJ=1,240
0116           MJOIND(1,MJ)=0
0117           MJOIND(2,MJ)=0
0118   110   CONTINUE
0119
0120 C----------------------------------------------------------------------
0121 C...MODE= 0: Generate a trial branching on interaction MINT(36) side
0122 C...MINT(30). Store if emission PT2 scale is largest so far.
0123 C...Also generate test joinings if MSTP(96)=1.
0124       ELSEIF(MODE.EQ.0) THEN
0125         IFAIL=-1
0126         MECOR=0
0127         ISUB=MINT(1)
0128         JS=MINT(30)
0129 C...No shower for structureless beam
0130         IF (MINT(44+JS).EQ.1) RETURN
0131         MI=MINT(36)
0132         SHAT=VINT(44)
0133 C...Absolute shower max scale = VINT(56)
0134         PT2=MIN(PT2NOW,VINT(56))
0135         IF (NISGEN(1,MI).EQ.0.AND.NISGEN(2,MI).EQ.0) SHTNOW(MI)=SHAT
0136 C...Define for which processes ME corrections have been implemented.
0137         IF(MSTP(68).EQ.1.OR.MSTP(68).EQ.3) THEN
0138           IF(ISUB.EQ.1.OR.ISUB.EQ.2.OR.ISUB.EQ.141.OR.ISUB.EQ
0139      &         .142.OR.ISUB.EQ.144) MECOR=1
0140           IF(ISUB.EQ.102.OR.ISUB.EQ.152.OR.ISUB.EQ.157) MECOR=2
0141           IF(ISUB.EQ.3.OR.ISUB.EQ.151.OR.ISUB.EQ.156) MECOR=3
0142 C...Calculate preweighting factor for ME-corrected processes.
0143           IF(MECOR.GE.1) CALL PYMEMX(MECOR,WTFF,WTGF,WTFG,WTGG)
0144         ENDIF
0145 C...Basic info on daughter for which to find mother.
0146         KFLB=K(IMI(JS,MI,1),2)
0147         KFLBA=IABS(KFLB)
0148 C...KSVCB: -1 for sea or first companion, 0 for valence or gluon, >1 for
0149 C...second companion.
0150         KSVCB=MAX(-1,IMI(JS,MI,2))
0151 C...Treat "first" companion of a pair like an ordinary sea quark
0152 C...(except that creation diagram is not allowed)
0153         IF(IMI(JS,MI,2).GT.IMISEP(MI)) KSVCB=-1
0154 C...X (rescaled to [0,1])
0155         XB=XMI(JS,MI)/VINT(142+JS)
0156 C...Massive quarks (use physical masses.)
0157         RMQ2=0D0
0158         MQMASS=0
0159         IF (KFLBA.EQ.4.OR.KFLBA.EQ.5) THEN
0160           RMQ2=RMC2
0161           IF (KFLBA.EQ.5) RMQ2=RMB2
0162 C...Special threshold treatment for non-photon beams
0163           IF (KFBEAM(JS).NE.22) MQMASS=KFLBA
0164         ENDIF
0165
0166 C...Flags for parton distribution calls.
0167         MINT(105)=MINT(102+JS)
0168         MINT(109)=MINT(106+JS)
0169         VINT(120)=VINT(2+JS)
0170
0171 C...Calculate initial parton distribution weights.
0172         IF(XB.GE.XMXC) THEN
0173           RETURN
0174         ELSEIF(MQMASS.EQ.0) THEN
0175           CALL PYPDFU(KFBEAM(JS),XB,PT2,XFB)
0176         ELSE
0177 C...Initialize massive quark PT2 dependent pdf underestimate.
0178           PT20=PT2
0179           CALL PYPDFU(KFBEAM(JS),XB,PT20,XFB)
0180 C.!.Tentative treatment of massive valence quarks.
0181           XQ0=MAX(1D-10,XPSVC(KFLB,KSVCB))
0182           XG0=XFB(21)
0183           TPM0=LOG(PT20/RMQ2)
0184           WPDF0=TPM0*XG0/XQ0
0185         ENDIF
0186         IF (KFLBA.LE.6) THEN
0187 C...For quarks, only include respective sea, val, or cmp part.
0188           IF (KSVCB.LE.0) THEN
0189             XFB(KFLB)=XPSVC(KFLB,KSVCB)
0190           ELSE
0191 C...Find companion's companion
0192             MISEA=0
0193   120       MISEA=MISEA+1
0194             IF (IMI(JS,MISEA,2).NE.IMI(JS,MI,1)) GOTO 120
0195             XS=XMI(JS,MISEA)
0196             XREM=VINT(142+JS)
0197             YS=XS/(XREM+XS)
0198 C...Momentum fraction of the companion quark.
0199 C...Rescale from XB = x/XREM to YB = x/(1-Sum_rest) -> factor (1-YS).
0200             YB=XB*(1D0-YS)
0201             XFB(KFLB)=PYFCMP(YB/VINT(140),YS/VINT(140),MSTP(87))
0202           ENDIF
0203         ENDIF
0204
0205 C...Determine overestimated z range: switch at c and b masses.
0206   130   IF (PT2.GT.TMIN*RMB2) THEN
0207           IZRG=3
0208           PT2MNE=MAX(TMIN*RMB2,PT2CUT)
0209           B0=23D0/6D0
0210           ALAM2=ALAM5**2
0211         ELSEIF(PT2.GT.TMIN*RMC2) THEN
0212           IZRG=2
0213           PT2MNE=MAX(TMIN*RMC2,PT2CUT)
0214           B0=25D0/6D0
0215           ALAM2=ALAM4**2
0216         ELSE
0217           IZRG=1
0218           PT2MNE=PT2CUT
0219           B0=27D0/6D0
0220           ALAM2=ALAM3**2
0221         ENDIF
0222 C...Divide Lambda by PARP(64) (equivalent to mult pT2 by PARP(64))
0223         ALAM2=ALAM2/PARP(64)
0224 C...Overestimated ZMAX:
0225         IF (MQMASS.EQ.0) THEN
0226 C...Massless
0227           ZMAX=1D0-0.5D0*(PT2MNE/SHTNOW(MI))*(SQRT(1D0+4D0*SHTNOW(MI)
0228      &         /PT2MNE)-1D0)
0229         ELSE
0230 C...Massive (limit for bremsstrahlung diagram > creation)
0231           FMQ=SQRT(RMQ2/SHTNOW(MI))
0232           ZMAX=1D0/(1D0+FMQ)
0233         ENDIF
0234         ZMIN=XB/XMXC
0235
0236 C...If kinematically impossible then do not evolve.
0237         IF(PT2.LT.PT2CUT.OR.ZMAX.LE.ZMIN) RETURN
0238
0239 C...Reset Altarelli-Parisi and PDF weights.
0240         DO 140 KFL=-5,5
0241           WTAP(KFL)=0D0
0242           WTPDF(KFL)=0D0
0243   140   CONTINUE
0244         WTAP(21)=0D0
0245         WTPDF(21)=0D0
0246 C...Zero joining weights and compute X(partner) and X(mother) values.
0247         IF (MSTP(96).NE.0) THEN
0248           NJN=0
0249           DO 150 MJ=1,MINT(31)
0250             WTAPJ(MJ)=0D0
0251             WTPDFJ(MJ)=0D0
0252             X1(MJ)=XMI(JS,MJ)/(VINT(142+JS)+XMI(JS,MJ))
0253             Y(MJ)=(XMI(JS,MI)+XMI(JS,MJ))/(VINT(142+JS)+XMI(JS,MJ)
0254      &           +XMI(JS,MI))
0255   150     CONTINUE
0256         ENDIF
0257
0258 C...Approximate Altarelli-Parisi weights (integrated AP dz).
0259 C...q -> q, g -> q or q -> q + gamma (already set which).
0260         IF(KFLBA.LE.5) THEN
0261 C...Val and cmp quarks get an extra sqrt(z) to smooth their bumps.
0262           IF (KSVCB.LT.0) THEN
0263             WTAP(KFLB)=(8D0/3D0)*LOG((1D0-ZMIN)/(1D0-ZMAX))
0264           ELSE
0265             RMIN=(1+SQRT(ZMIN))/(1-SQRT(ZMIN))
0266             RMAX=(1+SQRT(ZMAX))/(1-SQRT(ZMAX))
0267             WTAP(KFLB)=(8D0/3D0)*LOG(RMAX/RMIN)
0268           ENDIF
0269           WTAP(21)=0.5D0*(ZMAX-ZMIN)
0270           WTAPE=(2D0/9D0)*LOG((1D0-ZMIN)/(1D0-ZMAX))
0271           IF(MOD(KFLBA,2).EQ.0) WTAPE=4D0*WTAPE
0272           IF(MECOR.GE.1.AND.NISGEN(JS,MI).EQ.0) THEN
0273             WTAP(KFLB)=WTFF*WTAP(KFLB)
0274             WTAP(21)=WTGF*WTAP(21)
0275             WTAPE=WTFF*WTAPE
0276           ENDIF
0277           IF (KSVCB.GE.1) THEN
0278 C...Kill normal creation but add joining diagrams for cmp quark.
0279             WTAP(21)=0D0
0280             IF (KFLBA.EQ.4.OR.KFLBA.EQ.5) THEN
0281               CALL PYERRM(9,'(PYPTIS:) Sorry, I got a heavy companion'//
0282      &             " quark here. Not handled yet, giving up!")
0283               PT2=0D0
0284               MINT(51)=1
0285               RETURN
0286             ENDIF
0287 C...Check for possible joinings
0288             IF (MSTP(96).NE.0.AND.MJOIND(JS,MI).EQ.0) THEN
0289 C...Find companion's companion.
0290               MJ=0
0291   160         MJ=MJ+1
0292               IF (IMI(JS,MJ,2).NE.IMI(JS,MI,1)) GOTO 160
0293               IF (MJOIND(JS,MJ).EQ.0) THEN
0294                 Y(MI)=YB+YS
0295                 Z=YB/Y(MI)
0296                 WTAPJ(MJ)=Z*(1D0-Z)*0.5D0*(Z**2+(1D0-Z)**2)
0297                 IF (WTAPJ(MJ).GT.1D-6) THEN
0298                   NJN=1
0299                 ELSE
0300                   WTAPJ(MJ)=0D0
0301                 ENDIF
0302               ENDIF
0303 C...Add trial gluon joinings.
0304               DO 170 MJ=1,MINT(31)
0305                 KFLC=K(IMI(JS,MJ,1),2)
0306                 IF (KFLC.NE.21.OR.MJOIND(JS,MJ).NE.0) GOTO 170
0307                 Z=XMI(JS,MJ)/(XMI(JS,MI)+XMI(JS,MJ))
0308                 WTAPJ(MJ)=6D0*(Z**2+(1D0-Z)**2)
0309                 IF (WTAPJ(MJ).GT.1D-6) THEN
0310                   NJN=NJN+1
0311                 ELSE
0312                   WTAPJ(MJ)=0D0
0313                 ENDIF
0314   170         CONTINUE
0315             ENDIF
0316           ELSEIF (IMI(JS,MI,2).GE.0) THEN
0317 C...Kill creation diagram for val quarks and sea quarks with companions.
0318             WTAP(21)=0D0
0319           ELSEIF (MQMASS.EQ.0) THEN
0320 C...Extra safety factor for massless sea quark creation.
0321             WTAP(21)=WTAP(21)*1.25D0
0322           ENDIF
0323
0324 C...  q -> g, g -> g.
0325         ELSEIF(KFLB.EQ.21) THEN
0326 C...Here we decide later whether a quark picked up is valence or
0327 C...sea, so we maintain the extra factor sqrt(z) since we deal
0328 C...with the *sum* of sea and valence in this context.
0329           WTAPQ=(16D0/3D0)*(SQRT(1D0/ZMIN)-SQRT(1D0/ZMAX))
0330 C...new: do not allow backwards evol to pick up heavy flavour.
0331           DO 180 KFL=1,MIN(3,MSTP(58))
0332             WTAP(KFL)=WTAPQ
0333             WTAP(-KFL)=WTAPQ
0334   180     CONTINUE
0335           WTAP(21)=6D0*LOG(ZMAX*(1D0-ZMIN)/(ZMIN*(1D0-ZMAX)))
0336           IF(MECOR.GE.1.AND.NISGEN(JS,MI).EQ.0) THEN
0337             WTAPQ=WTFG*WTAPQ
0338             WTAP(21)=WTGG*WTAP(21)
0339           ENDIF
0340 C...Check for possible joinings (companions handled separately above)
0341           IF (MSTP(96).NE.0.AND.MINT(31).GE.2.AND.MJOIND(JS,MI).EQ.0)
0342      &         THEN
0343             DO 190 MJ=1,MINT(31)
0344               IF (MJ.EQ.MI.OR.MJOIND(JS,MJ).NE.0) GOTO 190
0345               KSVCC=IMI(JS,MJ,2)
0346               IF (IMI(JS,MJ,2).GT.IMISEP(MJ)) KSVCC=-1
0347               IF (KSVCC.GE.1) GOTO 190
0348               KFLC=K(IMI(JS,MJ,1),2)
0349 C...Only try g -> g + g once.
0350               IF (MJ.GT.MI.AND.KFLC.EQ.21) GOTO 190
0351               Z=XMI(JS,MJ)/(XMI(JS,MI)+XMI(JS,MJ))
0352               IF (KFLC.EQ.21) THEN
0353                 WTAPJ(MJ)=6D0*(Z**2+(1D0-Z)**2)
0354               ELSE
0355                 WTAPJ(MJ)=Z*4D0/3D0*(1D0+Z**2)
0356               ENDIF
0357               IF (WTAPJ(MJ).GT.1D-6) THEN
0358                 NJN=NJN+1
0359               ELSE
0360                 WTAPJ(MJ)=0D0
0361               ENDIF
0362   190       CONTINUE
0363           ENDIF
0364         ENDIF
0365
0366 C...Initialize massive quark evolution
0367         IF (MQMASS.NE.0) THEN
0368           RML=(RMQ2+VINT(18))/ALAM2
0369           TML=LOG(RML)
0370           TPL=LOG((PT2+VINT(18))/ALAM2)
0371           TPM=LOG((PT2+VINT(18))/RMQ2)
0372           WN=WTAP(21)*WPDF0/B0
0373         ENDIF
0374
0375
0376 C...Loopback point for iteration
0377         NTRY=0
0378         NTHRES=0
0379   200   NTRY=NTRY+1
0380         IF(NTRY.GT.500) THEN
0381           CALL PYERRM(9,'(PYPTIS:) failed to evolve shower.')
0382           MINT(51)=1
0383           RETURN
0384         ENDIF
0385
0386 C...  Calculate PDF weights and sum for evolution rate.
0387         WTSUM=0D0
0388         XFBO=MAX(1D-10,XFB(KFLB))
0389         DO 210 KFL=-5,5
0390           WTPDF(KFL)=XFB(KFL)/XFBO
0391           WTSUM=WTSUM+WTAP(KFL)*WTPDF(KFL)
0392   210   CONTINUE
0393 C...Only add gluon mother diagram for massless KFLB.
0394         IF(MQMASS.EQ.0) THEN
0395           WTPDF(21)=XFB(21)/XFBO
0396           WTSUM=WTSUM+WTAP(21)*WTPDF(21)
0397         ENDIF
0398         WTSUM=MAX(0.0001D0,WTSUM)
0399         WTSUMS=WTSUM
0400 C...Add joining diagrams where applicable.
0401         WTJOIN=0D0
0402         IF (MSTP(96).NE.0.AND.NJN.NE.0) THEN
0403           DO 220 MJ=1,MINT(31)
0404             IF (WTAPJ(MJ).LT.1D-3) GOTO 220
0405             WTPDFJ(MJ)=1D0/XFBO
0406 C...x and x*pdf (+ sea/val) for parton C.
0407             KFLC=K(IMI(JS,MJ,1),2)
0408             KFLCA=IABS(KFLC)
0409             KSVCC=MAX(-1,IMI(JS,MJ,2))
0410             IF (IMI(JS,MJ,2).GT.IMISEP(MJ)) KSVCC=-1
0411             MINT(30)=JS
0412             MINT(36)=MJ
0413             CALL PYPDFU(KFBEAM(JS),X1(MJ),PT2,XFJ)
0414             MINT(36)=MI
0415             IF (KFLCA.LE.6.AND.KSVCC.LE.0) THEN
0416               XFJ(KFLC)=XPSVC(KFLC,KSVCC)
0417             ELSEIF (KSVCC.GE.1) THEN
0418               print*, 'error! parton C is companion!'
0419             ENDIF
0420             WTPDFJ(MJ)=WTPDFJ(MJ)/XFJ(KFLC)
0421 C...x and x*pdf (+ sea/val) for parton A.
0422             KFLA=21
0423             KSVCA=0
0424             IF (KFLCA.EQ.21.AND.KFLBA.LE.5) THEN
0425               KFLA=KFLB
0426               KSVCA=KSVCB
0427             ELSEIF (KFLBA.EQ.21.AND.KFLCA.LE.5) THEN
0428               KFLA=KFLC
0429               KSVCA=KSVCC
0430             ENDIF
0431             MINT(30)=JS
0432             IF (KSVCA.LE.0) THEN
0433 C...Consider C the "evolved" parton if B is gluon. Val/sea
0434 C...counting will then be done correctly in PYPDFU.
0435               IF (KFLBA.EQ.21) MINT(36)=MJ
0436               CALL PYPDFU(KFBEAM(JS),Y(MJ),PT2,XFJ)
0437               MINT(36)=MI
0438               IF (IABS(KFLA).LE.6) XFJ(KFLA)=XPSVC(KFLA,KSVCA)
0439             ELSE
0440 C...If parton A is companion, use Y(MI) and YS in call to PYFCMP.
0441               XFJ(KFLA)=PYFCMP(Y(MI)/VINT(140),YS/VINT(140),MSTP(87))
0442             ENDIF
0443             WTPDFJ(MJ)=XFJ(KFLA)*WTPDFJ(MJ)
0444             WTJOIN=WTJOIN+WTAPJ(MJ)*WTPDFJ(MJ)
0445   220     CONTINUE
0446         ENDIF
0447
0448 C...Pick normal pT2 (in overestimated z range).
0449   230   PT2OLD=PT2
0450         WTSUM=WTSUMS
0451         PT2=ALAM2*((PT2+VINT(18))/ALAM2)**(PYR(0)**(B0/WTSUM))-VINT(18)
0452         KFLC=21
0453
0454 C...Evolve q -> q gamma separately, pick it if larger pT.
0455         IF(KFLBA.LE.5) THEN
0456           PT2QED=(PT2OLD+VINT(18))*PYR(0)**(1D0/(AEM2PI*WTAPE))-VINT(18)
0457           IF(PT2QED.GT.PT2) THEN
0458             PT2=PT2QED
0459             KFLC=22
0460             KFLA=KFLB
0461           ENDIF
0462         ENDIF
0463
0464 C...  Evolve massive quark creation separately.
0465         MCRQQ=0
0466         IF (MQMASS.NE.0) THEN
0467           PT2CR=(RMQ2+VINT(18))*(RML**(TPM/(TPL*PYR(0)**(-TML/WN)-TPM)))
0468      &         -VINT(18)
0469 C...  Ensure mininimum PT2CR and force creation near threshold.
0470           IF (PT2CR.LT.TMIN*RMQ2) THEN
0471             NTHRES=NTHRES+1
0472             IF (NTHRES.GT.50) THEN
0473               CALL PYERRM(9,'(PYPTIS:) no phase space left for '//
0474      &             'massive quark creation. Gave up trying.')
0475               MINT(51)=1
0476               RETURN
0477             ENDIF
0478             PT2=0D0
0479             PT2CR=TMIN*RMQ2
0480             MCRQQ=2
0481           ENDIF
0482 C...  Select largest PT2 (brems or creation):
0483           IF (PT2CR.GT.PT2) THEN
0484             MCRQQ=MAX(MCRQQ,1)
0485             WTSUM=0D0
0486             PT2=PT2CR
0487             KFLA=21
0488           ELSE
0489             MCRQQ=0
0490             KFLA=KFLB
0491           ENDIF
0492 C...  Compute logarithms for this PT2
0493           TPL=LOG((PT2+VINT(18))/ALAM2)
0494           TPM=LOG((PT2+VINT(18))/(RMQ2+VINT(18)))
0495           WTCRQQ=TPM/LOG(PT2/RMQ2)
0496         ENDIF
0497
0498 C...Evolve joining separately
0499         MJOIN=0
0500         IF (MSTP(96).NE.0.AND.NJN.NE.0) THEN
0501           PT2JN=ALAM2*((PT2OLD+VINT(18))/ALAM2)**(PYR(0)**(B0/WTJOIN))
0502      &         -VINT(18)
0503           IF (PT2JN.GE.PT2) THEN
0504             MJOIN=1
0505             PT2=PT2JN
0506           ENDIF
0507         ENDIF
0508
0509 C...Loopback if crossed c/b mass thresholds.
0510         IF(IZRG.EQ.3.AND.PT2.LT.RMB2) THEN
0511           PT2=RMB2
0512          GOTO 130
0513         ELSEIF(IZRG.EQ.2.AND.PT2.LT.RMC2) THEN
0514           PT2=RMC2
0515           GOTO 130
0516         ENDIF
0517
0518 C...Speed up shower. Skip if higher-PT acceptable branching
0519 C...already found somewhere else.
0520 C...Also finish if below lower cutoff.
0521
0522         IF (PT2.LT.PT2MX.OR.PT2.LT.PT2CUT) RETURN
0523
0524 C...Select parton A flavour (massive Q handled above.)
0525         IF (MQMASS.EQ.0.AND.KFLC.NE.22.AND.MJOIN.EQ.0) THEN
0526           WTRAN=PYR(0)*WTSUM
0527           KFLA=-6
0528   240     KFLA=KFLA+1
0529           WTRAN=WTRAN-WTAP(KFLA)*WTPDF(KFLA)
0530           IF(KFLA.LE.5.AND.WTRAN.GT.0D0) GOTO 240
0531           IF(KFLA.EQ.6) KFLA=21
0532         ELSEIF (MJOIN.EQ.1) THEN
0533 C...Tentative joining accept/reject.
0534           WTRAN=PYR(0)*WTJOIN
0535           MJ=0
0536   250     MJ=MJ+1
0537           WTRAN=WTRAN-WTAPJ(MJ)*WTPDFJ(MJ)
0538           IF(MJ.LE.MINT(31)-1.AND.WTRAN.GT.0D0) GOTO 250
0539           IF(MJOIND(JS,MJ).NE.0.OR.MJOIND(JS,MI).NE.0) THEN
0540             CALL PYERRM(9,'(PYPTIS:) Attempted double joining.'//
0541      &           ' Rejected.')
0542             GOTO 230
0543           ENDIF
0544 C...x*pdf (+ sea/val) at new pT2 for parton B.
0545           IF (KSVCB.LE.0) THEN
0546             MINT(30)=JS
0547             CALL PYPDFU(KFBEAM(JS),XB,PT2,XFB)
0548             IF (KFLBA.LE.6) XFB(KFLB)=XPSVC(KFLB,KSVCB)
0549           ELSE
0550 C...Companion distributions do not evolve.
0551             XFB(KFLB)=XFBO
0552           ENDIF
0553           WTVETO=1D0/WTPDFJ(MJ)/XFB(KFLB)
0554           KFLC=K(IMI(JS,MJ,1),2)
0555           KFLCA=IABS(KFLC)
0556           KSVCC=MAX(-1,IMI(JS,MJ,2))
0557           IF (KSVCB.GE.1) KSVCC=-1
0558 C...x*pdf (+ sea/val) at new pT2 for parton C.
0559           MINT(30)=JS
0560           MINT(36)=MJ
0561           CALL PYPDFU(KFBEAM(JS),X1(MJ),PT2,XFJ)
0562           MINT(36)=MI
0563           IF (KFLCA.LE.6.AND.KSVCC.LE.0) XFJ(KFLC)=XPSVC(KFLC,KSVCC)
0564           WTVETO=WTVETO/XFJ(KFLC)
0565 C...x and x*pdf (+ sea/val) at new pT2 for parton A.
0566           KFLA=21
0567           KSVCA=0
0568           IF (KFLCA.EQ.21.AND.KFLBA.LE.5) THEN
0569             KFLA=KFLB
0570             KSVCA=KSVCB
0571           ELSEIF (KFLBA.EQ.21.AND.KFLCA.LE.5) THEN
0572             KFLA=KFLC
0573             KSVCA=KSVCC
0574           ENDIF
0575           IF (KSVCA.LE.0) THEN
0576             MINT(30)=JS
0577             IF (KFLB.EQ.21) MINT(36)=MJ
0578             CALL PYPDFU(KFBEAM(JS),Y(MJ),PT2,XFJ)
0579             MINT(36)=MI
0580             IF (IABS(KFLA).LE.6) XFJ(KFLA)=XPSVC(KFLA,KSVCA)
0581           ELSE
0582             XFJ(KFLA)=PYFCMP(Y(MJ)/VINT(140),YS/VINT(140),MSTP(87))
0583           ENDIF
0584           WTVETO=WTVETO*XFJ(KFLA)
0585 C...Monte Carlo veto.
0586           IF (WTVETO.LT.PYR(0)) GOTO 200
0587 C...If accept, save PT2 of this joining.
0588           IF (PT2.GT.PT2MX) THEN
0589             PT2MX=PT2
0590             JSMX=2+JS
0591             MJN1MX=MJ
0592             MJN2MX=MI
0593             WTAPJ(MJ)=0D0
0594             NJN=0
0595           ENDIF
0596 C...Exit and continue evolution.
0597           GOTO 380
0598         ENDIF
0599         KFLAA=IABS(KFLA)
0600
0601 C...Choose z value (still in overestimated range) and corrective weight.
0602 C...Unphysical z will be rejected below when Q2 has is computed.
0603         WTZ=0D0
0604
0605 C...Note: ME and MQ>0 give corrections to overall weights, not shapes.
0606 C...q -> q + g or q -> q + gamma (already set which).
0607         IF (KFLAA.LE.5.AND.KFLBA.LE.5) THEN
0608           IF (KSVCB.LT.0) THEN
0609             Z=1D0-(1D0-ZMIN)*((1D0-ZMAX)/(1D0-ZMIN))**PYR(0)
0610           ELSE
0611             ZFAC=RMIN*(RMAX/RMIN)**PYR(0)
0612             Z=((1-ZFAC)/(1+ZFAC))**2
0613           ENDIF
0614           WTZ=0.5D0*(1D0+Z**2)
0615 C...Massive weight correction.
0616           IF (KFLBA.GE.4) WTZ=WTZ-Z*(1D0-Z)**2*RMQ2/PT2
0617 C...Valence quark weight correction (extra sqrt)
0618           IF (KSVCB.GE.0) WTZ=WTZ*SQRT(Z)
0619
0620 C...q -> g + q.
0621 C...NB: MQ>0 not yet implemented. Forced absent above.
0622         ELSEIF (KFLAA.LE.5.AND.KFLB.EQ.21) THEN
0623           KFLC=KFLA
0624           Z=ZMAX/(1D0+PYR(0)*(SQRT(ZMAX/ZMIN)-1D0))**2
0625           WTZ=0.5D0*(1D0+(1D0-Z)**2)*SQRT(Z)
0626
0627 C...g -> q + qbar.
0628         ELSEIF (KFLA.EQ.21.AND.KFLBA.LE.5) THEN
0629           KFLC=-KFLB
0630           Z=ZMIN+PYR(0)*(ZMAX-ZMIN)
0631           WTZ=Z**2+(1D0-Z)**2
0632 C...Massive correction
0633           IF (MQMASS.NE.0) THEN
0634             WTZ=WTZ+2D0*Z*(1D0-Z)*RMQ2/PT2
0635 C...Extra safety margin for light sea quark creation
0636           ELSEIF (KSVCB.LT.0) THEN
0637             WTZ=WTZ/1.25D0
0638           ENDIF
0639
0640 C...g -> g + g.
0641         ELSEIF (KFLA.EQ.21.AND.KFLB.EQ.21) THEN
0642           KFLC=21
0643           Z=1D0/(1D0+((1D0-ZMIN)/ZMIN)*((1D0-ZMAX)*ZMIN/
0644      &         (ZMAX*(1D0-ZMIN)))**PYR(0))
0645           WTZ=(1D0-Z*(1D0-Z))**2
0646         ENDIF
0647
0648 C...Derive Q2 from pT2.
0649         Q2B=PT2/(1D0-Z)
0650         IF (KFLBA.GE.4) Q2B=Q2B-RMQ2
0651
0652 C...Loopback if outside allowed z range for given pT2.
0653         RM2C=PYMASS(KFLC)**2
0654         PT2ADJ=Q2B-Z*(SHTNOW(MI)+Q2B)*(Q2B+RM2C)/SHTNOW(MI)
0655         IF (PT2ADJ.LT.1D-6) GOTO 230
0656
0657 C...Loopback if nonordered in angle/rapidity.
0658         IF (MSTP(62).GE.3.AND.NISGEN(JS,MI).GE.1) THEN
0659           IF(PT2.GT.((1D0-Z)/(Z*(1D0-ZSAV(JS,MI))))**2*PT2SAV(JS,MI))
0660      &         GOTO 230
0661         ENDIF
0662
0663 C...Select phi angle of branching at random.
0664         PHI=PARU(2)*PYR(0)
0665
0666 C...Matrix-element corrections for some processes.
0667         IF (MECOR.GE.1.AND.NISGEN(JS,MI).EQ.0) THEN
0668           IF (KFLAA.LE.20.AND.KFLBA.LE.20) THEN
0669             CALL PYMEWT(MECOR,1,Q2B*SHAT/SHTNOW(MI),Z,PHI,WTME)
0670             WTZ=WTZ*WTME/WTFF
0671           ELSEIF((KFLA.EQ.21.OR.KFLA.EQ.22).AND.KFLBA.LE.20) THEN
0672             CALL PYMEWT(MECOR,2,Q2B*SHAT/SHTNOW(MI),Z,PHI,WTME)
0673             WTZ=WTZ*WTME/WTGF
0674           ELSEIF(KFLAA.LE.20.AND.(KFLB.EQ.21.OR.KFLB.EQ.22)) THEN
0675             CALL PYMEWT(MECOR,3,Q2B*SHAT/SHTNOW(MI),Z,PHI,WTME)
0676             WTZ=WTZ*WTME/WTFG
0677           ELSEIF(KFLA.EQ.21.AND.KFLB.EQ.21) THEN
0678             CALL PYMEWT(MECOR,4,Q2B*SHAT/SHTNOW(MI),Z,PHI,WTME)
0679             WTZ=WTZ*WTME/WTGG
0680           ENDIF
0681         ENDIF
0682
0683 C...Parton distributions at new pT2 but old x.
0684         MINT(30)=JS
0685         CALL PYPDFU(KFBEAM(JS),XB,PT2,XFN)
0686 C...Treat val and cmp separately
0687         IF (KFLBA.LE.6.AND.KSVCB.LE.0) XFN(KFLB)=XPSVC(KFLB,KSVCB)
0688         IF (KSVCB.GE.1)
0689      &       XFN(KFLB)=PYFCMP(YB/VINT(140),YS/VINT(140),MSTP(87))
0690         XFBN=XFN(KFLB)
0691         IF(XFBN.LT.1D-20) THEN
0692           IF(KFLA.EQ.KFLB) THEN
0693             WTAP(KFLB)=0D0
0694             GOTO 200
0695           ELSE
0696             XFBN=1D-10
0697             XFN(KFLB)=XFBN
0698           ENDIF
0699         ENDIF
0700         DO 260 KFL=-5,5
0701           XFB(KFL)=XFN(KFL)
0702   260   CONTINUE
0703         XFB(21)=XFN(21)
0704
0705 C...Parton distributions at new pT2 and new x.
0706         XA=XB/Z
0707         MINT(30)=JS
0708         CALL PYPDFU(KFBEAM(JS),XA,PT2,XFA)
0709         IF (KFLBA.LE.5.AND.KFLAA.LE.5) THEN
0710 C...q -> q + g: only consider respective sea, val, or cmp content.
0711           IF (KSVCB.LE.0) THEN
0712             XFA(KFLA)=XPSVC(KFLA,KSVCB)
0713           ELSE
0714             YA=XA*(1D0-YS)
0715             XFA(KFLB)=PYFCMP(YA/VINT(140),YS/VINT(140),MSTP(87))
0716           ENDIF
0717         ENDIF
0718         XFAN=XFA(KFLA)
0719         IF(XFAN.LT.1D-20) THEN
0720           GOTO 200
0721         ENDIF
0722
0723 C...If weighting fails continue evolution.
0724         WTTOT=0D0
0725         IF (MCRQQ.EQ.0) THEN
0726           WTPDFA=1D0/WTPDF(KFLA)
0727           WTTOT=WTZ*XFAN/XFBN*WTPDFA
0728         ELSEIF(MCRQQ.EQ.1) THEN
0729           WTPDFA=TPM/WPDF0
0730           WTTOT=WTCRQQ*WTZ*XFAN/XFBN*WTPDFA
0731           XBEST=TPM/TPM0*XQ0
0732         ELSEIF(MCRQQ.EQ.2) THEN
0733 C...Force massive quark creation.
0734           WTTOT=1D0
0735         ENDIF
0736
0737 C...Loop back if trial emission fails.
0738         IF(WTTOT.GE.0D0.AND.WTTOT.LT.PYR(0)) GOTO 200
0739         WTACC=((1D0+PT2)/(0.25D0+PT2))**2
0740         IF(WTTOT.LT.0D0) THEN
0741           WRITE(CHWT,'(1P,E12.4)') WTTOT
0742           CALL PYERRM(19,'(PYPTIS:) Weight '//CHWT//' negative')
0743         ELSEIF(WTTOT.GT.WTACC) THEN
0744           WRITE(CHWT,'(1P,E12.4)') WTTOT
0745           IF (PT2.GT.PTEMAX.OR.WTTOT.GE.WTEMAX) THEN
0746 C...Too high weight: write out as error, but do not update error counter.
0747             IF(MSTU(29).EQ.0) MSTU(23)=MSTU(23)-1
0748             CALL PYERRM(19,
0749      &         '(PYPTIS:) Weight '//CHWT//' above unity')
0750             IF (PT2.GT.PTEMAX) PTEMAX=PT2
0751             IF (WTTOT.GT.WTEMAX) WTEMAX=WTTOT
0752           ELSE
0753             CALL PYERRM(9,
0754      &         '(PYPTIS:) Weight '//CHWT//' above unity')
0755           ENDIF
0756 C...Useful for debugging but commented out for distribution:
0757 C          print*, 'JS, MI',JS, MI
0758 C          print*, 'PT:',SQRT(PT2), ' MCRQQ',MCRQQ
0759 C          print*, 'A -> B C',KFLA, KFLB, KFLC
0760 C          XFAO=XFBO/WTPDFA
0761 C          print*, 'WT(Z,XFA,XFB)',WTZ, XFAN/XFAO, XFBO/XFBN
0762         ENDIF
0763
0764 C...Save acceptable branching.
0765         IF(PT2.GT.PT2MX) THEN
0766           MIMX=MINT(36)
0767           JSMX=JS
0768           PT2MX=PT2
0769           KFLAMX=KFLA
0770           KFLCMX=KFLC
0771           RM2CMX=RM2C
0772           Q2BMX=Q2B
0773           ZMX=Z
0774           PT2AMX=PT2ADJ
0775           PHIMX=PHI
0776         ENDIF
0777
0778 C----------------------------------------------------------------------
0779 C...MODE= 1: Accept stored shower branching. Update event record etc.
0780       ELSEIF (MODE.EQ.1) THEN
0781         MI=MIMX
0782         JS=JSMX
0783         SHAT=SHTNOW(MI)
0784         SIDE=3D0-2D0*JS
0785 C...Shift down rest of event record to make room for insertion.
0786         IT=IMISEP(MI)+1
0787         IM=IT+1
0788         IS=IMI(JS,MI,1)
0789         DO 280 I=N,IT,-1
0790           IF (K(I,3).GE.IT) K(I,3)=K(I,3)+2
0791           KT1=K(I,4)/MSTU(5)**2
0792           KT2=K(I,5)/MSTU(5)**2
0793           ID1=MOD(K(I,4),MSTU(5))
0794           ID2=MOD(K(I,5),MSTU(5))
0795           IM1=MOD(K(I,4)/MSTU(5),MSTU(5))
0796           IM2=MOD(K(I,5)/MSTU(5),MSTU(5))
0797           IF (ID1.GE.IT) ID1=ID1+2
0798           IF (ID2.GE.IT) ID2=ID2+2
0799           IF (IM1.GE.IT) IM1=IM1+2
0800           IF (IM2.GE.IT) IM2=IM2+2
0801           K(I,4)=KT1*MSTU(5)**2+IM1*MSTU(5)+ID1
0802           K(I,5)=KT2*MSTU(5)**2+IM2*MSTU(5)+ID2
0803           DO 270 IX=1,5
0804             K(I+2,IX)=K(I,IX)
0805             P(I+2,IX)=P(I,IX)
0806             V(I+2,IX)=V(I,IX)
0807   270     CONTINUE
0808           MCT(I+2,1)=MCT(I,1)
0809           MCT(I+2,2)=MCT(I,2)
0810   280   CONTINUE
0811         N=N+2
0812 C...Also update shifted-down pointers in IMI, IMISEP, and IPART.
0813         DO 290 JI=1,MINT(31)
0814           IF (IMI(1,JI,1).GE.IT) IMI(1,JI,1)=IMI(1,JI,1)+2
0815           IF (IMI(1,JI,2).GE.IT) IMI(1,JI,2)=IMI(1,JI,2)+2
0816           IF (IMI(2,JI,1).GE.IT) IMI(2,JI,1)=IMI(2,JI,1)+2
0817           IF (IMI(2,JI,2).GE.IT) IMI(2,JI,2)=IMI(2,JI,2)+2
0818           IF (JI.GE.MI) IMISEP(JI)=IMISEP(JI)+2
0819 C...Also update companion pointers to the present mother.
0820           IF (IMI(JS,JI,2).EQ.IS) IMI(JS,JI,2)=IM
0821   290   CONTINUE
0822         DO 300 IFS=1,NPART
0823           IF (IPART(IFS).GE.IT) IPART(IFS)=IPART(IFS)+2
0824   300   CONTINUE
0825 C...Zero entries dedicated for new timelike and mother partons.
0826         DO 320 I=IT,IT+1
0827           DO 310 J=1,5
0828             K(I,J)=0
0829             P(I,J)=0D0
0830             V(I,J)=0D0
0831   310     CONTINUE
0832           MCT(I,1)=0
0833           MCT(I,2)=0
0834   320   CONTINUE
0835
0836 C...Define timelike and new mother partons. History.
0837         K(IT,1)=3
0838         K(IT,2)=KFLCMX
0839         K(IM,1)=14
0840         K(IM,2)=KFLAMX
0841         K(IS,3)=IM
0842         K(IT,3)=IM
0843 C...Set mother origin = side.
0844         K(IM,3)=MINT(83)+JS+2
0845         IF(MI.GE.2) K(IM,3)=MINT(83)+JS
0846
0847 C...Define colour flow of branching.
0848         IM1=IM
0849         IM2=IM
0850 C...q -> q + gamma.
0851         IF(K(IT,2).EQ.22) THEN
0852           K(IT,1)=1
0853           ID1=IS
0854           ID2=IS
0855 C...q -> q + g.
0856         ELSEIF(K(IM,2).GT.0.AND.K(IM,2).LE.5.AND.K(IT,2).EQ.21) THEN
0857           ID1=IT
0858           ID2=IS
0859 C...q -> g + q.
0860         ELSEIF(K(IM,2).GT.0.AND.K(IM,2).LE.5) THEN
0861           ID1=IS
0862           ID2=IT
0863 C...qbar -> qbar + g.
0864         ELSEIF(K(IM,2).LT.0.AND.K(IM,2).GE.-5.AND.K(IT,2).EQ.21) THEN
0865           ID1=IS
0866           ID2=IT
0867 C...qbar -> g + qbar.
0868         ELSEIF(K(IM,2).LT.0.AND.K(IM,2).GE.-5) THEN
0869           ID1=IT
0870           ID2=IS
0871 C...g -> g + g; g -> q + qbar..
0872         ELSEIF((K(IT,2).EQ.21.AND.PYR(0).GT.0.5D0).OR.K(IT,2).LT.0) THEN
0873           ID1=IS
0874           ID2=IT
0875         ELSE
0876           ID1=IT
0877           ID2=IS
0878         ENDIF
0879         IF(IM1.EQ.IM) K(IM1,4)=K(IM1,4)+ID1
0880         IF(IM2.EQ.IM) K(IM2,5)=K(IM2,5)+ID2
0881         K(ID1,4)=K(ID1,4)+MSTU(5)*IM1
0882         K(ID2,5)=K(ID2,5)+MSTU(5)*IM2
0883         IF(ID1.NE.ID2) THEN
0884           K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
0885           K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
0886         ENDIF
0887         IF(K(IT,1).EQ.1) THEN
0888           K(IT,4)=0
0889           K(IT,5)=0
0890         ENDIF
0891 C...Update IMI and colour tag arrays.
0892         IMI(JS,MI,1)=IM
0893         DO 330 MC=1,2
0894           MCT(IT,MC)=0
0895           MCT(IM,MC)=0
0896   330   CONTINUE
0897         DO 340 JCS=4,5
0898           KCS=JCS
0899 C...If mother flag not yet set for spacelike parton, trace it.
0900           IF (K(IS,KCS)/MSTU(5)**2.LE.1) CALL PYCTTR(IS,-KCS,IM)
0901           IF(MINT(51).NE.0) RETURN
0902   340   CONTINUE
0903         DO 350 JCS=4,5
0904           KCS=JCS
0905 C...If mother flag not yet set for timelike parton, trace it.
0906           IF (K(IT,KCS)/MSTU(5)**2.LE.1) CALL PYCTTR(IT,KCS,IM)
0907           IF(MINT(51).NE.0) RETURN
0908   350   CONTINUE
0909
0910 C...Boost recoiling parton to compensate for Q2 scale.
0911         BETAZ=SIDE*(1D0-(1D0+Q2BMX/SHAT)**2)/
0912      &  (1D0+(1D0+Q2BMX/SHAT)**2)
0913         IR=IMI(3-JS,MI,1)
0914         CALL PYROBO(IR,IR,0D0,0D0,0D0,0D0,BETAZ)
0915
0916 C...Define system to be rotated and boosted
0917 C...(not including the 2 just added partons)
0918 C...(but including the docu lines for first interaction)
0919         IMIN=IMISEP(MI-1)+1
0920         IF (MI.EQ.1) IMIN=MINT(83)+5
0921         IMAX=IMISEP(MI)-2
0922
0923 C...Rotate back system in phi to compensate for subsequent rotation.
0924         CALL PYROBO(IMIN,IMAX,0D0,-PHIMX,0D0,0D0,0D0)
0925
0926 C...Define kinematics of new partons in old frame.
0927         IMAX=IMISEP(MI)
0928         P(IM,1)=SQRT(PT2AMX)*SHAT/(ZMX*(SHAT+Q2BMX))
0929         P(IM,3)=0.5D0*SQRT(SHAT)*((SHAT-Q2BMX)/((SHAT
0930      &       +Q2BMX)*ZMX)+(Q2BMX+RM2CMX)/SHAT)*SIDE
0931         P(IM,4)=SQRT(P(IM,1)**2+P(IM,3)**2)
0932         P(IT,1)=P(IM,1)
0933         P(IT,3)=P(IM,3)-0.5D0*(SHAT+Q2BMX)/SQRT(SHAT)*SIDE
0934         P(IT,4)=SQRT(P(IT,1)**2+P(IT,3)**2+RM2CMX)
0935         P(IT,5)=SQRT(RM2CMX)
0936
0937 C...Update internal line, now spacelike
0938         P(IS,1)=P(IM,1)-P(IT,1)
0939         P(IS,2)=P(IM,2)-P(IT,2)
0940         P(IS,3)=P(IM,3)-P(IT,3)
0941         P(IS,4)=P(IM,4)-P(IT,4)
0942         P(IS,5)=P(IS,4)**2-P(IS,1)**2-P(IS,2)**2-P(IS,3)**2
0943 C...Represent spacelike virtualities as -sqrt(abs(Q2)) .
0944         IF (P(IS,5).LT.0D0) THEN
0945           P(IS,5)=-SQRT(ABS(P(IS,5)))
0946         ELSE
0947           P(IS,5)=SQRT(P(IS,5))
0948         ENDIF
0949
0950 C...Boost entire system and rotate to new frame.
0951 C...(including docu lines)
0952         BETAX=(P(IM,1)+P(IR,1))/(P(IM,4)+P(IR,4))
0953         BETAZ=(P(IM,3)+P(IR,3))/(P(IM,4)+P(IR,4))
0954         IF(BETAX**2+BETAZ**2.GE.1D0) THEN
0955           CALL PYERRM(1,'(PYPTIS:) boost bigger than unity')
0956           MINT(51)=1
0957           IFAIL=-1
0958           RETURN
0959         ENDIF
0960         CALL PYROBO(IMIN,IMAX,0D0,0D0,-BETAX,0D0,-BETAZ)
0961         I1=IMI(1,MI,1)
0962         THETA=PYANGL(P(I1,3),P(I1,1))
0963         CALL PYROBO(IMIN,IMAX,-THETA,PHIMX,0D0,0D0,0D0)
0964
0965 C...Global statistics.
0966         MINT(352)=MINT(352)+1
0967         VINT(352)=VINT(352)+SQRT(P(IT,1)**2+P(IT,2)**2)
0968         IF (MINT(352).EQ.1) VINT(357)=SQRT(P(IT,1)**2+P(IT,2)**2)
0969
0970 C...Add parton with relevant pT scale for timelike shower.
0971         IF (K(IT,2).NE.22) THEN
0972           NPART=NPART+1
0973           IPART(NPART)=IT
0974           PTPART(NPART)=SQRT(PT2AMX)
0975         ENDIF
0976
0977 C...Update saved variables.
0978         SHTNOW(MIMX)=SHTNOW(MIMX)/ZMX
0979         NISGEN(JSMX,MIMX)=NISGEN(JSMX,MIMX)+1
0980         XMI(JSMX,MIMX)=XMI(JSMX,MIMX)/ZMX
0981         PT2SAV(JSMX,MIMX)=PT2MX
0982         ZSAV(JS,MIMX)=ZMX
0983
0984         KSA=IABS(K(IS,2))
0985         KMA=IABS(K(IM,2))
0986         IF (KSA.EQ.21.AND.KMA.GE.1.AND.KMA.LE.5) THEN
0987 C...Gluon reconstructs to quark.
0988 C...Decide whether newly created quark is valence or sea:
0989           MINT(30)=JS
0990           CALL PYPTMI(2,PT2NOW,PTDUM1,PTDUM2,IFAIL)
0991           IF(MINT(51).NE.0) RETURN
0992         ENDIF
0993         IF(KSA.GE.1.AND.KSA.LE.5.AND.KMA.EQ.21) THEN
0994 C...Quark reconstructs to gluon.
0995 C...Now some guy may have lost his companion. Check.
0996           ICMP=IMI(JS,MI,2)
0997           IF (ICMP.GT.0) THEN
0998             CALL PYERRM(9,'(PYPTIS:) Sorry, companion quark radiated'
0999      &           //' away. Cannot handle that yet. Giving up.')
1000             MINT(51)=1
1001             RETURN
1002           ELSEIF(ICMP.LT.0) THEN
1003 C...A sea quark with companion still in BR was reconstructed to a gluon.
1004 C...Companion should now be removed from the beam remnant.
1005 C...(Momentum integral is automatically updated in next call to PYPDFU.)
1006             ICMP=-ICMP
1007             IFL=-K(IS,2)
1008             DO 370 JCMP=ICMP,NVC(JS,IFL)-1
1009               XASSOC(JS,IFL,JCMP)=XASSOC(JS,IFL,JCMP+1)
1010               DO 360 JI=1,MINT(31)
1011                 KMI=-IMI(JS,JI,2)
1012                 JFL=-K(IMI(JS,JI,1),2)
1013                 IF (KMI.EQ.JCMP+1.AND.JFL.EQ.IFL) IMI(JS,JI,2)=IMI(JS,JI
1014      &               ,2)+1
1015   360         CONTINUE
1016   370       CONTINUE
1017             NVC(JS,IFL)=NVC(JS,IFL)-1
1018           ENDIF
1019 C...Set gluon IMI(JS,MI,2) = 0.
1020           IMI(JS,MI,2)=0
1021         ELSEIF(KSA.GE.1.AND.KSA.LE.5.AND.KMA.NE.21) THEN
1022 C...Quark reconstructing to quark. If sea with companion still in BR
1023 C...then update associated x value.
1024 C...(Momentum integral is automatically updated in next call to PYPDFU.)
1025           IF (IMI(JS,MI,2).LT.0) THEN
1026             ICMP=-IMI(JS,MI,2)
1027             IFL=-K(IS,2)
1028             XASSOC(JS,IFL,ICMP)=XMI(JSMX,MIMX)
1029           ENDIF
1030         ENDIF
1031
1032       ENDIF
1033
1034 C...If reached this point, normal exit.
1035   380 IFAIL=0
1036
1037       RETURN
1038       END