pythiaeRHIC/pythia/pyevnt.f

0001
0002 C*********************************************************************
0003
0004 C...PYEVNT
0005 C...Administers the generation of a high-pT event via calls to
0006 C...a number of subroutines.
0007
0008       SUBROUTINE PYEVNT
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 C...Commonblocks.
0015       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0016       COMMON/PYCTAG/NCT,MCT(4000,2)
0017       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0018       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0019       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0020       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0021       COMMON/PYINT1/MINT(400),VINT(400)
0022       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0023       COMMON/PYINT4/MWID(500),WIDS(500,5)
0024       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
0025       SAVE /PYJETS/,/PYDAT1/,/PYCTAG/,/PYDAT2/,/PYDAT3/,/PYPARS/,
0026      &/PYINT1/,/PYINT2/,/PYINT4/,/PYINT5/
0027 C...Local array.
0028       DIMENSION VTX(4)
0029
0030 C...Optionally let PYEVNW do the whole job.
0031       IF(MSTP(81).GE.20) THEN
0032         CALL PYEVNW
0033         RETURN
0034       ENDIF
0035
0036 C...Stop if no subprocesses on.
0037       IF(MINT(121).EQ.1.AND.MSTI(53).EQ.1) THEN
0038         WRITE(MSTU(11),5100)
0039         CALL PYSTOP(1)
0040       ENDIF
0041
0042 C...Initial values for some counters.
0043       MSTU(1)=0
0044       MSTU(2)=0
0045       N=0
0046       MINT(5)=MINT(5)+1
0047       MINT(7)=0
0048       MINT(8)=0
0049       MINT(30)=0
0050       MINT(83)=0
0051       MINT(84)=MSTP(126)
0052       MSTU(24)=0
0053       MSTU70=0
0054       MSTJ14=MSTJ(14)
0055 C...Normally, use K(I,4:5) colour info rather than /PYCTAG/.
0056       NCT=0
0057       MINT(33)=0
0058
0059 C...Let called routines know call is from PYEVNT (not PYEVNW).
0060       MINT(35)=1
0061       IF (MSTP(81).GE.10) MINT(35)=2
0062
0063 C...If variable energies: redo incoming kinematics and cross-section.
0064       MSTI(61)=0
0065       IF(MSTP(171).EQ.1) THEN
0066         CALL PYINKI(1)
0067         IF(MSTI(61).EQ.1) THEN
0068           MINT(5)=MINT(5)-1
0069           RETURN
0070         ENDIF
0071         IF(MINT(121).GT.1) CALL PYSAVE(3,1)
0072         CALL PYXTOT
0073       ENDIF
0074
0075 C...Loop over number of pileup events; check space left.
0076       IF(MSTP(131).LE.0) THEN
0077         NPILE=1
0078       ELSE
0079         CALL PYPILE(2)
0080         NPILE=MINT(81)
0081       ENDIF
0082       DO 270 IPILE=1,NPILE
0083         IF(MINT(84)+100.GE.MSTU(4)) THEN
0084           CALL PYERRM(11,
0085      &    '(PYEVNT:) no more space in PYJETS for pileup events')
0086           IF(MSTU(21).GE.1) GOTO 280
0087         ENDIF
0088         MINT(82)=IPILE
0089
0090 C...Generate variables of hard scattering.
0091         MINT(51)=0
0092         MSTI(52)=0
0093   100   CONTINUE
0094         IF(MINT(51).NE.0.OR.MSTU(24).NE.0) MSTI(52)=MSTI(52)+1
0095         MINT(31)=0
0096         MINT(39)=0
0097         MINT(51)=0
0098         MINT(57)=0
0099         CALL PYRAND
0100         IF(MSTI(61).EQ.1) THEN
0101           MINT(5)=MINT(5)-1
0102           RETURN
0103         ENDIF
0104         IF(MINT(51).EQ.2) RETURN
0105         ISUB=MINT(1)
0106         IF(MSTP(111).EQ.-1) GOTO 260
0107
0108 C...Loopback point if PYPREP fails, especially for junction topologies.
0109         NPREP=0
0110         MNT31S=MINT(31)
0111   110   NPREP=NPREP+1
0112         MINT(31)=MNT31S
0113
0114         IF((ISUB.LE.90.OR.ISUB.GE.95).AND.ISUB.NE.99) THEN
0115 C...Hard scattering (including low-pT):
0116 C...reconstruct kinematics and colour flow of hard scattering.
0117           MINT31=MINT(31)
0118   120     MINT(31)=MINT31
0119           MINT(51)=0
0120           CALL PYSCAT
0121           IF(MINT(51).EQ.1) GOTO 100
0122           IPU1=MINT(84)+1
0123           IPU2=MINT(84)+2
0124           IF(ISUB.EQ.95) GOTO 140
0125
0126 C...Reset statistics on activity in event.
0127         DO 130 J=351,359
0128           MINT(J)=0
0129           VINT(J)=0D0
0130   130   CONTINUE
0131
0132 C...Showering of initial state partons (optional).
0133           NFIN=N
0134           ALAMSV=PARJ(81)
0135           PARJ(81)=PARP(72)
0136           IF(MSTP(61).GE.1.AND.MINT(47).GE.2.AND.MINT(111).NE.12)
0137      &    CALL PYSSPA(IPU1,IPU2)
0138           PARJ(81)=ALAMSV
0139           IF(MINT(51).EQ.1) GOTO 100
0140
0141 C...Showering of final state partons (optional).
0142           ALAMSV=PARJ(81)
0143           PARJ(81)=PARP(72)
0144           IF(MSTP(71).GE.1.AND.ISET(ISUB).GE.2.AND.ISET(ISUB).LE.10)
0145      &    THEN
0146             IPU3=MINT(84)+3
0147             IPU4=MINT(84)+4
0148             IF(ISET(ISUB).EQ.5) IPU4=-3
0149             QMAX=VINT(55)
0150             IF(ISET(ISUB).EQ.2) QMAX=SQRT(PARP(71))*VINT(55)
0151             CALL PYSHOW(IPU3,IPU4,QMAX)
0152           ELSEIF(ISET(ISUB).EQ.11) THEN
0153             CALL PYADSH(NFIN)
0154           ENDIF
0155           PARJ(81)=ALAMSV
0156
0157 C...Allow possibility for user to abort event generation.
0158           IVETO=0
0159           IF(IPILE.EQ.1.AND.MSTP(143).EQ.1) CALL PYVETO(IVETO)
0160           IF(IVETO.EQ.1) GOTO 100
0161
0162 C...Decay of final state resonances.
0163           MINT(32)=0
0164           IF(MSTP(41).GE.1.AND.ISET(ISUB).LE.10) CALL PYRESD(0)
0165           IF(MINT(51).EQ.1) GOTO 100
0166           MINT(52)=N
0167
0168
0169 C...Multiple interactions - PYTHIA 6.3 intermediate style.
0170   140     IF(MSTP(81).GE.10.AND.MINT(50).EQ.1) THEN
0171             IF(ISUB.EQ.95) MINT(31)=MINT(31)+1
0172             CALL PYMIGN(6)
0173             IF(MINT(51).EQ.1) GOTO 100
0174             MINT(53)=N
0175
0176 C...Beam remnant flavour and colour assignments - new scheme.
0177             CALL PYMIHK
0178             IF(MINT(51).EQ.1.AND.MINT(57).GE.1.AND.MINT(57).LE.5)
0179      &      GOTO 120
0180             IF(MINT(51).EQ.1) GOTO 100
0181
0182 C...Primordial kT and beam remnant momentum sharing - new scheme.
0183             CALL PYMIRM
0184             IF(MINT(51).EQ.1.AND.MINT(57).GE.1.AND.MINT(57).LE.5)
0185      &      GOTO 120
0186             IF(MINT(51).EQ.1) GOTO 100
0187             IF(ISUB.EQ.95) MINT(31)=MINT(31)-1
0188
0189 C...Multiple interactions - PYTHIA 6.2 style.
0190           ELSEIF(MINT(111).NE.12) THEN
0191             IF (MSTP(81).GE.1.AND.MINT(50).EQ.1.AND.ISUB.NE.95) THEN
0192               CALL PYMULT(6)
0193               MINT(53)=N
0194             ENDIF
0195
0196 C...Hadron remnants and primordial kT.
0197             CALL PYREMN(IPU1,IPU2)
0198             IF(MINT(51).EQ.1.AND.MINT(57).GE.1.AND.MINT(57).LE.5) GOTO
0199      &           110
0200             IF(MINT(51).EQ.1) GOTO 100
0201           ENDIF
0202
0203         ELSEIF(ISUB.NE.99) THEN
0204 C...Diffractive and elastic scattering.
0205           CALL PYDIFF
0206
0207         ELSE
0208 C...DIS scattering (photon flux external).
0209           CALL PYDISG
0210           IF(MINT(51).EQ.1) GOTO 100
0211         ENDIF
0212
0213 C...Check that no odd resonance left undecayed.
0214         MINT(54)=N
0215         IF(MSTP(111).GE.1) THEN
0216           NFIX=N
0217           DO 150 I=MINT(84)+1,NFIX
0218             IF(K(I,1).GE.1.AND.K(I,1).LE.10.AND.K(I,2).NE.21.AND.
0219      &      K(I,2).NE.22) THEN
0220               KCA=PYCOMP(K(I,2))
0221               IF(MWID(KCA).NE.0.AND.MDCY(KCA,1).GE.1) THEN
0222                 CALL PYRESD(I)
0223                 IF(MINT(51).EQ.1) GOTO 100
0224               ENDIF
0225             ENDIF
0226   150     CONTINUE
0227         ENDIF
0228
0229 C...Boost hadronic subsystem to overall rest frame.
0230 C..(Only relevant when photon inside lepton beam.)
0231         IF(MINT(141).NE.0.OR.MINT(142).NE.0) CALL PYGAGA(4,WTGAGA)
0232
0233 C...Recalculate energies from momenta and masses (if desired).
0234         IF(MSTP(113).GE.1) THEN
0235           DO 160 I=MINT(83)+1,N
0236             IF(K(I,1).GT.0.AND.K(I,1).LE.10) P(I,4)=SQRT(P(I,1)**2+
0237      &      P(I,2)**2+P(I,3)**2+P(I,5)**2)
0238   160     CONTINUE
0239           NRECAL=N
0240         ENDIF
0241
0242 C...Colour reconnection before string formation
0243         IF (MSTP(95).GE.2) CALL PYFSCR(MINT(84)+1)
0244
0245 C...Rearrange partons along strings, check invariant mass cuts.
0246         MSTU(28)=0
0247         IF(MSTP(111).LE.0) MSTJ(14)=-1
0248         CALL PYPREP(MINT(84)+1)
0249         MSTJ(14)=MSTJ14
0250         IF(MINT(51).EQ.1.AND.MSTU(24).EQ.1) THEN
0251           MSTU(24)=0
0252           GOTO 100
0253         ENDIF
0254         IF (MINT(51).EQ.1.AND.NPREP.LE.5) GOTO 110
0255         IF (MINT(51).EQ.1) GOTO 100
0256         IF(MSTP(112).EQ.1.AND.MSTU(28).EQ.3) GOTO 100
0257         IF(MSTP(125).EQ.0.OR.MSTP(125).EQ.1) THEN
0258           DO 190 I=MINT(84)+1,N
0259             IF(K(I,2).EQ.94) THEN
0260               DO 180 I1=I+1,MIN(N,I+10)
0261                 IF(K(I1,3).EQ.I) THEN
0262                   K(I1,3)=MOD(K(I1,4)/MSTU(5),MSTU(5))
0263                   IF(K(I1,3).EQ.0) THEN
0264                     DO 170 II=MINT(84)+1,I-1
0265                         IF(K(II,2).EQ.K(I1,2)) THEN
0266                           IF(MOD(K(II,4),MSTU(5)).EQ.I1.OR.
0267      &                    MOD(K(II,5),MSTU(5)).EQ.I1) K(I1,3)=II
0268                         ENDIF
0269   170               CONTINUE
0270                     IF(K(I+1,3).EQ.0) K(I+1,3)=K(I,3)
0271                   ENDIF
0272                 ENDIF
0273   180         CONTINUE
0274             ENDIF
0275   190     CONTINUE
0276           CALL PYEDIT(12)
0277           CALL PYEDIT(14)
0278           IF(MSTP(125).EQ.0) CALL PYEDIT(15)
0279           IF(MSTP(125).EQ.0) MINT(4)=0
0280           DO 210 I=MINT(83)+1,N
0281             IF(K(I,1).EQ.11.AND.K(I,4).EQ.0.AND.K(I,5).EQ.0) THEN
0282               DO 200 I1=I+1,N
0283                 IF(K(I1,3).EQ.I.AND.K(I,4).EQ.0) K(I,4)=I1
0284                 IF(K(I1,3).EQ.I) K(I,5)=I1
0285   200         CONTINUE
0286             ENDIF
0287   210     CONTINUE
0288         ENDIF
0289
0290 C...Introduce separators between sections in PYLIST event listing.
0291         IF(IPILE.EQ.1.AND.MSTP(125).LE.0) THEN
0292           MSTU70=1
0293           MSTU(71)=N
0294         ELSEIF(IPILE.EQ.1) THEN
0295           MSTU70=3
0296           MSTU(71)=2
0297           MSTU(72)=MINT(4)
0298           MSTU(73)=N
0299         ENDIF
0300
0301 C...Go back to lab frame (needed for vertices, also in fragmentation).
0302         CALL PYFRAM(1)
0303
0304 C...Set nonvanishing production vertex (optional).
0305         IF(MSTP(151).EQ.1) THEN
0306           DO 220 J=1,4
0307             VTX(J)=PARP(150+J)*SQRT(-2D0*LOG(MAX(1D-10,PYR(0))))*
0308      &      SIN(PARU(2)*PYR(0))
0309   220     CONTINUE
0310           DO 240 I=MINT(83)+1,N
0311             DO 230 J=1,4
0312               V(I,J)=V(I,J)+VTX(J)
0313   230       CONTINUE
0314   240     CONTINUE
0315         ENDIF
0316
0317 C...Perform hadronization (if desired).
0318         IF(MSTP(111).GE.1) THEN
0319           CALL PYEXEC
0320           IF(MSTU(24).NE.0) GOTO 100
0321         ENDIF
0322         IF(MSTP(113).GE.1) THEN
0323           DO 250 I=NRECAL,N
0324             IF(P(I,5).GT.0D0) P(I,4)=SQRT(P(I,1)**2+
0325      &      P(I,2)**2+P(I,3)**2+P(I,5)**2)
0326   250     CONTINUE
0327         ENDIF
0328         IF(MSTP(125).EQ.0.OR.MSTP(125).EQ.1) CALL PYEDIT(14)
0329
0330 C...Store event information and calculate Monte Carlo estimates of
0331 C...subprocess cross-sections.
0332   260   IF(IPILE.EQ.1) CALL PYDOCU
0333
0334 C...Set counters for current pileup event and loop to next one.
0335         MSTI(41)=IPILE
0336         IF(IPILE.GE.2.AND.IPILE.LE.10) MSTI(40+IPILE)=ISUB
0337         IF(MSTU70.LT.10) THEN
0338           MSTU70=MSTU70+1
0339           MSTU(70+MSTU70)=N
0340         ENDIF
0341         MINT(83)=N
0342         MINT(84)=N+MSTP(126)
0343         IF(IPILE.LT.NPILE) CALL PYFRAM(2)
0344   270 CONTINUE
0345
0346 C...Generic information on pileup events. Reconstruct missing history.
0347       IF(MSTP(131).EQ.1.AND.MSTP(133).GE.1) THEN
0348         PARI(91)=VINT(132)
0349         PARI(92)=VINT(133)
0350         PARI(93)=VINT(134)
0351         IF(MSTP(133).GE.2) PARI(93)=PARI(93)*XSEC(0,3)/VINT(131)
0352       ENDIF
0353       CALL PYEDIT(16)
0354
0355 C...Transform to the desired coordinate frame.
0356   280 CALL PYFRAM(MSTP(124))
0357       MSTU(70)=MSTU70
0358       PARU(21)=VINT(1)
0359
0360 C...Error messages
0361  5100 FORMAT(1X,'Error: no subprocess switched on.'/
0362      &1X,'Execution stopped.')
0363
0364       RETURN
0365       END