pythiaeRHIC/pythia/pyupev.f

0001
0002 C*********************************************************************
0003
0004 C...PYUPEV
0005 C...Administers the hard-process generation required for output to the
0006 C...Les Houches event record.
0007
0008       SUBROUTINE PYUPEV
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014
0015 C...Commonblocks.
0016       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0017       COMMON/PYCTAG/NCT,MCT(4000,2)
0018       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0019       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0020       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0021       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0022       COMMON/PYINT1/MINT(400),VINT(400)
0023       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0024       COMMON/PYINT4/MWID(500),WIDS(500,5)
0025       SAVE /PYJETS/,/PYCTAG/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYPARS/,
0026      &/PYINT1/,/PYINT2/,/PYINT4/
0027
0028 C...HEPEUP for output.
0029       INTEGER MAXNUP
0030       PARAMETER (MAXNUP=500)
0031       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
0032       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
0033       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
0034      &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
0035      &VTIMUP(MAXNUP),SPINUP(MAXNUP)
0036       SAVE /HEPEUP/
0037
0038 C...Stop if no subprocesses on.
0039       IF(MINT(121).EQ.1.AND.MSTI(53).EQ.1) THEN
0040         WRITE(MSTU(11),5100)
0041         STOP
0042       ENDIF
0043
0044 C...Special flags for hard-process generation only.
0045       MSTP71=MSTP(71)
0046       MSTP(71)=0
0047       MST128=MSTP(128)
0048       MSTP(128)=1
0049
0050 C...Initial values for some counters.
0051       N=0
0052       MINT(5)=MINT(5)+1
0053       MINT(7)=0
0054       MINT(8)=0
0055       MINT(30)=0
0056       MINT(83)=0
0057       MINT(84)=MSTP(126)
0058       MSTU(24)=0
0059       MSTU70=0
0060       MSTJ14=MSTJ(14)
0061 C...Normally, use K(I,4:5) colour info rather than /PYCTAG/.
0062       MINT(33)=0
0063
0064 C...If variable energies: redo incoming kinematics and cross-section.
0065       MSTI(61)=0
0066       IF(MSTP(171).EQ.1) THEN
0067         CALL PYINKI(1)
0068         IF(MSTI(61).EQ.1) THEN
0069           MINT(5)=MINT(5)-1
0070           RETURN
0071         ENDIF
0072         IF(MINT(121).GT.1) CALL PYSAVE(3,1)
0073         CALL PYXTOT
0074       ENDIF
0075
0076 C...Do not allow pileup events.
0077       MINT(82)=1
0078
0079 C...Generate variables of hard scattering.
0080       MINT(51)=0
0081       MSTI(52)=0
0082   100 CONTINUE
0083       IF(MINT(51).NE.0.OR.MSTU(24).NE.0) MSTI(52)=MSTI(52)+1
0084       MINT(31)=0
0085       MINT(51)=0
0086       MINT(57)=0
0087       CALL PYRAND
0088       IF(MSTI(61).EQ.1) THEN
0089         MINT(5)=MINT(5)-1
0090         RETURN
0091       ENDIF
0092       IF(MINT(51).EQ.2) RETURN
0093       ISUB=MINT(1)
0094
0095       IF((ISUB.LE.90.OR.ISUB.GE.95).AND.ISUB.NE.99) THEN
0096 C...Hard scattering (including low-pT):
0097 C...reconstruct kinematics and colour flow of hard scattering.
0098         MINT31=MINT(31)
0099   110   MINT(31)=MINT31
0100         MINT(51)=0
0101         CALL PYSCAT
0102         IF(MINT(51).EQ.1) GOTO 100
0103         IPU1=MINT(84)+1
0104         IPU2=MINT(84)+2
0105
0106 C...Decay of final state resonances.
0107         MINT(32)=0
0108         IF(MSTP(41).GE.1.AND.ISET(ISUB).LE.10.AND.ISUB.NE.95)
0109      &  CALL PYRESD(0)
0110         IF(MINT(51).EQ.1) GOTO 100
0111         MINT(52)=N
0112
0113 C...Longitudinal boost of hard scattering.
0114         BETAZ=(VINT(41)-VINT(42))/(VINT(41)+VINT(42))
0115         CALL PYROBO(MINT(84)+1,N,0D0,0D0,0D0,0D0,BETAZ)
0116
0117       ELSEIF(ISUB.NE.99) THEN
0118 C...Diffractive and elastic scattering.
0119         CALL PYDIFF
0120
0121       ELSE
0122 C...DIS scattering (photon flux external).
0123         CALL PYDISG
0124         IF(MINT(51).EQ.1) GOTO 100
0125       ENDIF
0126
0127 C...Check that no odd resonance left undecayed.
0128       MINT(54)=N
0129       NFIX=N
0130       DO 120 I=MINT(84)+1,NFIX
0131         IF(K(I,1).GE.1.AND.K(I,1).LE.10.AND.K(I,2).NE.21.AND.
0132      &  K(I,2).NE.22) THEN
0133           KCA=PYCOMP(K(I,2))
0134           IF(MWID(KCA).NE.0.AND.MDCY(KCA,1).GE.1) THEN
0135             CALL PYRESD(I)
0136             IF(MINT(51).EQ.1) GOTO 100
0137           ENDIF
0138         ENDIF
0139   120 CONTINUE
0140
0141 C...Boost hadronic subsystem to overall rest frame.
0142 C..(Only relevant when photon inside lepton beam.)
0143       IF(MINT(141).NE.0.OR.MINT(142).NE.0) CALL PYGAGA(4,WTGAGA)
0144
0145 C...Store event information and calculate Monte Carlo estimates of
0146 C...subprocess cross-sections.
0147   130 CALL PYDOCU
0148
0149 C...Transform to the desired coordinate frame.
0150   140 CALL PYFRAM(MSTP(124))
0151       MSTU(70)=MSTU70
0152       PARU(21)=VINT(1)
0153
0154 C...Restore special flags for hard-process generation only.
0155       MSTP(71)=MSTP71
0156       MSTP(128)=MST128
0157
0158 C...Trace colour tags; convert to LHA style labels.
0159       NCT=100
0160       DO 150 I=MINT(84)+1,N
0161         MCT(I,1)=0
0162         MCT(I,2)=0
0163   150 CONTINUE
0164       DO 160 I=MINT(84)+1,N
0165         KQ=KCHG(PYCOMP(K(I,2)),2)*ISIGN(1,K(I,2))
0166         IF(K(I,1).EQ.3.OR.K(I,1).EQ.13.OR.K(I,1).EQ.14) THEN
0167           IF(K(I,4).NE.0.AND.(KQ.EQ.1.OR.KQ.EQ.2).AND.MCT(I,1).EQ.0)
0168      &    THEN
0169             IMO=MOD(K(I,4)/MSTU(5),MSTU(5))
0170             IDA=MOD(K(I,4),MSTU(5))
0171             IF(IMO.NE.0.AND.MOD(K(IMO,5)/MSTU(5),MSTU(5)).EQ.I.AND.
0172      &      MCT(IMO,2).NE.0) THEN
0173               MCT(I,1)=MCT(IMO,2)
0174             ELSEIF(IMO.NE.0.AND.MOD(K(IMO,4),MSTU(5)).EQ.I.AND.
0175      &      MCT(IMO,1).NE.0) THEN
0176               MCT(I,1)=MCT(IMO,1)
0177             ELSEIF(IDA.NE.0.AND.MOD(K(IDA,5),MSTU(5)).EQ.I.AND.
0178      &      MCT(IDA,2).NE.0) THEN
0179               MCT(I,1)=MCT(IDA,2)
0180             ELSE
0181               NCT=NCT+1
0182               MCT(I,1)=NCT
0183             ENDIF
0184           ENDIF
0185           IF(K(I,5).NE.0.AND.(KQ.EQ.-1.OR.KQ.EQ.2).AND.MCT(I,2).EQ.0)
0186      &    THEN
0187             IMO=MOD(K(I,5)/MSTU(5),MSTU(5))
0188             IDA=MOD(K(I,5),MSTU(5))
0189             IF(IMO.NE.0.AND.MOD(K(IMO,4)/MSTU(5),MSTU(5)).EQ.I.AND.
0190      &      MCT(IMO,1).NE.0) THEN
0191               MCT(I,2)=MCT(IMO,1)
0192             ELSEIF(IMO.NE.0.AND.MOD(K(IMO,5),MSTU(5)).EQ.I.AND.
0193      &      MCT(IMO,2).NE.0) THEN
0194               MCT(I,2)=MCT(IMO,2)
0195             ELSEIF(IDA.NE.0.AND.MOD(K(IDA,4),MSTU(5)).EQ.I.AND.
0196      &      MCT(IDA,1).NE.0) THEN
0197               MCT(I,2)=MCT(IDA,1)
0198             ELSE
0199               NCT=NCT+1
0200               MCT(I,2)=NCT
0201             ENDIF
0202           ENDIF
0203         ENDIF
0204   160 CONTINUE
0205
0206 C...Put event in HEPEUP commonblock.
0207       NUP=N-MINT(84)
0208       IDPRUP=MINT(1)
0209       XWGTUP=1D0
0210       SCALUP=VINT(53)
0211       AQEDUP=VINT(57)
0212       AQCDUP=VINT(58)
0213       DO 180 I=1,NUP
0214         IDUP(I)=K(I+MINT(84),2)
0215         IF(I.LE.2) THEN
0216           ISTUP(I)=-1
0217           MOTHUP(1,I)=0
0218           MOTHUP(2,I)=0
0219         ELSEIF(K(I+4,3).EQ.0) THEN
0220           ISTUP(I)=1
0221           MOTHUP(1,I)=1
0222           MOTHUP(2,I)=2
0223         ELSE
0224           ISTUP(I)=1
0225           MOTHUP(1,I)=K(I+MINT(84),3)-MINT(84)
0226           MOTHUP(2,I)=0
0227         ENDIF
0228         IF(I.GE.3.AND.K(I+MINT(84),3).GT.0)
0229      &  ISTUP(K(I+MINT(84),3)-MINT(84))=2
0230         ICOLUP(1,I)=MCT(I+MINT(84),1)
0231         ICOLUP(2,I)=MCT(I+MINT(84),2)
0232         DO 170 J=1,5
0233           PUP(J,I)=P(I+MINT(84),J)
0234   170   CONTINUE
0235         VTIMUP(I)=V(I,5)
0236         SPINUP(I)=9D0
0237   180 CONTINUE
0238
0239 C...Optionally write out event to disk. Minimal size for time/spin fields.
0240       IF(MSTP(162).GT.0) THEN
0241         WRITE(MSTP(162),5200) NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP
0242         DO 190 I=1,NUP
0243           IF(VTIMUP(I).EQ.0D0) THEN
0244             WRITE(MSTP(162),5300) IDUP(I),ISTUP(I),MOTHUP(1,I),
0245      &      MOTHUP(2,I),ICOLUP(1,I),ICOLUP(2,I),(PUP(J,I),J=1,5),
0246      &      ' 0. 9.'
0247           ELSE
0248             WRITE(MSTP(162),5400) IDUP(I),ISTUP(I),MOTHUP(1,I),
0249      &      MOTHUP(2,I),ICOLUP(1,I),ICOLUP(2,I),(PUP(J,I),J=1,5),
0250      &      VTIMUP(I),' 9.'
0251           ENDIF
0252   190   CONTINUE
0253
0254 C...Optional extra line with parton-density information.
0255         IF(MSTP(165).GE.1) WRITE(MSTP(162),5500) MSTI(15),MSTI(16),
0256      &  PARI(33),PARI(34),PARI(23),PARI(29),PARI(30)
0257       ENDIF
0258
0259 C...Error messages and other print formats.
0260  5100 FORMAT(1X,'Error: no subprocess switched on.'/
0261      &1X,'Execution stopped.')
0262  5200 FORMAT(1P,2I6,4E14.6)
0263  5300 FORMAT(1P,I8,5I5,5E18.10,A6)
0264  5400 FORMAT(1P,I8,5I5,5E18.10,E12.4,A3)
0265  5500 FORMAT(1P,'#pdf ',2I5,5E18.10)
0266
0267       RETURN
0268       END