hijing/src/luexec.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUEXEC
0005
0006 C...Purpose: to administrate the fragmentation and decay chain.
0007       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0008       SAVE /LUJETS/
0009       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0010       SAVE /LUDAT1/
0011       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0012       SAVE /LUDAT2/
0013       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
0014       SAVE /LUDAT3/
0015       DIMENSION PS(2,6)
0016
0017 C...Initialize and reset.
0018       MSTU(24)=0
0019       IF(MSTU(12).GE.1) CALL LULIST(0)
0020       MSTU(31)=MSTU(31)+1
0021       MSTU(1)=0
0022       MSTU(2)=0
0023       MSTU(3)=0
0024       MCONS=1
0025
0026 C...Sum up momentum, energy and charge for starting entries.
0027       NSAV=N
0028       DO 100 I=1,2
0029       DO 100 J=1,6
0030   100 PS(I,J)=0.
0031       DO 120 I=1,N
0032       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 120
0033       DO 110 J=1,4
0034   110 PS(1,J)=PS(1,J)+P(I,J)
0035       PS(1,6)=PS(1,6)+LUCHGE(K(I,2))
0036   120 CONTINUE
0037       PARU(21)=PS(1,4)
0038
0039 C...Prepare system for subsequent fragmentation/decay.
0040       CALL LUPREP(0)
0041
0042 C...Loop through jet fragmentation and particle decays.
0043       MBE=0
0044   130 MBE=MBE+1
0045       IP=0
0046   140 IP=IP+1
0047       KC=0
0048       IF(K(IP,1).GT.0.AND.K(IP,1).LE.10) KC=LUCOMP(K(IP,2))
0049       IF(KC.EQ.0) THEN
0050
0051 C...Particle decay if unstable and allowed. Save long-lived particle
0052 C...decays until second pass after Bose-Einstein effects.
0053       ELSEIF(KCHG(KC,2).EQ.0) THEN
0054         IF(MSTJ(21).GE.1.AND.MDCY(KC,1).GE.1.AND.(MSTJ(51).LE.0.OR.MBE.
0055      &  EQ.2.OR.PMAS(KC,2).GE.PARJ(91).OR.IABS(K(IP,2)).EQ.311))
0056      &  CALL LUDECY(IP)
0057
0058 C...Decay products may develop a shower.
0059         IF(MSTJ(92).GT.0) THEN
0060           IP1=MSTJ(92)
0061           QMAX=SQRT(MAX(0.,(P(IP1,4)+P(IP1+1,4))**2-(P(IP1,1)+P(IP1+1,
0062      &    1))**2-(P(IP1,2)+P(IP1+1,2))**2-(P(IP1,3)+P(IP1+1,3))**2))
0063           CALL LUSHOW(IP1,IP1+1,QMAX)
0064           CALL LUPREP(IP1)
0065           MSTJ(92)=0
0066         ELSEIF(MSTJ(92).LT.0) THEN
0067           IP1=-MSTJ(92)
0068           CALL LUSHOW(IP1,-3,P(IP,5))
0069           CALL LUPREP(IP1)
0070           MSTJ(92)=0
0071         ENDIF
0072
0073 C...Jet fragmentation: string or independent fragmentation.
0074       ELSEIF(K(IP,1).EQ.1.OR.K(IP,1).EQ.2) THEN
0075         MFRAG=MSTJ(1)
0076         IF(MFRAG.GE.1.AND.K(IP,1).EQ.1) MFRAG=2
0077         IF(MSTJ(21).GE.2.AND.K(IP,1).EQ.2.AND.N.GT.IP) THEN
0078           IF(K(IP+1,1).EQ.1.AND.K(IP+1,3).EQ.K(IP,3).AND.
0079      &    K(IP,3).GT.0.AND.K(IP,3).LT.IP) THEN
0080             IF(KCHG(LUCOMP(K(K(IP,3),2)),2).EQ.0) MFRAG=MIN(1,MFRAG)
0081           ENDIF
0082         ENDIF
0083         IF(MFRAG.EQ.1) CALL LUSTRF(IP)
0084         IF(MFRAG.EQ.2) CALL LUINDF(IP)
0085         IF(MFRAG.EQ.2.AND.K(IP,1).EQ.1) MCONS=0
0086         IF(MFRAG.EQ.2.AND.(MSTJ(3).LE.0.OR.MOD(MSTJ(3),5).EQ.0)) MCONS=0
0087       ENDIF
0088
0089 C...Loop back if enough space left in LUJETS and no error abort.
0090       IF(MSTU(24).NE.0.AND.MSTU(21).GE.2) THEN
0091       ELSEIF(IP.LT.N.AND.N.LT.MSTU(4)-20-MSTU(32)) THEN
0092         GOTO 140
0093       ELSEIF(IP.LT.N) THEN
0094         CALL LUERRM(11,'(LUEXEC:) no more memory left in LUJETS')
0095       ENDIF
0096
0097 C...Include simple Bose-Einstein effect parametrization if desired.
0098       IF(MBE.EQ.1.AND.MSTJ(51).GE.1) THEN
0099         CALL LUBOEI(NSAV)
0100         GOTO 130
0101       ENDIF
0102
0103 C...Check that momentum, energy and charge were conserved.
0104       DO 160 I=1,N
0105       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 160
0106       DO 150 J=1,4
0107   150 PS(2,J)=PS(2,J)+P(I,J)
0108       PS(2,6)=PS(2,6)+LUCHGE(K(I,2))
0109   160 CONTINUE
0110       PDEV=(ABS(PS(2,1)-PS(1,1))+ABS(PS(2,2)-PS(1,2))+ABS(PS(2,3)-
0111      &PS(1,3))+ABS(PS(2,4)-PS(1,4)))/(1.+ABS(PS(2,4))+ABS(PS(1,4)))
0112       IF(MCONS.EQ.1.AND.PDEV.GT.PARU(11)) CALL LUERRM(15,
0113      &'(LUEXEC:) four-momentum was not conserved')
0114       IF(MCONS.EQ.1.AND.ABS(PS(2,6)-PS(1,6)).GT.0.1) CALL LUERRM(15,
0115      &'(LUEXEC:) charge was not conserved')
0116
0117       RETURN
0118       END