PYTHIA/jetset/luexec.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE LUEXEC
   5
   6 C...Purpose: to administrate the fragmentation and decay chain.
   7       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
   8       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
   9       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  10       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
  11       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/
  12       DIMENSION PS(2,6)
  13
  14 C...Initialize and reset.
  15       MSTU(24)=0
  16       IF(MSTU(12).GE.1) CALL LULIST(0)
  17       MSTU(31)=MSTU(31)+1
  18       MSTU(1)=0
  19       MSTU(2)=0
  20       MSTU(3)=0
  21       IF(MSTU(17).LE.0) MSTU(90)=0
  22       MCONS=1
  23
  24 C...Sum up momentum, energy and charge for starting entries.
  25       NSAV=N
  26       DO 110 I=1,2
  27       DO 100 J=1,6
  28       PS(I,J)=0.
  29   100 CONTINUE
  30   110 CONTINUE
  31       DO 130 I=1,N
  32       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 130
  33       DO 120 J=1,4
  34       PS(1,J)=PS(1,J)+P(I,J)
  35   120 CONTINUE
  36       PS(1,6)=PS(1,6)+LUCHGE(K(I,2))
  37   130 CONTINUE
  38       PARU(21)=PS(1,4)
  39
  40 C...Prepare system for subsequent fragmentation/decay.
  41       CALL LUPREP(0)
  42
  43 C...Loop through jet fragmentation and particle decays.
  44       MBE=0
  45   140 MBE=MBE+1
  46       IP=0
  47   150 IP=IP+1
  48       KC=0
  49       IF(K(IP,1).GT.0.AND.K(IP,1).LE.10) KC=LUCOMP(K(IP,2))
  50       IF(KC.EQ.0) THEN
  51
  52 C...Particle decay if unstable and allowed. Save long-lived particle
  53 C...decays until second pass after Bose-Einstein effects.
  54       ELSEIF(KCHG(KC,2).EQ.0) THEN
  55         IF(MSTJ(21).GE.1.AND.MDCY(KC,1).GE.1.AND.(MSTJ(51).LE.0.OR.MBE
  56      &  .EQ.2.OR.PMAS(KC,2).GE.PARJ(91).OR.IABS(K(IP,2)).EQ.311))
  57      &  CALL LUDECY(IP)
  58
  59 C...Decay products may develop a shower.
  60         IF(MSTJ(92).GT.0) THEN
  61           IP1=MSTJ(92)
  62           QMAX=SQRT(MAX(0.,(P(IP1,4)+P(IP1+1,4))**2-(P(IP1,1)+P(IP1+1,
  63      &    1))**2-(P(IP1,2)+P(IP1+1,2))**2-(P(IP1,3)+P(IP1+1,3))**2))
  64           CALL LUSHOW(IP1,IP1+1,QMAX)
  65           CALL LUPREP(IP1)
  66           MSTJ(92)=0
  67         ELSEIF(MSTJ(92).LT.0) THEN
  68           IP1=-MSTJ(92)
  69           CALL LUSHOW(IP1,-3,P(IP,5))
  70           CALL LUPREP(IP1)
  71           MSTJ(92)=0
  72         ENDIF
  73
  74 C...Jet fragmentation: string or independent fragmentation.
  75       ELSEIF(K(IP,1).EQ.1.OR.K(IP,1).EQ.2) THEN
  76         MFRAG=MSTJ(1)
  77         IF(MFRAG.GE.1.AND.K(IP,1).EQ.1) MFRAG=2
  78         IF(MSTJ(21).GE.2.AND.K(IP,1).EQ.2.AND.N.GT.IP) THEN
  79           IF(K(IP+1,1).EQ.1.AND.K(IP+1,3).EQ.K(IP,3).AND.
  80      &    K(IP,3).GT.0.AND.K(IP,3).LT.IP) THEN
  81             IF(KCHG(LUCOMP(K(K(IP,3),2)),2).EQ.0) MFRAG=MIN(1,MFRAG)
  82           ENDIF
  83         ENDIF
  84         IF(MFRAG.EQ.1) CALL LUSTRF(IP)
  85         IF(MFRAG.EQ.2) CALL LUINDF(IP)
  86         IF(MFRAG.EQ.2.AND.K(IP,1).EQ.1) MCONS=0
  87         IF(MFRAG.EQ.2.AND.(MSTJ(3).LE.0.OR.MOD(MSTJ(3),5).EQ.0)) MCONS=0
  88       ENDIF
  89
  90 C...Loop back if enough space left in LUJETS and no error abort.
  91       IF(MSTU(24).NE.0.AND.MSTU(21).GE.2) THEN
  92       ELSEIF(IP.LT.N.AND.N.LT.MSTU(4)-20-MSTU(32)) THEN
  93         GOTO 150
  94       ELSEIF(IP.LT.N) THEN
  95         CALL LUERRM(11,'(LUEXEC:) no more memory left in LUJETS')
  96       ENDIF
  97
  98 C...Include simple Bose-Einstein effect parametrization if desired.
  99       IF(MBE.EQ.1.AND.MSTJ(51).GE.1) THEN
 100         CALL LUBOEI(NSAV)
 101         GOTO 140
 102       ENDIF
 103
 104 C...Check that momentum, energy and charge were conserved.
 105       DO 170 I=1,N
 106       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 170
 107       DO 160 J=1,4
 108       PS(2,J)=PS(2,J)+P(I,J)
 109   160 CONTINUE
 110       PS(2,6)=PS(2,6)+LUCHGE(K(I,2))
 111   170 CONTINUE
 112       PDEV=(ABS(PS(2,1)-PS(1,1))+ABS(PS(2,2)-PS(1,2))+ABS(PS(2,3)-
 113      &PS(1,3))+ABS(PS(2,4)-PS(1,4)))/(1.+ABS(PS(2,4))+ABS(PS(1,4)))
 114       IF(MCONS.EQ.1.AND.PDEV.GT.PARU(11)) CALL LUERRM(15,
 115      &'(LUEXEC:) four-momentum was not conserved')
 116       IF(MCONS.EQ.1.AND.ABS(PS(2,6)-PS(1,6)).GT.0.1) CALL LUERRM(15,
 117      &'(LUEXEC:) charge was not conserved')
 118
 119       RETURN
 120       END