pythiaeRHIC/pythia/pyexec.f

0001
0002 C*********************************************************************
0003
0004 C...PYEXEC
0005 C...Administrates the fragmentation and decay chain.
0006
0007       SUBROUTINE PYEXEC
0008
0009 C...Double precision and integer declarations.
0010       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0011       IMPLICIT INTEGER(I-N)
0012       INTEGER PYK,PYCHGE,PYCOMP
0013 C...Commonblocks.
0014       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0015       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0016       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0017       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0018       COMMON/PYINT1/MINT(400),VINT(400)
0019       COMMON/PYINT4/MWID(500),WIDS(500,5)
0020       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYINT1/,/PYINT4/
0021 C...Local array.
0022       DIMENSION PS(2,6),IJOIN(100)
0023
0024 C...Initialize and reset.
0025       MSTU(24)=0
0026       IF(MSTU(12).NE.12345) CALL PYLIST(0)
0027       MSTU(29)=0
0028       MSTU(31)=MSTU(31)+1
0029       MSTU(1)=0
0030       MSTU(2)=0
0031       MSTU(3)=0
0032       IF(MSTU(17).LE.0) MSTU(90)=0
0033       MCONS=1
0034
0035 C...Sum up momentum, energy and charge for starting entries.
0036       NSAV=N
0037       DO 110 I=1,2
0038         DO 100 J=1,6
0039           PS(I,J)=0D0
0040   100   CONTINUE
0041   110 CONTINUE
0042       DO 130 I=1,N
0043         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 130
0044         DO 120 J=1,4
0045           PS(1,J)=PS(1,J)+P(I,J)
0046   120   CONTINUE
0047         PS(1,6)=PS(1,6)+PYCHGE(K(I,2))
0048   130 CONTINUE
0049       PARU(21)=PS(1,4)
0050
0051 C...Start by all decays of coloured resonances involved in shower.
0052       NORIG=N
0053       DO 140 I=1,NORIG
0054         IF(K(I,1).EQ.3) THEN
0055           KC=PYCOMP(K(I,2))
0056           IF(MWID(KC).NE.0.AND.KCHG(KC,2).NE.0) CALL PYRESD(I)
0057         ENDIF
0058   140 CONTINUE
0059
0060 C...Prepare system for subsequent fragmentation/decay.
0061       CALL PYPREP(0)
0062       IF(MINT(51).NE.0) RETURN
0063
0064 C...Loop through jet fragmentation and particle decays.
0065       MBE=0
0066   150 MBE=MBE+1
0067       IP=0
0068   160 IP=IP+1
0069       KC=0
0070       IF(K(IP,1).GT.0.AND.K(IP,1).LE.10) KC=PYCOMP(K(IP,2))
0071       IF(KC.EQ.0) THEN
0072
0073 C...Deal with any remaining undecayed resonance
0074 C...(normally the task of PYEVNT, so seldom used).
0075       ELSEIF(MWID(KC).NE.0) THEN
0076         IBEG=IP
0077         IF(KCHG(KC,2).NE.0.AND.K(I,1).NE.3) THEN
0078           IBEG=IP+1
0079   170     IBEG=IBEG-1
0080           IF(IBEG.GE.2.AND.K(IBEG,1).EQ.2) GOTO 170
0081           IF(K(IBEG,1).NE.2) IBEG=IBEG+1
0082           IEND=IP-1
0083   180     IEND=IEND+1
0084           IF(IEND.LT.N.AND.K(IEND,1).EQ.2) GOTO 180
0085           IF(IEND.LT.N.AND.KCHG(PYCOMP(K(IEND,2)),2).EQ.0) GOTO 180
0086           NJOIN=0
0087           DO 190 I=IBEG,IEND
0088             IF(KCHG(PYCOMP(K(IEND,2)),2).NE.0) THEN
0089               NJOIN=NJOIN+1
0090               IJOIN(NJOIN)=I
0091             ENDIF
0092   190     CONTINUE
0093         ENDIF
0094         CALL PYRESD(IP)
0095         CALL PYPREP(IBEG)
0096         IF(MINT(51).NE.0) RETURN
0097
0098 C...Particle decay if unstable and allowed. Save long-lived particle
0099 C...decays until second pass after Bose-Einstein effects.
0100       ELSEIF(KCHG(KC,2).EQ.0) THEN
0101         IF(MSTJ(21).GE.1.AND.MDCY(KC,1).GE.1.AND.(MSTJ(51).LE.0.OR.MBE
0102      &  .EQ.2.OR.PMAS(KC,2).GE.PARJ(91).OR.IABS(K(IP,2)).EQ.311))
0103      &  CALL PYDECY(IP)
0104
0105 C...Decay products may develop a shower.
0106         IF(MSTJ(92).GT.0) THEN
0107           IP1=MSTJ(92)
0108           QMAX=SQRT(MAX(0D0,(P(IP1,4)+P(IP1+1,4))**2-(P(IP1,1)+P(IP1+1,
0109      &    1))**2-(P(IP1,2)+P(IP1+1,2))**2-(P(IP1,3)+P(IP1+1,3))**2))
0110           MINT(33)=0
0111           CALL PYSHOW(IP1,IP1+1,QMAX)
0112           CALL PYPREP(IP1)
0113           IF(MINT(51).NE.0) RETURN
0114           MSTJ(92)=0
0115         ELSEIF(MSTJ(92).LT.0) THEN
0116           IP1=-MSTJ(92)
0117           MINT(33)=0
0118           CALL PYSHOW(IP1,-3,P(IP,5))
0119           CALL PYPREP(IP1)
0120           IF(MINT(51).NE.0) RETURN
0121           MSTJ(92)=0
0122         ENDIF
0123
0124 C...Jet fragmentation: string or independent fragmentation.
0125       ELSEIF(K(IP,1).EQ.1.OR.K(IP,1).EQ.2) THEN
0126         MFRAG=MSTJ(1)
0127         IF(MFRAG.GE.1.AND.K(IP,1).EQ.1) MFRAG=2
0128         IF(MSTJ(21).GE.2.AND.K(IP,1).EQ.2.AND.N.GT.IP) THEN
0129           IF(K(IP+1,1).EQ.1.AND.K(IP+1,3).EQ.K(IP,3).AND.
0130      &    K(IP,3).GT.0.AND.K(IP,3).LT.IP) THEN
0131             IF(KCHG(PYCOMP(K(K(IP,3),2)),2).EQ.0) MFRAG=MIN(1,MFRAG)
0132           ENDIF
0133         ENDIF
0134         IF(MFRAG.EQ.1) CALL PYSTRF(IP)
0135         IF(MFRAG.EQ.2) CALL PYINDF(IP)
0136         IF(MFRAG.EQ.2.AND.K(IP,1).EQ.1) MCONS=0
0137         IF(MFRAG.EQ.2.AND.(MSTJ(3).LE.0.OR.MOD(MSTJ(3),5).EQ.0)) MCONS=0
0138       ENDIF
0139
0140 C...Loop back if enough space left in PYJETS and no error abort.
0141       IF(MSTU(24).NE.0.AND.MSTU(21).GE.2) THEN
0142       ELSEIF(IP.LT.N.AND.N.LT.MSTU(4)-20-MSTU(32)) THEN
0143         GOTO 160
0144       ELSEIF(IP.LT.N) THEN
0145         CALL PYERRM(11,'(PYEXEC:) no more memory left in PYJETS')
0146       ENDIF
0147
0148 C...Include simple Bose-Einstein effect parametrization if desired.
0149       IF(MBE.EQ.1.AND.MSTJ(51).GE.1) THEN
0150         CALL PYBOEI(NSAV)
0151         GOTO 150
0152       ENDIF
0153
0154 C...Check that momentum, energy and charge were conserved.
0155       DO 210 I=1,N
0156         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 210
0157         DO 200 J=1,4
0158           PS(2,J)=PS(2,J)+P(I,J)
0159   200   CONTINUE
0160         PS(2,6)=PS(2,6)+PYCHGE(K(I,2))
0161   210 CONTINUE
0162       PDEV=(ABS(PS(2,1)-PS(1,1))+ABS(PS(2,2)-PS(1,2))+ABS(PS(2,3)-
0163      &PS(1,3))+ABS(PS(2,4)-PS(1,4)))/(1D0+ABS(PS(2,4))+ABS(PS(1,4)))
0164       IF(MCONS.EQ.1.AND.PDEV.GT.PARU(11)) CALL PYERRM(15,
0165      &'(PYEXEC:) four-momentum was not conserved')
0166       IF(MCONS.EQ.1.AND.ABS(PS(2,6)-PS(1,6)).GT.0.1D0) CALL PYERRM(15,
0167      &'(PYEXEC:) charge was not conserved')
0168
0169       RETURN
0170       END