pythiaeRHIC/pythia/pyptmi.f

0001
0002 C*********************************************************************
0003
0004 C...PYPTMI
0005 C...Handles the generation of additional interactions in the new
0006 C...multiple interactions framework.
0007 C...MODE=-1 : Initalize MI from scratch.
0008 C...MODE= 0 : Generate trial interaction. Start at PT2NOW, solve
0009 C...         Sudakov for PT2, abort if below PT2CUT.
0010 C...MODE= 1 : Accept interaction at PT2NOW and store variables.
0011 C...MODE= 2 : Decide sea/val/cmp for kicked-out quark at PT2NOW
0012 C...PT2NOW  : Starting (max) PT2 scale for evolution.
0013 C...PT2CUT  : Lower limit for evolution.
0014 C...PT2     : Result of evolution. Generated PT2 for trial interaction.
0015 C...IFAIL   : Status return code.
0016 C...         = 0: All is well.
0017 C...         < 0: Phase space exhausted, generation to be terminated.
0018 C...         > 0: Additional interaction vetoed, but continue evolution.
0019
0020       SUBROUTINE PYPTMI(MODE,PT2NOW,PT2CUT,PT2,IFAIL)
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/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0033       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0034       COMMON/PYINT1/MINT(400),VINT(400)
0035       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0036       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
0037       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
0038       COMMON/PYINT7/SIGT(0:6,0:6,0:5)
0039       COMMON/PYINTM/KFIVAL(2,3),NMI(2),IMI(2,800,2),NVC(2,-6:6),
0040      &     XASSOC(2,-6:6,240),XPSVC(-6:6,-1:240),PVCTOT(2,-1:1),
0041      &     XMI(2,240),PT2MI(240),IMISEP(0:240)
0042       COMMON/PYISMX/MIMX,JSMX,KFLAMX,KFLCMX,KFBEAM(2),NISGEN(2,240),
0043      &     PT2MX,PT2AMX,ZMX,RM2CMX,Q2BMX,PHIMX
0044       COMMON/PYCTAG/NCT,MCT(4000,2)
0045 C...Local arrays and saved variables.
0046       DIMENSION WDTP(0:400),WDTE(0:400,0:5),XPQ(-25:25)
0047
0048       SAVE /PYPART/,/PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYPARS/,
0049      &     /PYINT1/,/PYINT2/,/PYINT3/,/PYINT5/,/PYINT7/,/PYINTM/,
0050      &     /PYISMX/,/PYCTAG/
0051       SAVE XT2FAC,SIGS
0052
0053       IFAIL=0
0054 C...Set MI subprocess = QCD 2 -> 2.
0055       ISUB=96
0056
0057 C----------------------------------------------------------------------
0058 C...MODE=-1: Initialize from scratch
0059       IF (MODE.EQ.-1) THEN
0060 C...Initialize PT2 array.
0061         PT2MI(1)=VINT(54)
0062 C...Initialize list of incoming beams and partons from two sides.
0063         DO 110 JS=1,2
0064           DO 100 MI=1,240
0065             IMI(JS,MI,1)=0
0066             IMI(JS,MI,2)=0
0067   100     CONTINUE
0068           NMI(JS)=1
0069           IMI(JS,1,1)=MINT(84)+JS
0070           IMI(JS,1,2)=0
0071           XMI(JS,1)=VINT(40+JS)
0072 C...Rescale x values to fractions of photon energy.
0073           IF(MINT(18+JS).EQ.1) XMI(JS,1)=VINT(40+JS)/VINT(154+JS)
0074 C...Hard reset: hard interaction initiators motherless by definition.
0075           K(MINT(84)+JS,3)=2+JS
0076           K(MINT(84)+JS,4)=MOD(K(MINT(84)+JS,4),MSTU(5))
0077           K(MINT(84)+JS,5)=MOD(K(MINT(84)+JS,5),MSTU(5))
0078   110   CONTINUE
0079         IMISEP(0)=MINT(84)
0080         IMISEP(1)=N
0081         IF (MOD(MSTP(81),10).GE.1) THEN
0082           IF(MSTP(82).LE.1) THEN
0083             SIGRAT=XSEC(ISUB,1)/MAX(1D-10,VINT(315)*VINT(316)*SIGT(0,0
0084      &           ,5))
0085             IF(MINT(141).NE.0.OR.MINT(142).NE.0) SIGRAT=SIGRAT*
0086      &           VINT(317)/(VINT(318)*VINT(320))
0087             XT2FAC=SIGRAT*VINT(149)/(1D0-VINT(149))
0088           ELSE
0089             XT2FAC=VINT(146)*VINT(148)*XSEC(ISUB,1)/
0090      &           MAX(1D-10,SIGT(0,0,5))*VINT(149)*(1D0+VINT(149))
0091           ENDIF
0092         ENDIF
0093 C...Zero entries relating to scatterings beyond the first.
0094         DO 120 MI=2,240
0095           IMI(1,MI,1)=0
0096           IMI(2,MI,1)=0
0097           IMI(1,MI,2)=0
0098           IMI(2,MI,2)=0
0099           IMISEP(MI)=IMISEP(1)
0100           PT2MI(MI)=0D0
0101           XMI(1,MI)=0D0
0102           XMI(2,MI)=0D0
0103   120   CONTINUE
0104 C...Initialize factors for PDF reshaping.
0105         DO 140 JS=1,2
0106           KFBEAM(JS)=MINT(10+JS)
0107           IF(MINT(18+JS).EQ.1) KFBEAM(JS)=22
0108           KFABM=IABS(KFBEAM(JS))
0109           KFSBM=ISIGN(1,KFBEAM(JS))
0110
0111 C...Zero flavour content of incoming beam particle.
0112           KFIVAL(JS,1)=0
0113           KFIVAL(JS,2)=0
0114           KFIVAL(JS,3)=0
0115 C...  Flavour content of baryon.
0116           IF(KFABM.GT.1000) THEN
0117             KFIVAL(JS,1)=KFSBM*MOD(KFABM/1000,10)
0118             KFIVAL(JS,2)=KFSBM*MOD(KFABM/100,10)
0119             KFIVAL(JS,3)=KFSBM*MOD(KFABM/10,10)
0120 C...  Flavour content of pi+-, K+-.
0121           ELSEIF(KFABM.EQ.211) THEN
0122             KFIVAL(JS,1)=KFSBM*2
0123             KFIVAL(JS,2)=-KFSBM
0124           ELSEIF(KFABM.EQ.321) THEN
0125             KFIVAL(JS,1)=-KFSBM*3
0126             KFIVAL(JS,2)=KFSBM*2
0127 C...  Flavour content of pi0, gamma, K0S, K0L not defined yet.
0128           ENDIF
0129
0130 C...Zero initial valence and companion content.
0131           DO 130 IFL=-6,6
0132             NVC(JS,IFL)=0
0133   130     CONTINUE
0134   140   CONTINUE
0135 C...Set up colour line tags starting from hard interaction initiators.
0136         NCT=0
0137 C...Reset colour tag array and colour processing flags.
0138         DO 150 I=IMISEP(0)+1,N
0139           MCT(I,1)=0
0140           MCT(I,2)=0
0141           K(I,4)=MOD(K(I,4),MSTU(5)**2)
0142           K(I,5)=MOD(K(I,5),MSTU(5)**2)
0143   150   CONTINUE
0144 C...  Consider each side in turn.
0145         DO 170 JS=1,2
0146           I1=IMI(JS,1,1)
0147           I2=IMI(3-JS,1,1)
0148           DO 160 JCS=4,5
0149             IF (K(I1,2).NE.21.AND.(9-2*JCS).NE.ISIGN(1,K(I1,2)))
0150      &           GOTO 160
0151             IF (K(I1,JCS)/MSTU(5)**2.NE.0) GOTO 160
0152             KCS=JCS
0153             CALL PYCTTR(I1,KCS,I2)
0154             IF(MINT(51).NE.0) RETURN
0155   160     CONTINUE
0156   170   CONTINUE
0157
0158 C...Range checking for companion quark pdf large-x param.
0159         IF (MSTP(87).LT.0) THEN
0160           CALL PYERRM(19,'(PYPTMI:) MSTP(87) out of range. Forced'//
0161      &         ' MSTP(87)=0')
0162           MSTP(87)=0
0163         ELSEIF (MSTP(87).GT.4) THEN
0164           CALL PYERRM(19,'(PYPTMI:) MSTP(87) out of range. Forced'//
0165      &         ' MSTP(87)=4')
0166           MSTP(87)=4
0167         ENDIF
0168
0169 C----------------------------------------------------------------------
0170 C...MODE=0: Generate trial interaction. Return codes:
0171 C...IFAIL < 0: Phase space exhausted, generation to be terminated.
0172 C...IFAIL = 0: Additional interaction generated at PT2.
0173 C...IFAIL > 0: Additional interaction vetoed, but continue evolution.
0174       ELSEIF (MODE.EQ.0) THEN
0175 C...Abolute MI max scale = VINT(62)
0176         XT2=4D0*MIN(PT2NOW,VINT(62))/VINT(2)
0177   180   IF(MSTP(82).LE.1) THEN
0178           XT2=XT2FAC*XT2/(XT2FAC-XT2*LOG(PYR(0)))
0179           IF(XT2.LT.VINT(149)) IFAIL=-2
0180         ELSE
0181           IF(XT2.LE.0.01001D0*VINT(149)) THEN
0182             IFAIL=-3
0183           ELSE
0184             XT2=XT2FAC*(XT2+VINT(149))/(XT2FAC-(XT2+VINT(149))*
0185      &           LOG(PYR(0)))-VINT(149)
0186           ENDIF
0187         ENDIF
0188 C...Also exit if below lower limit or if higher trial branching
0189 C...already found.
0190         PT2=0.25D0*VINT(2)*XT2
0191         IF (PT2.LE.PT2CUT) IFAIL=-4
0192         IF (PT2.LE.PT2MX) IFAIL=-5
0193         IF (IFAIL.NE.0) THEN
0194           PT2=0D0
0195           RETURN
0196         ENDIF
0197         IF(MSTP(82).GE.2) PT2=MAX(0.25D0*VINT(2)*0.01D0*VINT(149),PT2)
0198         VINT(25)=4D0*PT2/VINT(2)
0199         XT2=VINT(25)
0200
0201 C...Choose tau and y*. Calculate cos(theta-hat).
0202         IF(PYR(0).LE.COEF(ISUB,1)) THEN
0203           TAUT=(2D0*(1D0+SQRT(1D0-XT2))/XT2-1D0)**PYR(0)
0204           TAU=XT2*(1D0+TAUT)**2/(4D0*TAUT)
0205         ELSE
0206           TAU=XT2*(1D0+TAN(PYR(0)*ATAN(SQRT(1D0/XT2-1D0)))**2)
0207         ENDIF
0208         VINT(21)=TAU
0209 C...New: require shat > 1.
0210         IF(TAU*VINT(2).LT.1D0) GOTO 180
0211         CALL PYKLIM(2)
0212         RYST=PYR(0)
0213         MYST=1
0214         IF(RYST.GT.COEF(ISUB,8)) MYST=2
0215         IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
0216         CALL PYKMAP(2,MYST,PYR(0))
0217         VINT(23)=SQRT(MAX(0D0,1D0-XT2/TAU))*(-1)**INT(1.5D0+PYR(0))
0218
0219 C...Check that x not used up. Accept or reject kinematical variables.
0220         X1M=SQRT(TAU)*EXP(VINT(22))
0221         X2M=SQRT(TAU)*EXP(-VINT(22))
0222         IF(VINT(143)-X1M.LT.0.01D0.OR.VINT(144)-X2M.LT.0.01D0) GOTO 180
0223         VINT(71)=0.5D0*VINT(1)*SQRT(XT2)
0224         CALL PYSIGH(NCHN,SIGS)
0225         IF(MINT(141).NE.0.OR.MINT(142).NE.0) SIGS=SIGS*VINT(320)
0226         IF(SIGS.LT.XSEC(ISUB,1)*PYR(0)) GOTO 180
0227         IF(MINT(141).NE.0.OR.MINT(142).NE.0) SIGS=SIGS/VINT(320)
0228
0229 C...Save if highest PT so far.
0230         IF (PT2.GT.PT2MX) THEN
0231           JSMX=0
0232           MIMX=MINT(31)+1
0233           PT2MX=PT2
0234         ENDIF
0235
0236 C----------------------------------------------------------------------
0237 C...MODE=1: Generate and save accepted scattering.
0238       ELSEIF (MODE.EQ.1) THEN
0239         PT2=PT2NOW
0240 C...Reset K, P, V, and MCT vectors.
0241         DO 200 I=N+1,N+4
0242           DO 190 J=1,5
0243             K(I,J)=0
0244             P(I,J)=0D0
0245             V(I,J)=0D0
0246   190     CONTINUE
0247           MCT(I,1)=0
0248           MCT(I,2)=0
0249   200   CONTINUE
0250
0251         NTRY=0
0252 C...Choose flavour of reacting partons (and subprocess).
0253   210   NTRY=NTRY+1
0254         IF (NTRY.GT.50) THEN
0255           CALL PYERRM(9,'(PYPTMI:) Unable to generate additional '
0256      &               //'interaction. Giving up!')
0257           MINT(51)=1
0258           RETURN
0259         ENDIF
0260         RSIGS=SIGS*PYR(0)
0261         DO 220 ICHN=1,NCHN
0262           KFL1=ISIG(ICHN,1)
0263           KFL2=ISIG(ICHN,2)
0264           ICONMI=ISIG(ICHN,3)
0265           RSIGS=RSIGS-SIGH(ICHN)
0266           IF(RSIGS.LE.0D0) GOTO 230
0267   220   CONTINUE
0268
0269 C...Reassign to appropriate process codes.
0270   230   ISUBMI=ICONMI/10
0271         ICONMI=MOD(ICONMI,10)
0272
0273 C...Choose new quark flavour for annihilation graphs
0274         IF(ISUBMI.EQ.12.OR.ISUBMI.EQ.53) THEN
0275           SH=VINT(21)*VINT(2)
0276           CALL PYWIDT(21,SH,WDTP,WDTE)
0277   240     RKFL=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))*PYR(0)
0278           DO 250 I=1,MDCY(21,3)
0279             KFLF=KFDP(I+MDCY(21,2)-1,1)
0280             RKFL=RKFL-(WDTE(I,1)+WDTE(I,2)+WDTE(I,4))
0281             IF(RKFL.LE.0D0) GOTO 260
0282   250     CONTINUE
0283   260     IF(ISUBMI.EQ.53.AND.ICONMI.LE.2) THEN
0284             IF(KFLF.GE.4) GOTO 240
0285           ELSEIF(ISUBMI.EQ.53.AND.ICONMI.LE.4) THEN
0286             KFLF=4
0287             ICONMI=ICONMI-2
0288           ELSEIF(ISUBMI.EQ.53) THEN
0289             KFLF=5
0290             ICONMI=ICONMI-4
0291           ENDIF
0292         ENDIF
0293
0294 C...Final state flavours and colour flow: default values
0295         JS=1
0296         KFL3=KFL1
0297         KFL4=KFL2
0298         KCC=20
0299         KCS=ISIGN(1,KFL1)
0300
0301         IF(ISUBMI.EQ.11) THEN
0302 C...f + f' -> f + f' (g exchange); th = (p(f)-p(f))**2
0303           KCC=ICONMI
0304           IF(KFL1*KFL2.LT.0) KCC=KCC+2
0305
0306         ELSEIF(ISUBMI.EQ.12) THEN
0307 C...f + fbar -> f' + fbar'; th = (p(f)-p(f'))**2
0308           KFL3=ISIGN(KFLF,KFL1)
0309           KFL4=-KFL3
0310           KCC=4
0311
0312         ELSEIF(ISUBMI.EQ.13) THEN
0313 C...f + fbar -> g + g; th arbitrary
0314           KFL3=21
0315           KFL4=21
0316           KCC=ICONMI+4
0317
0318         ELSEIF(ISUBMI.EQ.28) THEN
0319 C...f + g -> f + g; th = (p(f)-p(f))**2
0320           IF(KFL1.EQ.21) JS=2
0321           KCC=ICONMI+6
0322           IF(KFL1.EQ.21) KCC=KCC+2
0323           IF(KFL1.NE.21) KCS=ISIGN(1,KFL1)
0324           IF(KFL2.NE.21) KCS=ISIGN(1,KFL2)
0325
0326         ELSEIF(ISUBMI.EQ.53) THEN
0327 C...g + g -> f + fbar; th arbitrary
0328           KCS=(-1)**INT(1.5D0+PYR(0))
0329           KFL3=ISIGN(KFLF,KCS)
0330           KFL4=-KFL3
0331           KCC=ICONMI+10
0332
0333         ELSEIF(ISUBMI.EQ.68) THEN
0334 C...g + g -> g + g; th arbitrary
0335           KCC=ICONMI+12
0336           KCS=(-1)**INT(1.5D0+PYR(0))
0337         ENDIF
0338
0339 C...Check that massive sea quarks have non-zero phase space for g -> Q Q
0340         IF (IABS(KFL3).EQ.4.OR.IABS(KFL4).EQ.4.OR.IABS(KFL3).EQ.5
0341      &       .OR.IABS(KFL4).EQ.5) THEN
0342           RMMAX2=MAX(PMAS(PYCOMP(KFL3),1),PMAS(PYCOMP(KFL4),1))**2
0343           IF (PT2.LE.1.05*RMMAX2) THEN
0344             IF (NTRY.EQ.1) CALL PYERRM(9,'(PYPTMI:) Heavy quarks'
0345      &           //' created below threshold. Rejected.')
0346             GOTO 210
0347           ENDIF
0348         ENDIF
0349
0350 C...Store flavours of scattering.
0351         MINT(13)=KFL1
0352         MINT(14)=KFL2
0353         MINT(15)=KFL1
0354         MINT(16)=KFL2
0355         MINT(21)=KFL3
0356         MINT(22)=KFL4
0357
0358 C...Set flavours and mothers of scattering partons.
0359         K(N+1,1)=14
0360         K(N+2,1)=14
0361         K(N+3,1)=3
0362         K(N+4,1)=3
0363         K(N+1,2)=KFL1
0364         K(N+2,2)=KFL2
0365         K(N+3,2)=KFL3
0366         K(N+4,2)=KFL4
0367         K(N+1,3)=MINT(83)+1
0368         K(N+2,3)=MINT(83)+2
0369         K(N+3,3)=N+1
0370         K(N+4,3)=N+2
0371
0372 C...Store colour connection indices.
0373         DO 270 J=1,2
0374           JC=J
0375           IF(KCS.EQ.-1) JC=3-J
0376           IF(ICOL(KCC,1,JC).NE.0) K(N+1,J+3)=N+ICOL(KCC,1,JC)
0377           IF(ICOL(KCC,2,JC).NE.0) K(N+2,J+3)=N+ICOL(KCC,2,JC)
0378           IF(ICOL(KCC,3,JC).NE.0) K(N+3,J+3)=MSTU(5)*(N+ICOL(KCC,3,JC))
0379           IF(ICOL(KCC,4,JC).NE.0) K(N+4,J+3)=MSTU(5)*(N+ICOL(KCC,4,JC))
0380   270   CONTINUE
0381
0382 C...Store incoming and outgoing partons in their CM-frame.
0383         SHR=SQRT(VINT(21))*VINT(1)
0384         P(N+1,3)=0.5D0*SHR
0385         P(N+1,4)=0.5D0*SHR
0386         P(N+2,3)=-0.5D0*SHR
0387         P(N+2,4)=0.5D0*SHR
0388         P(N+3,5)=PYMASS(K(N+3,2))
0389         P(N+4,5)=PYMASS(K(N+4,2))
0390         IF(P(N+3,5)+P(N+4,5).GE.SHR) THEN
0391           IFAIL=1
0392           RETURN
0393         ENDIF
0394         P(N+3,4)=0.5D0*(SHR+(P(N+3,5)**2-P(N+4,5)**2)/SHR)
0395         P(N+3,3)=SQRT(MAX(0D0,P(N+3,4)**2-P(N+3,5)**2))
0396         P(N+4,4)=SHR-P(N+3,4)
0397         P(N+4,3)=-P(N+3,3)
0398
0399 C...Rotate outgoing partons using cos(theta)=(th-uh)/lam(sh,sqm3,sqm4)
0400         PHI=PARU(2)*PYR(0)
0401         CALL PYROBO(N+3,N+4,ACOS(VINT(23)),PHI,0D0,0D0,0D0)
0402
0403 C...Global statistics.
0404         MINT(351)=MINT(351)+1
0405         VINT(351)=VINT(351)+SQRT(P(N+3,1)**2+P(N+3,2)**2)
0406         IF (MINT(351).EQ.1) VINT(356)=SQRT(P(N+3,1)**2+P(N+3,2)**2)
0407
0408 C...Keep track of loose colour ends and information on scattering.
0409         MINT(31)=MINT(31)+1
0410         MINT(36)=MINT(31)
0411         PT2MI(MINT(36))=PT2
0412         IMISEP(MINT(31))=N+4
0413         DO 280 JS=1,2
0414           IMI(JS,MINT(31),1)=N+JS
0415           IMI(JS,MINT(31),2)=0
0416           XMI(JS,MINT(31))=VINT(40+JS)
0417           NMI(JS)=NMI(JS)+1
0418 C...Update cumulative counters
0419           VINT(142+JS)=VINT(142+JS)-VINT(40+JS)
0420           VINT(150+JS)=VINT(150+JS)+VINT(40+JS)
0421   280   CONTINUE
0422
0423 C...Add to list of final state partons
0424         IPART(NPART+1)=N+3
0425         IPART(NPART+2)=N+4
0426         PTPART(NPART+1)=SQRT(PT2)
0427         PTPART(NPART+2)=SQRT(PT2)
0428         NPART=NPART+2
0429
0430 C...Initialize ISR
0431         NISGEN(1,MINT(31))=0
0432         NISGEN(2,MINT(31))=0
0433
0434 C...Update ER
0435         N=N+4
0436         IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
0437           CALL PYERRM(11,'(PYMIGN:) no more memory left in PYJETS')
0438           MINT(51)=1
0439           RETURN
0440         ENDIF
0441
0442 C...Finally, assign colour tags to new partons
0443         DO 300 JS=1,2
0444           I1=IMI(JS,MINT(31),1)
0445           I2=IMI(3-JS,MINT(31),1)
0446           DO 290 JCS=4,5
0447             IF (K(I1,2).NE.21.AND.(9-2*JCS).NE.ISIGN(1,K(I1,2)))
0448      &           GOTO 290
0449             IF (K(I1,JCS)/MSTU(5)**2.NE.0) GOTO 290
0450             KCS=JCS
0451             CALL PYCTTR(I1,KCS,I2)
0452             IF(MINT(51).NE.0) RETURN
0453   290     CONTINUE
0454   300   CONTINUE
0455
0456 C----------------------------------------------------------------------
0457 C...MODE=2: Decide whether quarks in last scattering were valence,
0458 C...companion, or sea.
0459       ELSEIF (MODE.EQ.2) THEN
0460         JS=MINT(30)
0461         MI=MINT(36)
0462         PT2=PT2NOW
0463         KFSBM=ISIGN(1,MINT(10+JS))
0464         IFL=K(IMI(JS,MI,1),2)
0465         IMI(JS,MI,2)=0
0466         IF (IABS(IFL).GE.6) THEN
0467           IF (IABS(IFL).EQ.6) THEN
0468             CALL PYERRM(29,'(PYPTMI:) top in initial state!')
0469           ENDIF
0470           RETURN
0471         ENDIF
0472 C...Get PDFs at X(rescaled) and PT2 of the current initiator.
0473 C...(Do not include the parton itself in the X rescaling.)
0474         X=XMI(JS,MI)
0475         XRSC=X/(VINT(142+JS)+X)
0476 C...Note: XPSVC = x*pdf.
0477         MINT(30)=JS
0478         CALL PYPDFU(KFBEAM(JS),XRSC,PT2,XPQ)
0479         SEA=XPSVC(IFL,-1)
0480         VAL=XPSVC(IFL,0)
0481         CMP=0D0
0482         DO 310 IVC=1,NVC(JS,IFL)
0483           CMP=CMP+XPSVC(IFL,IVC)
0484   310   CONTINUE
0485
0486 C...Decide (Extra factor x cancels in the dvision).
0487   320   RVCS=PYR(0)*(SEA+VAL+CMP)
0488         IVNOW=1
0489   330   IF (RVCS.LE.VAL.AND.IVNOW.GE.1) THEN
0490 C...Safety check that valence present; pi0/gamma/K0S/K0L special cases.
0491           IVNOW=0
0492           IF(KFIVAL(JS,1).EQ.IFL) IVNOW=IVNOW+1
0493           IF(KFIVAL(JS,2).EQ.IFL) IVNOW=IVNOW+1
0494           IF(KFIVAL(JS,3).EQ.IFL) IVNOW=IVNOW+1
0495           IF(KFIVAL(JS,1).EQ.0) THEN
0496             IF(KFBEAM(JS).EQ.111.AND.IABS(IFL).LE.2) IVNOW=1
0497             IF(KFBEAM(JS).EQ.22.AND.IABS(IFL).LE.5) IVNOW=1
0498             IF((KFBEAM(JS).EQ.130.OR.KFBEAM(JS).EQ.310).AND.
0499      &           (IABS(IFL).EQ.1.OR.IABS(IFL).EQ.3)) IVNOW=1
0500           ELSE
0501 C...Count down valence remaining. Do not count current scattering.
0502             DO 340 I1=1,NMI(JS)
0503               IF (I1.EQ.MINT(36)) GOTO 340
0504               IF (K(IMI(JS,I1,1),2).EQ.IFL.AND.IMI(JS,I1,2).EQ.0)
0505      &             IVNOW=IVNOW-1
0506   340       CONTINUE
0507           ENDIF
0508           IF(IVNOW.EQ.0) GOTO 330
0509 C...Mark valence.
0510           IMI(JS,MI,2)=0
0511 C...Sets valence content of gamma, pi0, K0S, K0L if not done.
0512           IF(KFIVAL(JS,1).EQ.0) THEN
0513             IF(KFBEAM(JS).EQ.111.OR.KFBEAM(JS).EQ.22) THEN
0514               KFIVAL(JS,1)=IFL
0515               KFIVAL(JS,2)=-IFL
0516             ELSEIF(KFBEAM(JS).EQ.130.OR.KFBEAM(JS).EQ.310) THEN
0517               KFIVAL(JS,1)=IFL
0518               IF(IABS(IFL).EQ.1) KFIVAL(JS,2)=ISIGN(3,-IFL)
0519               IF(IABS(IFL).NE.1) KFIVAL(JS,2)=ISIGN(1,-IFL)
0520             ENDIF
0521           ENDIF
0522
0523         ELSEIF (RVCS.LE.VAL+SEA) THEN
0524 C...If sea, add opposite sign companion parton. Store X and I.
0525           NVC(JS,-IFL)=NVC(JS,-IFL)+1
0526           XASSOC(JS,-IFL,NVC(JS,-IFL))=XMI(JS,MI)
0527 C...Set pointer to companion
0528           IMI(JS,MI,2)=-NVC(JS,-IFL)
0529
0530         ELSE
0531 C...If companion, decide which one.
0532           IF (NVC(JS,IFL).EQ.0) THEN
0533             CMP=0D0
0534             CALL PYERRM(9,'(PYPTMI:) No cmp quark, but pdf != 0!')
0535             GOTO 320
0536           ENDIF
0537           CMPSUM=VAL+SEA
0538           ISEL=0
0539   350     ISEL=ISEL+1
0540           CMPSUM=CMPSUM+XPSVC(IFL,ISEL)
0541           IF (RVCS.GT.CMPSUM.AND.ISEL.LT.NVC(JS,IFL)) GOTO 350
0542 C...Find original sea (anti-)quark. Do not consider current scattering.
0543           IASSOC=0
0544           DO 360 I1=1,NMI(JS)
0545             IF (I1.EQ.MINT(36)) GOTO 360
0546             IF (K(IMI(JS,I1,1),2).NE.-IFL) GOTO 360
0547             IF (-IMI(JS,I1,2).EQ.ISEL) THEN
0548               IMI(JS,MI,2)=IMI(JS,I1,1)
0549               IMI(JS,I1,2)=IMI(JS,MI,1)
0550             ENDIF
0551   360     CONTINUE
0552 C...Mark companion "out-kicked".
0553           XASSOC(JS,IFL,ISEL)=-XASSOC(JS,IFL,ISEL)
0554         ENDIF
0555
0556       ENDIF
0557       RETURN
0558       END