Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:21:14 

0001  
0002 C*********************************************************************
0003  
0004 C...PYRAND
0005 C...Generates quantities characterizing the high-pT scattering at the
0006 C...parton level according to the matrix elements. Chooses incoming,
0007 C...reacting partons, their momentum fractions and one of the possible
0008 C...subprocesses.
0009  
0010       SUBROUTINE PYRAND
0011  
0012 C...Double precision and integer declarations.
0013       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0014       IMPLICIT INTEGER(I-N)
0015       INTEGER PYK,PYCHGE,PYCOMP
0016 C...Parameter statement to help give large particle numbers.
0017       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0018      &KEXCIT=4000000,KDIMEN=5000000)
0019  
0020 C...User process initialization and event commonblocks.
0021       INTEGER MAXPUP
0022       PARAMETER (MAXPUP=100)
0023       INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
0024       DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
0025       COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
0026      &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
0027      &LPRUP(MAXPUP)
0028       INTEGER MAXNUP
0029       PARAMETER (MAXNUP=500)
0030       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
0031       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
0032       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
0033      &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
0034      &VTIMUP(MAXNUP),SPINUP(MAXNUP)
0035       SAVE /HEPRUP/,/HEPEUP/
0036  
0037 C...Commonblocks.
0038       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0039       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0040       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0041       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0042       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0043       COMMON/PYINT1/MINT(400),VINT(400)
0044       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0045       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
0046       COMMON/PYINT4/MWID(500),WIDS(500,5)
0047       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
0048       COMMON/PYINT7/SIGT(0:6,0:6,0:5)
0049       COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
0050       SAVE /PYDAT1/,/PYDAT2/,/PYDAT3/,/PYSUBS/,/PYPARS/,/PYINT1/,
0051      &/PYINT2/,/PYINT3/,/PYINT4/,/PYINT5/,/PYINT7/,/PYMSSM/
0052 C...Local arrays.
0053       DIMENSION XPQ(-25:25),PMM(2),PDIF(4),BHAD(4),PMMN(2)
0054  
0055 C...Parameters and data used in elastic/diffractive treatment.
0056       DATA EPS/0.0808D0/, ALP/0.25D0/, CRES/2D0/, PMRC/1.062D0/,
0057      &SMP/0.880D0/, BHAD/2.3D0,1.4D0,1.4D0,0.23D0/
0058  
0059 C...Initial values, specifically for (first) semihard interaction.
0060       MINT(10)=0
0061       MINT(17)=0
0062       MINT(18)=0
0063       VINT(143)=1D0
0064       VINT(144)=1D0
0065       VINT(157)=0D0
0066       VINT(158)=0D0
0067       MFAIL=0
0068       IF(MSTP(171).EQ.1.AND.MSTP(172).EQ.2) MFAIL=1
0069       ISUB=0
0070       ISTSB=0
0071       LOOP=0
0072   100 LOOP=LOOP+1
0073       MINT(51)=0
0074       MINT(143)=1
0075       VINT(97)=1D0
0076  
0077 C...Start by assuming incoming photon is entering subprocess.
0078       IF(MINT(11).EQ.22) THEN
0079          MINT(15)=22
0080          VINT(307)=VINT(3)**2
0081       ENDIF
0082       IF(MINT(12).EQ.22) THEN
0083          MINT(16)=22
0084          VINT(308)=VINT(4)**2
0085       ENDIF
0086       MINT(103)=MINT(11)
0087       MINT(104)=MINT(12)
0088  
0089 C...Choice of process type - first event of pileup.
0090       INMULT=0
0091       IF(MINT(82).EQ.1.AND.ISUB.GE.91.AND.ISUB.LE.96) THEN
0092       ELSEIF(MINT(82).EQ.1) THEN
0093  
0094 C...For gamma-p or gamma-gamma first pick between alternatives.
0095         IGA=0
0096         IF(MINT(121).GT.1) CALL PYSAVE(4,IGA)
0097         MINT(122)=IGA
0098  
0099 C...For real gamma + gamma with different nature, flip at random.
0100         IF(MINT(11).EQ.22.AND.MINT(12).EQ.22.AND.MINT(123).GE.4.AND.
0101      &  MSTP(14).LE.10.AND.PYR(0).GT.0.5D0) THEN
0102           MINTSV=MINT(41)
0103           MINT(41)=MINT(42)
0104           MINT(42)=MINTSV
0105           MINTSV=MINT(45)
0106           MINT(45)=MINT(46)
0107           MINT(46)=MINTSV
0108           MINTSV=MINT(107)
0109           MINT(107)=MINT(108)
0110           MINT(108)=MINTSV
0111           IF(MINT(47).EQ.2.OR.MINT(47).EQ.3) MINT(47)=5-MINT(47)
0112         ENDIF
0113  
0114 C...Pick process type, possibly by user process machinery.
0115 C...(If the latter, also event will be picked here.)
0116         IF(MINT(111).GE.11.AND.IABS(IDWTUP).EQ.2.AND.LOOP.GE.2) THEN
0117           CALL UPEVNT
0118           CALL PYUPRE
0119         ELSEIF(MINT(111).GE.11.AND.IABS(IDWTUP).GE.3) THEN
0120           CALL UPEVNT
0121           CALL PYUPRE
0122           ISUB=0
0123   110     ISUB=ISUB+1
0124           IF((ISET(ISUB).NE.11.OR.KFPR(ISUB,2).NE.IDPRUP).AND.
0125      &    ISUB.LT.500) GOTO 110
0126         ELSE
0127           RSUB=XSEC(0,1)*PYR(0)
0128           DO 120 I=1,500
0129             IF(MSUB(I).NE.1.OR.I.EQ.96) GOTO 120
0130             ISUB=I
0131             RSUB=RSUB-XSEC(I,1)
0132             IF(RSUB.LE.0D0) GOTO 130
0133   120     CONTINUE
0134   130     IF(ISUB.EQ.95) ISUB=96
0135           IF(ISUB.EQ.96) INMULT=1
0136           IF(ISET(ISUB).EQ.11) THEN
0137             IDPRUP=KFPR(ISUB,2)
0138             CALL UPEVNT
0139             CALL PYUPRE
0140           ENDIF
0141         ENDIF
0142  
0143 C...Choice of inclusive process type - pileup events.
0144       ELSEIF(MINT(82).GE.2.AND.ISUB.EQ.0) THEN
0145         RSUB=VINT(131)*PYR(0)
0146         ISUB=96
0147         IF(RSUB.GT.SIGT(0,0,5)) ISUB=94
0148         IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)) ISUB=93
0149         IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)+SIGT(0,0,3)) ISUB=92
0150         IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)+SIGT(0,0,3)+SIGT(0,0,2))
0151      &  ISUB=91
0152         IF(ISUB.EQ.96) INMULT=1
0153       ENDIF
0154  
0155 C...Choice of photon energy and flux factor inside lepton.
0156       IF(MINT(141).NE.0.OR.MINT(142).NE.0) THEN
0157         IF (MSTP(199).EQ.1) THEN
0158           CALL PYGAGA(5,WTGAGA)
0159         ELSE
0160           CALL PYGAGA(3,WTGAGA)
0161         ENDIF
0162         IF(ISUB.GE.131.AND.ISUB.LE.140) THEN
0163           CKIN(3)=MAX(VINT(285),VINT(154))
0164           CKIN(1)=2D0*CKIN(3)
0165         ENDIF
0166 C...When necessary set direct/resolved photon by hand.
0167       ELSEIF(MINT(15).EQ.22.OR.MINT(16).EQ.22) THEN
0168         IF(MINT(15).EQ.22.AND.MINT(41).EQ.2) MINT(15)=0
0169         IF(MINT(16).EQ.22.AND.MINT(42).EQ.2) MINT(16)=0
0170       ENDIF
0171  
0172 C...Restrict direct*resolved processes to pTmin >= Q,
0173 C...to avoid doublecounting  with DIS.
0174       IF(MSTP(18).EQ.3.AND.ISUB.GE.131.AND.ISUB.LE.136) THEN
0175         IF(MINT(15).EQ.22) THEN
0176           CKIN(3)=MAX(VINT(285),VINT(154),ABS(VINT(3)))
0177         ELSE
0178           CKIN(3)=MAX(VINT(285),VINT(154),ABS(VINT(4)))
0179         ENDIF
0180         CKIN(1)=2D0*CKIN(3)
0181       ENDIF
0182  
0183 C...Set up for multiple interactions (may include impact parameter).
0184       IF(INMULT.EQ.1) THEN
0185         IF(MINT(35).LE.1) CALL PYMULT(2)
0186         IF(MINT(35).GE.2) CALL PYMIGN(2)
0187       ENDIF
0188  
0189 C...Loopback point for minimum bias in photon physics.
0190       LOOP2=0
0191   140 LOOP2=LOOP2+1
0192       IF(MINT(82).EQ.1) NGEN(0,1)=NGEN(0,1)+MINT(143)
0193       IF(MINT(82).EQ.1) NGEN(ISUB,1)=NGEN(ISUB,1)+MINT(143)
0194       IF(ISUB.EQ.96.AND.LOOP2.EQ.1.AND.MINT(82).EQ.1)
0195      &NGEN(97,1)=NGEN(97,1)+MINT(143)
0196       MINT(1)=ISUB
0197       ISTSB=ISET(ISUB)
0198  
0199 C...Random choice of flavour for some SUSY processes.
0200       IF(ISUB.GE.201.AND.ISUB.LE.301) THEN
0201 C...~e_L ~nu_e or ~mu_L ~nu_mu.
0202         IF(ISUB.EQ.210) THEN
0203           KFPR(ISUB,1)=KSUSY1+11+2*INT(0.5D0+PYR(0))
0204           KFPR(ISUB,2)=KFPR(ISUB,1)+1
0205 C...~nu_e ~nu_e(bar) or ~nu_mu ~nu_mu(bar).
0206         ELSEIF(ISUB.EQ.213) THEN
0207           KFPR(ISUB,1)=KSUSY1+12+2*INT(0.5D0+PYR(0))
0208           KFPR(ISUB,2)=KFPR(ISUB,1)
0209 C...~q ~chi/~g; ~q = ~d, ~u, ~s, ~c or ~b.
0210         ELSEIF(ISUB.GE.246.AND.ISUB.LE.259.AND.ISUB.NE.255.AND.
0211      &  ISUB.NE.257) THEN
0212           IF(ISUB.GE.258) THEN
0213             RKF=4D0
0214           ELSE
0215             RKF=5D0
0216           ENDIF
0217           IF(MOD(ISUB,2).EQ.0) THEN
0218             KFPR(ISUB,1)=KSUSY1+1+INT(RKF*PYR(0))
0219           ELSE
0220             KFPR(ISUB,1)=KSUSY2+1+INT(RKF*PYR(0))
0221           ENDIF
0222 C...~q1 ~q2; ~q = ~d, ~u, ~s, or ~c.
0223         ELSEIF(ISUB.GE.271.AND.ISUB.LE.276) THEN
0224           IF(ISUB.EQ.271.OR.ISUB.EQ.274) THEN
0225             KSU1=KSUSY1
0226             KSU2=KSUSY1
0227           ELSEIF(ISUB.EQ.272.OR.ISUB.EQ.275) THEN
0228             KSU1=KSUSY2
0229             KSU2=KSUSY2
0230           ELSEIF(PYR(0).LT.0.5D0) THEN
0231             KSU1=KSUSY1
0232             KSU2=KSUSY2
0233           ELSE
0234             KSU1=KSUSY2
0235             KSU2=KSUSY1
0236           ENDIF
0237           KFPR(ISUB,1)=KSU1+1+INT(4D0*PYR(0))
0238           KFPR(ISUB,2)=KSU2+1+INT(4D0*PYR(0))
0239 C...~q ~q(bar);  ~q = ~d, ~u, ~s, or ~c.
0240         ELSEIF(ISUB.EQ.277.OR.ISUB.EQ.279) THEN
0241           KFPR(ISUB,1)=KSUSY1+1+INT(4D0*PYR(0))
0242           KFPR(ISUB,2)=KFPR(ISUB,1)
0243         ELSEIF(ISUB.EQ.278.OR.ISUB.EQ.280) THEN
0244           KFPR(ISUB,1)=KSUSY2+1+INT(4D0*PYR(0))
0245           KFPR(ISUB,2)=KFPR(ISUB,1)
0246 C...~q1 ~q2; ~q = ~d, ~u, ~s, or ~c.
0247         ELSEIF(ISUB.GE.281.AND.ISUB.LE.286) THEN
0248           IF(ISUB.EQ.281.OR.ISUB.EQ.284) THEN
0249             KSU1=KSUSY1
0250             KSU2=KSUSY1
0251           ELSEIF(ISUB.EQ.282.OR.ISUB.EQ.285) THEN
0252             KSU1=KSUSY2
0253             KSU2=KSUSY2
0254           ELSEIF(PYR(0).LT.0.5D0) THEN
0255             KSU1=KSUSY1
0256             KSU2=KSUSY2
0257           ELSE
0258             KSU1=KSUSY2
0259             KSU2=KSUSY1
0260           ENDIF
0261           IF(ISUB.EQ.281.OR.ISUB.LE.283) THEN
0262             RKF=5D0
0263           ELSE
0264             RKF=4D0
0265           ENDIF
0266           KFPR(ISUB,2)=KSU2+1+INT(RKF*PYR(0))
0267         ENDIF
0268       ENDIF
0269  
0270 C...Find resonances (explicit or implicit in cross-section).
0271       MINT(72)=0
0272       KFR1=0
0273       IF(ISTSB.EQ.1.OR.ISTSB.EQ.3.OR.ISTSB.EQ.5) THEN
0274         KFR1=KFPR(ISUB,1)
0275       ELSEIF(ISUB.EQ.24.OR.ISUB.EQ.25.OR.ISUB.EQ.110.OR.ISUB.EQ.165.OR.
0276      &  ISUB.EQ.171.OR.ISUB.EQ.176) THEN
0277         KFR1=23
0278       ELSEIF(ISUB.EQ.23.OR.ISUB.EQ.26.OR.ISUB.EQ.166.OR.ISUB.EQ.172.OR.
0279      &  ISUB.EQ.177) THEN
0280         KFR1=24
0281       ELSEIF(ISUB.GE.71.AND.ISUB.LE.77) THEN
0282         KFR1=25
0283         IF(MSTP(46).EQ.5) THEN
0284           KFR1=89
0285           PMAS(89,1)=PARP(45)
0286           PMAS(89,2)=PARP(45)**3/(96D0*PARU(1)*PARP(47)**2)
0287         ENDIF
0288       ELSEIF(ISUB.EQ.194) THEN
0289         KFR1=KTECHN+113
0290       ELSEIF(ISUB.EQ.195) THEN
0291         KFR1=KTECHN+213
0292       ELSEIF(ISUB.GE.361.AND.ISUB.LE.368) THEN
0293         KFR1=KTECHN+113
0294       ELSEIF(ISUB.GE.370.AND.ISUB.LE.377) THEN
0295         KFR1=KTECHN+213
0296       ENDIF
0297       CKMX=CKIN(2)
0298       IF(CKMX.LE.0D0) CKMX=VINT(1)
0299       KCR1=PYCOMP(KFR1)
0300       IF(KFR1.NE.0) THEN
0301         IF(CKIN(1).GT.PMAS(KCR1,1)+20D0*PMAS(KCR1,2).OR.
0302      &  CKMX.LT.PMAS(KCR1,1)-20D0*PMAS(KCR1,2)) KFR1=0
0303       ENDIF
0304       IF(KFR1.NE.0) THEN
0305         TAUR1=PMAS(KCR1,1)**2/VINT(2)
0306         IF(KFR1.EQ.KTECHN+113) THEN
0307           CALL PYTECM(S1,S2)
0308           TAUR1=S1/VINT(2)
0309         ENDIF
0310         GAMR1=PMAS(KCR1,1)*PMAS(KCR1,2)/VINT(2)
0311         MINT(72)=1
0312         MINT(73)=KFR1
0313         VINT(73)=TAUR1
0314         VINT(74)=GAMR1
0315       ENDIF
0316       IF(ISUB.EQ.141.OR.ISUB.EQ.194.OR.(ISUB.GE.364.AND.ISUB.LE.368))
0317      $THEN
0318         KFR2=23
0319         IF(ISUB.EQ.194) THEN
0320           KFR2=KTECHN+223
0321         ELSEIF(ISUB.GE.364.AND.ISUB.LE.368) THEN
0322           KFR2=KTECHN+223
0323         ENDIF
0324         KCR2=PYCOMP(KFR2)
0325         TAUR2=PMAS(KCR2,1)**2/VINT(2)
0326         IF(KFR2.EQ.KTECHN+223) THEN
0327           CALL PYTECM(S1,S2)
0328           TAUR2=S2/VINT(2)
0329         ENDIF
0330         GAMR2=PMAS(KCR2,1)*PMAS(KCR2,2)/VINT(2)
0331         IF(CKIN(1).GT.PMAS(KCR2,1)+20D0*PMAS(KCR2,2).OR.
0332      &  CKMX.LT.PMAS(KCR2,1)-20D0*PMAS(KCR2,2)) KFR2=0
0333         IF(KFR2.NE.0.AND.KFR1.NE.0) THEN
0334           MINT(72)=2
0335           MINT(74)=KFR2
0336           VINT(75)=TAUR2
0337           VINT(76)=GAMR2
0338         ELSEIF(KFR2.NE.0) THEN
0339           KFR1=KFR2
0340           TAUR1=TAUR2
0341           GAMR1=GAMR2
0342           MINT(72)=1
0343           MINT(73)=KFR1
0344           VINT(73)=TAUR1
0345           VINT(74)=GAMR1
0346         ENDIF
0347       ENDIF
0348  
0349 C...Find product masses and minimum pT of process,
0350 C...optionally with broadening according to a truncated Breit-Wigner.
0351       VINT(63)=0D0
0352       VINT(64)=0D0
0353       MINT(71)=0
0354       VINT(71)=CKIN(3)
0355       IF(MINT(82).GE.2) VINT(71)=0D0
0356       VINT(80)=1D0
0357       IF(ISTSB.EQ.2.OR.ISTSB.EQ.4) THEN
0358         NBW=0
0359         DO 160 I=1,2
0360           PMMN(I)=0D0
0361           IF(KFPR(ISUB,I).EQ.0) THEN
0362           ELSEIF(MSTP(42).LE.0.OR.PMAS(PYCOMP(KFPR(ISUB,I)),2).LT.
0363      &      PARP(41)) THEN
0364             VINT(62+I)=PMAS(PYCOMP(KFPR(ISUB,I)),1)**2
0365           ELSE
0366             NBW=NBW+1
0367 C...This prevents SUSY/t particles from becoming too light.
0368             KFLW=KFPR(ISUB,I)
0369             IF(KFLW/KSUSY1.EQ.1.OR.KFLW/KSUSY1.EQ.2) THEN
0370               KCW=PYCOMP(KFLW)
0371               PMMN(I)=PMAS(KCW,1)
0372               DO 150 IDC=MDCY(KCW,2),MDCY(KCW,2)+MDCY(KCW,3)-1
0373                 IF(MDME(IDC,1).GT.0.AND.BRAT(IDC).GT.1E-4) THEN
0374                   PMSUM=PMAS(PYCOMP(KFDP(IDC,1)),1)+
0375      &            PMAS(PYCOMP(KFDP(IDC,2)),1)
0376                   IF(KFDP(IDC,3).NE.0) PMSUM=PMSUM+
0377      &            PMAS(PYCOMP(KFDP(IDC,3)),1)
0378                   PMMN(I)=MIN(PMMN(I),PMSUM)
0379                 ENDIF
0380   150         CONTINUE
0381             ELSEIF(KFLW.EQ.6) THEN
0382               PMMN(I)=PMAS(24,1)+PMAS(5,1)
0383             ENDIF
0384           ENDIF
0385   160   CONTINUE
0386         IF(NBW.GE.1) THEN
0387           CKIN41=CKIN(41)
0388           CKIN43=CKIN(43)
0389           CKIN(41)=MAX(PMMN(1),CKIN(41))
0390           CKIN(43)=MAX(PMMN(2),CKIN(43))
0391           CALL PYOFSH(4,0,KFPR(ISUB,1),KFPR(ISUB,2),0D0,PQM3,PQM4)
0392           CKIN(41)=CKIN41
0393           CKIN(43)=CKIN43
0394           IF(MINT(51).EQ.1) THEN
0395             IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
0396             IF(MFAIL.EQ.1) THEN
0397               MSTI(61)=1
0398               RETURN
0399             ENDIF
0400             GOTO 100
0401           ENDIF
0402           VINT(63)=PQM3**2
0403           VINT(64)=PQM4**2
0404         ENDIF
0405         IF(MIN(VINT(63),VINT(64)).LT.CKIN(6)**2) MINT(71)=1
0406         IF(MINT(71).EQ.1) VINT(71)=MAX(CKIN(3),CKIN(5))
0407       ENDIF
0408  
0409 C...Prepare for additional variable choices in 2 -> 3.
0410       IF(ISTSB.EQ.5) THEN
0411         VINT(201)=0D0
0412         IF(KFPR(ISUB,2).GT.0) VINT(201)=PMAS(PYCOMP(KFPR(ISUB,2)),1)
0413         VINT(206)=VINT(201)
0414         IF(ISUB.EQ.401.OR.ISUB.EQ.402) VINT(206)=PMAS(5,1)
0415         VINT(204)=PMAS(23,1)
0416         IF(ISUB.EQ.124.OR.ISUB.EQ.174.OR.ISUB.EQ.179.OR.ISUB.EQ.351)
0417      &   VINT(204)=PMAS(24,1) 
0418         IF(ISUB.EQ.352) VINT(204)=PMAS(PYCOMP(9900024),1)
0419         IF(ISUB.EQ.121.OR.ISUB.EQ.122.OR.ISUB.EQ.181.OR.ISUB.EQ.182.OR.
0420      &    ISUB.EQ.186.OR.ISUB.EQ.187.OR.ISUB.EQ.401.OR.ISUB.EQ.402)
0421      &         VINT(204)=VINT(201)
0422         VINT(209)=VINT(204)
0423           IF(ISUB.EQ.401.OR.ISUB.EQ.402) VINT(209)=VINT(206)
0424       ENDIF
0425  
0426 C...Select incoming VDM particle (rho/omega/phi/J/psi).
0427       IF(ISTSB.NE.0.AND.(MINT(101).GE.2.OR.MINT(102).GE.2).AND.
0428      &(MINT(123).EQ.2.OR.MINT(123).EQ.3.OR.MINT(123).EQ.7)) THEN
0429         VRN=PYR(0)*SIGT(0,0,5)
0430         IF(MINT(101).LE.1) THEN
0431           I1MN=0
0432           I1MX=0
0433         ELSE
0434           I1MN=1
0435           I1MX=MINT(101)
0436         ENDIF
0437         IF(MINT(102).LE.1) THEN
0438           I2MN=0
0439           I2MX=0
0440         ELSE
0441           I2MN=1
0442           I2MX=MINT(102)
0443         ENDIF
0444         DO 180 I1=I1MN,I1MX
0445           KFV1=110*I1+3
0446           DO 170 I2=I2MN,I2MX
0447             KFV2=110*I2+3
0448             VRN=VRN-SIGT(I1,I2,5)
0449             IF(VRN.LE.0D0) GOTO 190
0450   170     CONTINUE
0451   180   CONTINUE
0452   190   IF(MINT(101).GE.2) MINT(103)=KFV1
0453         IF(MINT(102).GE.2) MINT(104)=KFV2
0454       ENDIF
0455  
0456       IF(ISTSB.EQ.0) THEN
0457 C...Elastic scattering or single or double diffractive scattering.
0458  
0459 C...Select incoming particle (rho/omega/phi/J/psi for VDM) and mass.
0460         MINT(103)=MINT(11)
0461         MINT(104)=MINT(12)
0462         PMM(1)=VINT(3)
0463         PMM(2)=VINT(4)
0464         IF(MINT(101).GE.2.OR.MINT(102).GE.2) THEN
0465           JJ=ISUB-90
0466           VRN=PYR(0)*SIGT(0,0,JJ)
0467           IF(MINT(101).LE.1) THEN
0468             I1MN=0
0469             I1MX=0
0470           ELSE
0471             I1MN=1
0472             I1MX=MINT(101)
0473           ENDIF
0474           IF(MINT(102).LE.1) THEN
0475             I2MN=0
0476             I2MX=0
0477           ELSE
0478             I2MN=1
0479             I2MX=MINT(102)
0480           ENDIF
0481           DO 210 I1=I1MN,I1MX
0482             KFV1=110*I1+3
0483             DO 200 I2=I2MN,I2MX
0484               KFV2=110*I2+3
0485               VRN=VRN-SIGT(I1,I2,JJ)
0486               IF(VRN.LE.0D0) GOTO 220
0487   200       CONTINUE
0488   210     CONTINUE
0489   220     IF(MINT(101).GE.2) THEN
0490             MINT(103)=KFV1
0491             PMM(1)=PYMASS(KFV1)
0492           ENDIF
0493           IF(MINT(102).GE.2) THEN
0494             MINT(104)=KFV2
0495             PMM(2)=PYMASS(KFV2)
0496           ENDIF
0497         ENDIF
0498         VINT(67)=PMM(1)
0499         VINT(68)=PMM(2)
0500  
0501 C...Select mass for GVMD states (rejecting previous assignment).
0502         Q0S=4D0*PARP(15)**2
0503         Q1S=4D0*VINT(154)**2
0504         LOOP3=0
0505   230   LOOP3=LOOP3+1
0506         DO 240 JT=1,2
0507           IF(MINT(106+JT).EQ.3) THEN
0508             PS=VINT(2+JT)**2
0509             PMM(JT)=(Q0S+PS)*(Q1S+PS)/
0510      &      (Q0S+PYR(0)*(Q1S-Q0S)+PS)-PS
0511             IF(MINT(102+JT).GE.333) PMM(JT)=PMM(JT)-
0512      &      PMAS(PYCOMP(113),1)+PMAS(PYCOMP(MINT(102+JT)),1)
0513           ENDIF
0514   240   CONTINUE
0515         IF(PMM(1)+PMM(2)+PARP(104).GE.VINT(1)) THEN
0516           IF(LOOP3.LT.100.AND.(MINT(107).EQ.3.OR.MINT(108).EQ.3))
0517      &    GOTO 230
0518           GOTO 100
0519         ENDIF
0520  
0521 C...Side/sides of diffractive system.
0522         MINT(17)=0
0523         MINT(18)=0
0524         IF(ISUB.EQ.92.OR.ISUB.EQ.94) MINT(17)=1
0525         IF(ISUB.EQ.93.OR.ISUB.EQ.94) MINT(18)=1
0526  
0527 C...Find masses of particles and minimal masses of diffractive states.
0528         DO 250 JT=1,2
0529           PDIF(JT)=PMM(JT)
0530           VINT(68+JT)=PDIF(JT)
0531           IF(MINT(16+JT).EQ.1) PDIF(JT)=PDIF(JT)+PARP(102)
0532   250   CONTINUE
0533         SH=VINT(2)
0534         if (MINT(11).eq.22) then
0535           SQM1 = -1.*(VINT(3)**2)
0536         else
0537           SQM1=PMM(1)**2
0538         endif
0539         SQM2=PMM(2)**2
0540         SQM3=PDIF(1)**2
0541         SQM4=PDIF(2)**2
0542         SMRES1=(PMM(1)+PMRC)**2
0543         SMRES2=(PMM(2)+PMRC)**2
0544  
0545 C...Find elastic slope and lower limit diffractive slope.
0546         IHA=MAX(2,IABS(MINT(103))/110)
0547         IF(IHA.GE.5) IHA=1
0548         IHB=MAX(2,IABS(MINT(104))/110)
0549         IF(IHB.GE.5) IHB=1
0550         IF(ISUB.EQ.91) THEN
0551           IF(IHA.eq.2) then
0552             PMVIRT=PMAS(PYCOMP(113),1)
0553             BMN=5.84/(1+(PARP(165))*(VINT(307)/(PMVIRT**2))**PARP(166))
0554      &          +4.5
0555           ELSEIF(IHA.eq.3) then
0556             PMVIRT=PMAS(PYCOMP(333),1)
0557             BMN=5.84/(1+(PARP(165))*(VINT(307)/(PMVIRT**2))**PARP(166))
0558      &          +4.5
0559           elseif(IHA.eq.4) then
0560             PMVIRT=PMAS(PYCOMP(443),1)
0561             BMN=4.5
0562           else
0563             BMN=2D0*BHAD(IHA)+2D0*BHAD(IHB)+4D0*SH**EPS-4.2D0
0564           endif
0565         ELSEIF(ISUB.EQ.92) THEN
0566           BMN=MAX(2D0,2D0*BHAD(IHB))
0567         ELSEIF(ISUB.EQ.93) THEN
0568           BMN=MAX(2D0,2D0*BHAD(IHA))
0569         ELSEIF(ISUB.EQ.94) THEN
0570           BMN=2D0*ALP*4D0
0571         ENDIF
0572  
0573 C...Determine maximum possible t range and coefficient of generation.
0574         SQLA12=(SH-SQM1-SQM2)**2-4D0*SQM1*SQM2
0575         SQLA34=(SH-SQM3-SQM4)**2-4D0*SQM3*SQM4
0576         THA=SH-(SQM1+SQM2+SQM3+SQM4)+(SQM1-SQM2)*(SQM3-SQM4)/SH
0577         THB=SQRT(MAX(0D0,SQLA12))*SQRT(MAX(0D0,SQLA34))/SH
0578         THC=(SQM3-SQM1)*(SQM4-SQM2)+(SQM1+SQM4-SQM2-SQM3)*
0579      &  (SQM1*SQM4-SQM2*SQM3)/SH
0580         THL=-0.5D0*(THA+THB)
0581         THU=THC/THL
0582         THRND=EXP(MAX(-50D0,BMN*(THL-THU)))-1D0
0583  
0584 C...Select diffractive mass/masses according to dm^2/m^2.
0585         LOOP3=0
0586   260   LOOP3=LOOP3+1
0587         DO 270 JT=1,2
0588           IF(MINT(16+JT).EQ.0) THEN
0589             PDIF(2+JT)=PDIF(JT)
0590           ELSE
0591             PMMIN=PDIF(JT)
0592             PMMAX=MAX(VINT(2+JT),VINT(1)-PDIF(3-JT))
0593             PDIF(2+JT)=PMMIN*(PMMAX/PMMIN)**PYR(0)
0594           ENDIF
0595   270   CONTINUE
0596         SQM3=PDIF(3)**2
0597         SQM4=PDIF(4)**2
0598  
0599 C..Additional mass factors, including resonance enhancement.
0600         IF(PDIF(3)+PDIF(4).GE.VINT(1)) THEN
0601           IF(LOOP3.LT.100) GOTO 260
0602           GOTO 100
0603         ENDIF
0604         IF(ISUB.EQ.92) THEN
0605           FSD=(1D0-SQM3/SH)*(1D0+CRES*SMRES1/(SMRES1+SQM3))
0606           IF(FSD.LT.PYR(0)*(1D0+CRES)) GOTO 260
0607         ELSEIF(ISUB.EQ.93) THEN
0608           FSD=(1D0-SQM4/SH)*(1D0+CRES*SMRES2/(SMRES2+SQM4))
0609           IF(FSD.LT.PYR(0)*(1D0+CRES)) GOTO 260
0610         ELSEIF(ISUB.EQ.94) THEN
0611           FDD=(1D0-(PDIF(3)+PDIF(4))**2/SH)*(SH*SMP/
0612      &    (SH*SMP+SQM3*SQM4))*(1D0+CRES*SMRES1/(SMRES1+SQM3))*
0613      &    (1D0+CRES*SMRES2/(SMRES2+SQM4))
0614           IF(FDD.LT.PYR(0)*(1D0+CRES)**2) GOTO 260
0615         ENDIF
0616  
0617 C...Select t according to exp(Bmn*t) and correct to right slope.
0618         TH=THU+LOG(1D0+THRND*PYR(0))/BMN
0619 C        if (PMVIRT.eq.PMAS(PYCOMP(443),1)) then
0620 C        write(*,*) BMN, PMVIRT, VINT(307), MINT(101), PMM(1), PMM(2), 
0621 C     &             SQM1, SQM2, "THA: ", THA, "THB: ", THB, "THC: ", THC,
0622 C     &             "THL: ", THL, "THU: ", THU, "THRND: ", THRND, 
0623 C     &              "TH: ", TH    
0624 C        write(*,*) "======================="
0625 C        endif
0626         IF(ISUB.GE.92) THEN
0627           IF(ISUB.EQ.92) THEN
0628             BADD=2D0*ALP*LOG(SH/SQM3)
0629             IF(BHAD(IHB).LT.1D0) BADD=MAX(0D0,BADD+2D0*BHAD(IHB)-2D0)
0630           ELSEIF(ISUB.EQ.93) THEN
0631             BADD=2D0*ALP*LOG(SH/SQM4)
0632             IF(BHAD(IHA).LT.1D0) BADD=MAX(0D0,BADD+2D0*BHAD(IHA)-2D0)
0633           ELSEIF(ISUB.EQ.94) THEN
0634             BADD=2D0*ALP*(LOG(EXP(4D0)+SH/(ALP*SQM3*SQM4))-4D0)
0635           ENDIF
0636           IF(EXP(MAX(-50D0,BADD*(TH-THU))).LT.PYR(0)) GOTO 260
0637         ENDIF
0638  
0639 C...Check whether m^2 and t choices are consistent.
0640         SQLA34=(SH-SQM3-SQM4)**2-4D0*SQM3*SQM4
0641         THA=SH-(SQM1+SQM2+SQM3+SQM4)+(SQM1-SQM2)*(SQM3-SQM4)/SH
0642         THB=SQRT(MAX(0D0,SQLA12))*SQRT(MAX(0D0,SQLA34))/SH
0643         IF(THB.LE.1D-8) GOTO 260
0644         THC=(SQM3-SQM1)*(SQM4-SQM2)+(SQM1+SQM4-SQM2-SQM3)*
0645      &  (SQM1*SQM4-SQM2*SQM3)/SH
0646         THLM=-0.5D0*(THA+THB)
0647         THUM=THC/THLM
0648         IF(TH.LT.THLM.OR.TH.GT.THUM) GOTO 260
0649  
0650 C...Information to output.
0651         VINT(21)=1D0
0652         VINT(22)=0D0
0653         VINT(23)=MIN(1D0,MAX(-1D0,(THA+2D0*TH)/THB))
0654         VINT(45)=TH
0655         VINT(59)=2D0*SQRT(MAX(0D0,-(THC+THA*TH+TH**2)))/THB
0656         VINT(63)=PDIF(3)**2
0657         VINT(64)=PDIF(4)**2
0658         VINT(283)=PMM(1)**2/4D0
0659         VINT(284)=PMM(2)**2/4D0
0660  
0661 C...Note: in the following, by In is meant the integral over the
0662 C...quantity multiplying coefficient cn.
0663 C...Choose tau according to h1(tau)/tau, where
0664 C...h1(tau) = c1 + I1/I2*c2*1/tau + I1/I3*c3*1/(tau+tau_R) +
0665 C...I1/I4*c4*tau/((s*tau-m^2)^2+(m*Gamma)^2) +
0666 C...I1/I5*c5*1/(tau+tau_R') +
0667 C...I1/I6*c6*tau/((s*tau-m'^2)^2+(m'*Gamma')^2) +
0668 C...I1/I7*c7*tau/(1.-tau), and
0669 C...c1 + c2 + c3 + c4 + c5 + c6 + c7 = 1.
0670       ELSEIF(ISTSB.GE.1.AND.ISTSB.LE.5) THEN
0671         CALL PYKLIM(1)
0672         IF(MINT(51).NE.0) THEN
0673           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
0674           IF(MFAIL.EQ.1) THEN
0675             MSTI(61)=1
0676             RETURN
0677           ENDIF
0678           GOTO 100
0679         ENDIF
0680         RTAU=PYR(0)
0681         MTAU=1
0682         IF(RTAU.GT.COEF(ISUB,1)) MTAU=2
0683         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)) MTAU=3
0684         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)) MTAU=4
0685         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4))
0686      &  MTAU=5
0687         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+
0688      &  COEF(ISUB,5)) MTAU=6
0689         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+
0690      &  COEF(ISUB,5)+COEF(ISUB,6)) MTAU=7
0691         CALL PYKMAP(1,MTAU,PYR(0))
0692  
0693 C...2 -> 3, 4 processes:
0694 C...Choose tau' according to h4(tau,tau')/tau', where
0695 C...h4(tau,tau') = c1 + I1/I2*c2*(1 - tau/tau')^3/tau' +
0696 C...I1/I3*c3*1/(1 - tau'), and c1 + c2 + c3 = 1.
0697         IF(ISTSB.GE.3.AND.ISTSB.LE.5) THEN
0698           CALL PYKLIM(4)
0699           IF(MINT(51).NE.0) THEN
0700             IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
0701             IF(MFAIL.EQ.1) THEN
0702               MSTI(61)=1
0703               RETURN
0704             ENDIF
0705             GOTO 100
0706           ENDIF
0707           RTAUP=PYR(0)
0708           MTAUP=1
0709           IF(RTAUP.GT.COEF(ISUB,18)) MTAUP=2
0710           IF(RTAUP.GT.COEF(ISUB,18)+COEF(ISUB,19)) MTAUP=3
0711           CALL PYKMAP(4,MTAUP,PYR(0))
0712         ENDIF
0713  
0714 C...Choose y* according to h2(y*), where
0715 C...h2(y*) = I0/I1*c1*(y*-y*min) + I0/I2*c2*(y*max-y*) +
0716 C...I0/I3*c3*1/cosh(y*) + I0/I4*c4*1/(1-exp(y*-y*max)) +
0717 C...I0/I5*c5*1/(1-exp(-y*-y*min)), I0 = y*max-y*min,
0718 C...and c1 + c2 + c3 + c4 + c5 = 1.
0719         CALL PYKLIM(2)
0720         IF(MINT(51).NE.0) THEN
0721           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
0722           IF(MFAIL.EQ.1) THEN
0723             MSTI(61)=1
0724             RETURN
0725           ENDIF
0726           GOTO 100
0727         ENDIF
0728         RYST=PYR(0)
0729         MYST=1
0730         IF(RYST.GT.COEF(ISUB,8)) MYST=2
0731         IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
0732         IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)+COEF(ISUB,10)) MYST=4
0733         IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)+COEF(ISUB,10)+
0734      &  COEF(ISUB,11)) MYST=5
0735         CALL PYKMAP(2,MYST,PYR(0))
0736  
0737 C...2 -> 2 processes:
0738 C...Choose cos(theta-hat) (cth) according to h3(cth), where
0739 C...h3(cth) = c0 + I0/I1*c1*1/(A - cth) + I0/I2*c2*1/(A + cth) +
0740 C...I0/I3*c3*1/(A - cth)^2 + I0/I4*c4*1/(A + cth)^2,
0741 C...A = 1 + 2*(m3*m4/sh)^2 (= 1 for massless products),
0742 C...and c0 + c1 + c2 + c3 + c4 = 1.
0743         CALL PYKLIM(3)
0744         IF(MINT(51).NE.0) THEN
0745           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
0746           IF(MFAIL.EQ.1) THEN
0747             MSTI(61)=1
0748             RETURN
0749           ENDIF
0750           GOTO 100
0751         ENDIF
0752         IF(ISTSB.EQ.2.OR.ISTSB.EQ.4) THEN
0753           RCTH=PYR(0)
0754           MCTH=1
0755           IF(RCTH.GT.COEF(ISUB,13)) MCTH=2
0756           IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)) MCTH=3
0757           IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)+COEF(ISUB,15)) MCTH=4
0758           IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)+COEF(ISUB,15)+
0759      &    COEF(ISUB,16)) MCTH=5
0760           CALL PYKMAP(3,MCTH,PYR(0))
0761         ENDIF
0762  
0763 C...2 -> 3 : select pT1, phi1, pT2, phi2, y3 for 3 outgoing.
0764         IF(ISTSB.EQ.5) THEN
0765           CALL PYKMAP(5,0,0D0)
0766           IF(MINT(51).NE.0) THEN
0767             IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
0768             IF(MFAIL.EQ.1) THEN
0769               MSTI(61)=1
0770               RETURN
0771             ENDIF
0772             GOTO 100
0773           ENDIF
0774         ENDIF
0775  
0776 C...DIS as f + gamma* -> f process: set dummy values.
0777       ELSEIF(ISTSB.EQ.8) THEN
0778         VINT(21)=0.9D0
0779         VINT(22)=0D0
0780         VINT(23)=0D0
0781         VINT(47)=0D0
0782         VINT(48)=0D0
0783  
0784 C...Low-pT or multiple interactions (first semihard interaction).
0785       ELSEIF(ISTSB.EQ.9) THEN
0786         IF(MINT(35).LE.1) CALL PYMULT(3)
0787         IF(MINT(35).GE.2) CALL PYMIGN(3)
0788         ISUB=MINT(1)
0789  
0790 C...Study user-defined process: kinematics plus weight.
0791       ELSEIF(ISTSB.EQ.11) THEN
0792         IF(IDWTUP.GT.0.AND.XWGTUP.LT.0D0) CALL
0793      &  PYERRM(26,'(PYRAND:) Negative XWGTUP for user process')
0794         MSTI(51)=0
0795         IF(NUP.LE.0) THEN
0796           MINT(51)=2
0797           MSTI(51)=1
0798           IF(MINT(82).EQ.1) THEN
0799             NGEN(0,1)=NGEN(0,1)-1
0800             NGEN(ISUB,1)=NGEN(ISUB,1)-1
0801           ENDIF
0802           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
0803           RETURN
0804         ENDIF
0805  
0806 C...Extract cross section event weight.
0807         IF(IABS(IDWTUP).EQ.1.OR.IABS(IDWTUP).EQ.4) THEN
0808           SIGS=1D-9*XWGTUP
0809         ELSE
0810           SIGS=1D-9*XSECUP(KFPR(ISUB,1))
0811         ENDIF
0812         IF(IABS(IDWTUP).GE.1.AND.IABS(IDWTUP).LE.3) THEN
0813           VINT(97)=SIGN(1D0,XWGTUP)
0814         ELSE
0815           VINT(97)=1D-9*XWGTUP
0816         ENDIF
0817  
0818 C...Construct 'trivial' kinematical variables needed.
0819         KFL1=IDUP(1)
0820         KFL2=IDUP(2)
0821         VINT(41)=PUP(4,1)/EBMUP(1)
0822         VINT(42)=PUP(4,2)/EBMUP(2)
0823         VINT(21)=VINT(41)*VINT(42)
0824         VINT(22)=0.5D0*LOG(VINT(41)/VINT(42))
0825         VINT(44)=VINT(21)*VINT(2)
0826         VINT(43)=SQRT(MAX(0D0,VINT(44)))
0827         VINT(55)=SCALUP
0828         IF(SCALUP.LE.0D0) VINT(55)=VINT(43)
0829         VINT(56)=VINT(55)**2
0830         VINT(57)=AQEDUP
0831         VINT(58)=AQCDUP
0832  
0833 C...Construct other kinematical variables needed (approximately).
0834         VINT(23)=0D0
0835         VINT(26)=VINT(21)
0836         VINT(45)=-0.5D0*VINT(44)
0837         VINT(46)=-0.5D0*VINT(44)
0838         VINT(49)=VINT(43)
0839         VINT(50)=VINT(44)
0840         VINT(51)=VINT(55)
0841         VINT(52)=VINT(56)
0842         VINT(53)=VINT(55)
0843         VINT(54)=VINT(56)
0844         VINT(25)=0D0
0845         VINT(48)=0D0
0846         IF(ISTUP(1).NE.-1.OR.ISTUP(2).NE.-1) CALL PYERRM(26,
0847      &  '(PYRAND:) unacceptable ISTUP code for incoming particles')
0848         DO 280 IUP=3,NUP
0849           IF(ISTUP(IUP).LT.1.OR.ISTUP(IUP).GT.3) CALL PYERRM(26,
0850      &    '(PYRAND:) unacceptable ISTUP code for particles')
0851           IF(ISTUP(IUP).EQ.1) VINT(25)=VINT(25)+2D0*(PUP(5,IUP)**2+
0852      &    PUP(1,IUP)**2+PUP(2,IUP)**2)/VINT(2)
0853           IF(ISTUP(IUP).EQ.1) VINT(48)=VINT(48)+0.5D0*(PUP(1,IUP)**2+
0854      &    PUP(2,IUP)**2)
0855   280   CONTINUE
0856         VINT(47)=SQRT(VINT(48))
0857       ENDIF
0858  
0859 C...Choose azimuthal angle.
0860       VINT(24)=0D0
0861       IF(ISTSB.NE.11) VINT(24)=PARU(2)*PYR(0)
0862  
0863 C...Check against user cuts on kinematics at parton level.
0864       MINT(51)=0
0865       IF((ISUB.LE.90.OR.ISUB.GT.100).AND.ISTSB.LE.10) CALL PYKLIM(0)
0866       IF(MINT(51).NE.0) THEN
0867         IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
0868         IF(MFAIL.EQ.1) THEN
0869           MSTI(61)=1
0870           RETURN
0871         ENDIF
0872         GOTO 100
0873       ENDIF
0874       IF(MINT(82).EQ.1.AND.MSTP(141).GE.1.AND.ISTSB.LE.10) THEN
0875         MCUT=0
0876         IF(MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+MSUB(95).EQ.0)
0877      &  CALL PYKCUT(MCUT)
0878         IF(MCUT.NE.0) THEN
0879           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
0880           IF(MFAIL.EQ.1) THEN
0881             MSTI(61)=1
0882             RETURN
0883           ENDIF
0884           GOTO 100
0885         ENDIF
0886       ENDIF
0887  
0888 C...Calculate differential cross-section for different subprocesses.
0889       IF(ISTSB.LE.10) CALL PYSIGH(NCHN,SIGS)
0890       SIGSOR=SIGS
0891       SIGLPT=SIGT(0,0,5)*VINT(315)*VINT(316)
0892  
0893 C...Multiply cross section by lepton -> photon flux factor.
0894       IF(MINT(141).NE.0.OR.MINT(142).NE.0) THEN
0895         SIGS=WTGAGA*SIGS
0896         DO 290 ICHN=1,NCHN
0897           SIGH(ICHN)=WTGAGA*SIGH(ICHN)
0898   290   CONTINUE
0899         SIGLPT=WTGAGA*SIGLPT
0900       ENDIF
0901  
0902 C...Multiply cross-section by user-defined weights.
0903       IF(MSTP(173).EQ.1) THEN
0904         SIGS=PARP(173)*SIGS
0905         DO 300 ICHN=1,NCHN
0906           SIGH(ICHN)=PARP(173)*SIGH(ICHN)
0907   300   CONTINUE
0908         SIGLPT=PARP(173)*SIGLPT
0909       ENDIF
0910       WTXS=1D0
0911       SIGSWT=SIGS
0912       VINT(99)=1D0
0913       VINT(100)=1D0
0914       IF(MINT(82).EQ.1.AND.MSTP(142).GE.1) THEN
0915         IF(ISUB.NE.96.AND.MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+
0916      &  MSUB(95).EQ.0) CALL PYEVWT(WTXS)
0917         SIGSWT=WTXS*SIGS
0918         VINT(99)=WTXS
0919         IF(MSTP(142).EQ.1) VINT(100)=1D0/WTXS
0920       ENDIF
0921  
0922 C...Calculations for Monte Carlo estimate of all cross-sections.
0923       IF(MINT(82).EQ.1.AND.ISUB.LE.90.OR.ISUB.GE.96) THEN
0924         IF(MSTP(142).LE.1) THEN
0925           XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS
0926         ELSE
0927           XSEC(ISUB,2)=XSEC(ISUB,2)+SIGSWT
0928         ENDIF
0929       ELSEIF(MINT(82).EQ.1) THEN
0930         XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS
0931       ENDIF
0932       IF((ISUB.EQ.95.OR.ISUB.EQ.96).AND.LOOP2.EQ.1.AND.
0933      &MINT(82).EQ.1) XSEC(97,2)=XSEC(97,2)+SIGLPT
0934  
0935 C...Multiple interactions: store results of cross-section calculation.
0936       IF(MINT(50).EQ.1.AND.MSTP(82).GE.3) THEN
0937         VINT(153)=SIGSOR
0938         IF(MINT(35).LE.1) CALL PYMULT(4)
0939         IF(MINT(35).GE.2) CALL PYMIGN(4)
0940       ENDIF
0941  
0942 C...Ratio of actual to maximum cross section.
0943       IF(ISTSB.NE.11) THEN
0944         VIOL=SIGSWT/XSEC(ISUB,1)
0945         IF(ISUB.EQ.96.AND.MSTP(173).EQ.1) VIOL=VIOL/PARP(174)
0946       ELSEIF(IDWTUP.EQ.1.OR.IDWTUP.EQ.2) THEN
0947         VIOL=XWGTUP/XMAXUP(KFPR(ISUB,1))
0948       ELSEIF(IDWTUP.EQ.-1.OR.IDWTUP.EQ.-2) THEN
0949         VIOL=ABS(XWGTUP)/ABS(XMAXUP(KFPR(ISUB,1)))
0950       ELSE
0951         VIOL=1D0
0952       ENDIF
0953  
0954 C...Check that weight not negative.
0955       IF(MSTP(123).LE.0) THEN
0956         IF(VIOL.LT.-1D-3) THEN
0957           WRITE(MSTU(11),5000) VIOL,NGEN(0,3)+1
0958           IF(MSTP(122).GE.1) WRITE(MSTU(11),5100) ISUB,VINT(21),
0959      &    VINT(22),VINT(23),VINT(26)
0960           CALL PYSTOP(2)
0961         ENDIF
0962       ELSE
0963         IF(VIOL.LT.MIN(-1D-3,VINT(109))) THEN
0964           VINT(109)=VIOL
0965           IF(MSTP(123).LE.2) WRITE(MSTU(11),5200) VIOL,NGEN(0,3)+1
0966           IF(MSTP(122).GE.1) WRITE(MSTU(11),5100) ISUB,VINT(21),
0967      &    VINT(22),VINT(23),VINT(26)
0968         ENDIF
0969       ENDIF
0970  
0971 C...Weighting using estimate of maximum of differential cross-section.
0972       RATND=1D0
0973       IF(MFAIL.EQ.0.AND.ISUB.NE.95.AND.ISUB.NE.96) THEN
0974         IF(VIOL.LT.PYR(0)) THEN
0975           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
0976           IF(ISUB.GE.91.AND.ISUB.LE.94) ISUB=0
0977           GOTO 100
0978         ENDIF
0979       ELSEIF(MFAIL.EQ.0) THEN
0980         RATND=SIGLPT/XSEC(95,1)
0981         VIOL=VIOL/RATND
0982         IF(LOOP2.EQ.1.AND.RATND.LT.PYR(0)) THEN
0983           IF(VIOL.GT.PYR(0).AND.MINT(82).EQ.1.AND.MSUB(95).EQ.1.AND.
0984      &    (ISUB.LE.90.OR.ISUB.GE.95)) NGEN(95,1)=NGEN(95,1)+MINT(143)
0985           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
0986           ISUB=0
0987           GOTO 100
0988         ENDIF
0989         IF(VIOL.LT.PYR(0)) THEN
0990           GOTO 140
0991         ENDIF
0992       ELSEIF(ISUB.NE.95.AND.ISUB.NE.96) THEN
0993         IF(VIOL.LT.PYR(0)) THEN
0994           MSTI(61)=1
0995           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
0996           RETURN
0997         ENDIF
0998       ELSE
0999         RATND=SIGLPT/XSEC(95,1)
1000         IF(LOOP.EQ.1.AND.RATND.LT.PYR(0)) THEN
1001           MSTI(61)=1
1002           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
1003           RETURN
1004         ENDIF
1005         VIOL=VIOL/RATND
1006         IF(VIOL.LT.PYR(0)) THEN
1007           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
1008           GOTO 100
1009         ENDIF
1010       ENDIF
1011  
1012 C...Check for possible violation of estimated maximum of differential
1013 C...cross-section used in weighting.
1014       IF(MSTP(123).LE.0) THEN
1015         IF(VIOL.GT.1D0) THEN
1016           WRITE(MSTU(11),5300) VIOL,NGEN(0,3)+1
1017           IF(MSTP(122).GE.2) WRITE(MSTU(11),5100) ISUB,VINT(21),
1018      &    VINT(22),VINT(23),VINT(26)
1019           CALL PYSTOP(2)
1020         ENDIF
1021       ELSEIF(MSTP(123).EQ.1) THEN
1022         IF(VIOL.GT.VINT(108)) THEN
1023           VINT(108)=VIOL
1024           IF(VIOL.GT.1.0001D0) THEN
1025             MINT(10)=1
1026             WRITE(MSTU(11),5400) VIOL,NGEN(0,3)+1
1027             IF(MSTP(122).GE.2) WRITE(MSTU(11),5100) ISUB,VINT(21),
1028      &      VINT(22),VINT(23),VINT(26)
1029           ENDIF
1030         ENDIF
1031       ELSEIF(VIOL.GT.VINT(108)) THEN
1032         VINT(108)=VIOL
1033         IF(VIOL.GT.1D0) THEN
1034           MINT(10)=1
1035           IF(MSTP(123).EQ.2) WRITE(MSTU(11),5400) VIOL,NGEN(0,3)+1
1036           IF(ISTSB.EQ.11.AND.(IABS(IDWTUP).EQ.1.OR.IABS(IDWTUP).EQ.2))
1037      &    THEN
1038             XMAXUP(KFPR(ISUB,1))=VIOL*XMAXUP(KFPR(ISUB,1))
1039             IF(KFPR(ISUB,1).LE.9) THEN
1040               IF(MSTP(123).EQ.2) WRITE(MSTU(11),5800) KFPR(ISUB,1),
1041      &        XMAXUP(KFPR(ISUB,1))
1042             ELSEIF(KFPR(ISUB,1).LE.99) THEN
1043               IF(MSTP(123).EQ.2) WRITE(MSTU(11),5900) KFPR(ISUB,1),
1044      &        XMAXUP(KFPR(ISUB,1))
1045             ELSE
1046               IF(MSTP(123).EQ.2) WRITE(MSTU(11),6000) KFPR(ISUB,1),
1047      &        XMAXUP(KFPR(ISUB,1))
1048             ENDIF
1049           ENDIF
1050           IF(ISTSB.NE.11.OR.IABS(IDWTUP).EQ.1) THEN
1051             XDIF=XSEC(ISUB,1)*(VIOL-1D0)
1052             XSEC(ISUB,1)=XSEC(ISUB,1)+XDIF
1053             IF(MSUB(ISUB).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96))
1054      &      XSEC(0,1)=XSEC(0,1)+XDIF
1055             IF(MSTP(122).GE.2) WRITE(MSTU(11),5100) ISUB,VINT(21),
1056      &      VINT(22),VINT(23),VINT(26)
1057             IF(ISUB.LE.9) THEN
1058               IF(MSTP(123).EQ.2) WRITE(MSTU(11),5500) ISUB,XSEC(ISUB,1)
1059             ELSEIF(ISUB.LE.99) THEN
1060               IF(MSTP(123).EQ.2) WRITE(MSTU(11),5600) ISUB,XSEC(ISUB,1)
1061             ELSE
1062               IF(MSTP(123).EQ.2) WRITE(MSTU(11),5700) ISUB,XSEC(ISUB,1)
1063             ENDIF
1064           ENDIF
1065           VINT(108)=1D0
1066         ENDIF
1067       ENDIF
1068  
1069 C...Multiple interactions: choose impact parameter (if not already done).
1070       IF(MINT(39).EQ.0) VINT(148)=1D0
1071       IF(MINT(50).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GE.96).AND.
1072      &MSTP(82).GE.3) THEN
1073         IF(MINT(35).LE.1) CALL PYMULT(5)
1074         IF(MINT(35).GE.2) CALL PYMIGN(5)
1075         IF(VINT(150).LT.PYR(0)) THEN
1076           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
1077           IF(MFAIL.EQ.1) THEN
1078             MSTI(61)=1
1079             RETURN
1080           ENDIF
1081           GOTO 100
1082         ENDIF
1083       ENDIF
1084       IF(MINT(82).EQ.1) NGEN(0,2)=NGEN(0,2)+1
1085       IF(MINT(82).EQ.1.AND.MSUB(95).EQ.1) THEN
1086         IF(ISUB.LE.90.OR.ISUB.GE.95) NGEN(95,1)=NGEN(95,1)+MINT(143)
1087         IF(ISUB.LE.90.OR.ISUB.GE.96) NGEN(96,2)=NGEN(96,2)+1
1088       ENDIF
1089       IF(ISUB.LE.90.OR.ISUB.GE.96) MINT(31)=MINT(31)+1
1090  
1091 C...Choose flavour of reacting partons (and subprocess).
1092       IF(ISTSB.GE.11) GOTO 320
1093       RSIGS=SIGS*PYR(0)
1094       QT2=VINT(48)
1095       RQQBAR=PARP(87)*(1D0-(QT2/(QT2+(PARP(88)*PARP(82)*
1096      &(VINT(1)/PARP(89))**PARP(90))**2))**2)
1097       IF(ISUB.NE.95.AND.(ISUB.NE.96.OR.MSTP(82).LE.1.OR.
1098      &PYR(0).GT.RQQBAR)) THEN
1099         DO 310 ICHN=1,NCHN
1100           KFL1=ISIG(ICHN,1)
1101           KFL2=ISIG(ICHN,2)
1102           MINT(2)=ISIG(ICHN,3)
1103           RSIGS=RSIGS-SIGH(ICHN)
1104           IF(RSIGS.LE.0D0) GOTO 320
1105   310   CONTINUE
1106  
1107 C...Multiple interactions: choose qqbar preferentially at small pT.
1108       ELSEIF(ISUB.EQ.96) THEN
1109         MINT(105)=MINT(103)
1110         MINT(109)=MINT(107)
1111         CALL PYSPLI(MINT(11),21,KFL1,KFLDUM)
1112         MINT(105)=MINT(104)
1113         MINT(109)=MINT(108)
1114         CALL PYSPLI(MINT(12),21,KFL2,KFLDUM)
1115         MINT(1)=11
1116         MINT(2)=1
1117         IF(KFL1.EQ.KFL2.AND.PYR(0).LT.0.5D0) MINT(2)=2
1118  
1119 C...Low-pT: choose string drawing configuration.
1120       ELSE
1121         KFL1=21
1122         KFL2=21
1123         RSIGS=6D0*PYR(0)
1124         MINT(2)=1
1125         IF(RSIGS.GT.1D0) MINT(2)=2
1126         IF(RSIGS.GT.2D0) MINT(2)=3
1127       ENDIF
1128  
1129 C...Reassign QCD process. Partons before initial state radiation.
1130   320 IF(MINT(2).GT.10) THEN
1131         MINT(1)=MINT(2)/10
1132         MINT(2)=MOD(MINT(2),10)
1133       ENDIF
1134       IF(MINT(82).EQ.1.AND.MSTP(111).GE.0) NGEN(MINT(1),2)=
1135      &NGEN(MINT(1),2)+1
1136       MINT(15)=KFL1
1137       MINT(16)=KFL2
1138       MINT(13)=MINT(15)
1139       MINT(14)=MINT(16)
1140       VINT(141)=VINT(41)
1141       VINT(142)=VINT(42)
1142       VINT(151)=0D0
1143       VINT(152)=0D0
1144  
1145 C...Calculate x value of photon for parton inside photon inside e.
1146       DO 350 JT=1,2
1147         MINT(18+JT)=0
1148         VINT(154+JT)=0D0
1149         MSPLI=0
1150         IF(JT.EQ.1.AND.MINT(43).LE.2) MSPLI=1
1151         IF(JT.EQ.2.AND.MOD(MINT(43),2).EQ.1) MSPLI=1
1152         IF(IABS(MINT(14+JT)).LE.8.OR.MINT(14+JT).EQ.21) MSPLI=MSPLI+1
1153         IF(MSPLI.EQ.2) THEN
1154           KFLH=MINT(14+JT)
1155           XHRD=VINT(140+JT)
1156           Q2HRD=VINT(54)
1157           MINT(105)=MINT(102+JT)
1158           MINT(109)=MINT(106+JT)
1159           VINT(120)=VINT(2+JT)
1160           IF(MSTP(57).LE.1) THEN
1161             CALL PYPDFU(22,XHRD,Q2HRD,XPQ)
1162           ELSE
1163             CALL PYPDFL(22,XHRD,Q2HRD,XPQ)
1164           ENDIF
1165           WTMX=4D0*XPQ(KFLH)
1166           IF(MSTP(13).EQ.2) THEN
1167             Q2PMS=Q2HRD/PMAS(11,1)**2
1168             WTMX=WTMX*LOG(MAX(2D0,Q2PMS*(1D0-XHRD)/XHRD**2))
1169           ENDIF
1170   330     XE=XHRD**PYR(0)
1171           XG=MIN(1D0-1D-10,XHRD/XE)
1172           IF(MSTP(57).LE.1) THEN
1173             CALL PYPDFU(22,XG,Q2HRD,XPQ)
1174           ELSE
1175             CALL PYPDFL(22,XG,Q2HRD,XPQ)
1176           ENDIF
1177           WT=(1D0+(1D0-XE)**2)*XPQ(KFLH)
1178           IF(MSTP(13).EQ.2) WT=WT*LOG(MAX(2D0,Q2PMS*(1D0-XE)/XE**2))
1179           IF(WT.LT.PYR(0)*WTMX) GOTO 330
1180           MINT(18+JT)=1
1181           VINT(154+JT)=XE
1182           DO 340 KFLS=-25,25
1183             XSFX(JT,KFLS)=XPQ(KFLS)
1184   340     CONTINUE
1185         ENDIF
1186   350 CONTINUE
1187  
1188 C...Pick scale where photon is resolved.
1189       Q0S=PARP(15)**2
1190       Q1S=VINT(154)**2
1191       VINT(283)=0D0
1192       IF(MINT(107).EQ.3) THEN
1193         IF(MSTP(66).EQ.1) THEN
1194           VINT(283)=Q0S*(VINT(54)/Q0S)**PYR(0)
1195         ELSEIF(MSTP(66).EQ.2) THEN
1196           PS=VINT(3)**2
1197           Q2EFF=VINT(54)*((Q0S+PS)/(VINT(54)+PS))*
1198      &    EXP(PS*(VINT(54)-Q0S)/((VINT(54)+PS)*(Q0S+PS)))
1199           Q2INT=SQRT(Q0S*Q2EFF)
1200           VINT(283)=Q2INT*(VINT(54)/Q2INT)**PYR(0)
1201         ELSEIF(MSTP(66).EQ.3) THEN
1202           VINT(283)=Q0S*(Q1S/Q0S)**PYR(0)
1203         ELSEIF(MSTP(66).GE.4) THEN
1204           PS=0.25D0*VINT(3)**2
1205           VINT(283)=(Q0S+PS)*(Q1S+PS)/
1206      &    (Q0S+PYR(0)*(Q1S-Q0S)+PS)-PS
1207         ENDIF
1208       ENDIF
1209       VINT(284)=0D0
1210       IF(MINT(108).EQ.3) THEN
1211         IF(MSTP(66).EQ.1) THEN
1212           VINT(284)=Q0S*(VINT(54)/Q0S)**PYR(0)
1213         ELSEIF(MSTP(66).EQ.2) THEN
1214           PS=VINT(4)**2
1215           Q2EFF=VINT(54)*((Q0S+PS)/(VINT(54)+PS))*
1216      &    EXP(PS*(VINT(54)-Q0S)/((VINT(54)+PS)*(Q0S+PS)))
1217           Q2INT=SQRT(Q0S*Q2EFF)
1218           VINT(284)=Q2INT*(VINT(54)/Q2INT)**PYR(0)
1219         ELSEIF(MSTP(66).EQ.3) THEN
1220           VINT(284)=Q0S*(Q1S/Q0S)**PYR(0)
1221         ELSEIF(MSTP(66).GE.4) THEN
1222           PS=0.25D0*VINT(4)**2
1223           VINT(284)=(Q0S+PS)*(Q1S+PS)/
1224      &    (Q0S+PYR(0)*(Q1S-Q0S)+PS)-PS
1225         ENDIF
1226       ENDIF
1227       IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
1228  
1229 C...Format statements for differential cross-section maximum violations.
1230  5000 FORMAT(/1X,'Error: negative cross-section fraction',1P,D11.3,1X,
1231      &'in event',1X,I7,'D0'/1X,'Execution stopped!')
1232  5100 FORMAT(1X,'ISUB = ',I3,'; Point of violation:'/1X,'tau =',1P,
1233      &D11.3,', y* =',D11.3,', cthe = ',0P,F11.7,', tau'' =',1P,D11.3)
1234  5200 FORMAT(/1X,'Warning: negative cross-section fraction',1P,D11.3,1X,
1235      &'in event',1X,I7)
1236  5300 FORMAT(/1X,'Error: maximum violated by',1P,D11.3,1X,
1237      &'in event',1X,I7,'D0'/1X,'Execution stopped!')
1238  5400 FORMAT(/1X,'Advisory warning: maximum violated by',1P,D11.3,1X,
1239      &'in event',1X,I7)
1240  5500 FORMAT(1X,'XSEC(',I1,',1) increased to',1P,D11.3)
1241  5600 FORMAT(1X,'XSEC(',I2,',1) increased to',1P,D11.3)
1242  5700 FORMAT(1X,'XSEC(',I3,',1) increased to',1P,D11.3)
1243  5800 FORMAT(1X,'XMAXUP(',I1,') increased to',1P,D11.3)
1244  5900 FORMAT(1X,'XMAXUP(',I2,') increased to',1P,D11.3)
1245  6000 FORMAT(1X,'XMAXUP(',I3,') increased to',1P,D11.3)
1246  
1247       RETURN
1248       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!