hijing/src/lutabu.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUTABU(MTABU)
0005
0006 C...Purpose: to evaluate various properties of an event, with
0007 C...statistics accumulated during the course of the run and
0008 C...printed at the end.
0009       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0010       SAVE /LUJETS/
0011       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0012       SAVE /LUDAT1/
0013       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0014       SAVE /LUDAT2/
0015       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
0016       SAVE /LUDAT3/
0017       DIMENSION KFIS(100,2),NPIS(100,0:10),KFFS(400),NPFS(400,4),
0018      &FEVFM(10,4),FM1FM(3,10,4),FM2FM(3,10,4),FMOMA(4),FMOMS(4),
0019      &FEVEE(50),FE1EC(50),FE2EC(50),FE1EA(25),FE2EA(25),
0020      &KFDM(8),KFDC(200,0:8),NPDC(200)
0021       SAVE NEVIS,NKFIS,KFIS,NPIS,NEVFS,NPRFS,NFIFS,NCHFS,NKFFS,
0022      &KFFS,NPFS,NEVFM,NMUFM,FM1FM,FM2FM,NEVEE,FE1EC,FE2EC,FE1EA,
0023      &FE2EA,NEVDC,NKFDC,NREDC,KFDC,NPDC
0024       CHARACTER CHAU*16,CHIS(2)*12,CHDC(8)*12
0025       DATA NEVIS/0/,NKFIS/0/,NEVFS/0/,NPRFS/0/,NFIFS/0/,NCHFS/0/,
0026      &NKFFS/0/,NEVFM/0/,NMUFM/0/,FM1FM/120*0./,FM2FM/120*0./,
0027      &NEVEE/0/,FE1EC/50*0./,FE2EC/50*0./,FE1EA/25*0./,FE2EA/25*0./,
0028      &NEVDC/0/,NKFDC/0/,NREDC/0/
0029
0030 C...Reset statistics on initial parton state.
0031       IF(MTABU.EQ.10) THEN
0032         NEVIS=0
0033         NKFIS=0
0034
0035 C...Identify and order flavour content of initial state.
0036       ELSEIF(MTABU.EQ.11) THEN
0037         NEVIS=NEVIS+1
0038         KFM1=2*IABS(MSTU(161))
0039         IF(MSTU(161).GT.0) KFM1=KFM1-1
0040         KFM2=2*IABS(MSTU(162))
0041         IF(MSTU(162).GT.0) KFM2=KFM2-1
0042         KFMN=MIN(KFM1,KFM2)
0043         KFMX=MAX(KFM1,KFM2)
0044         DO 100 I=1,NKFIS
0045         IF(KFMN.EQ.KFIS(I,1).AND.KFMX.EQ.KFIS(I,2)) THEN
0046           IKFIS=-I
0047           GOTO 110
0048         ELSEIF(KFMN.LT.KFIS(I,1).OR.(KFMN.EQ.KFIS(I,1).AND.
0049      &  KFMX.LT.KFIS(I,2))) THEN
0050           IKFIS=I
0051           GOTO 110
0052         ENDIF
0053   100   CONTINUE
0054         IKFIS=NKFIS+1
0055   110   IF(IKFIS.LT.0) THEN
0056           IKFIS=-IKFIS
0057         ELSE
0058           IF(NKFIS.GE.100) RETURN
0059           DO 120 I=NKFIS,IKFIS,-1
0060           KFIS(I+1,1)=KFIS(I,1)
0061           KFIS(I+1,2)=KFIS(I,2)
0062           DO 120 J=0,10
0063   120     NPIS(I+1,J)=NPIS(I,J)
0064           NKFIS=NKFIS+1
0065           KFIS(IKFIS,1)=KFMN
0066           KFIS(IKFIS,2)=KFMX
0067           DO 130 J=0,10
0068   130     NPIS(IKFIS,J)=0
0069         ENDIF
0070         NPIS(IKFIS,0)=NPIS(IKFIS,0)+1
0071
0072 C...Count number of partons in initial state.
0073         NP=0
0074         DO 150 I=1,N
0075         IF(K(I,1).LE.0.OR.K(I,1).GT.12) THEN
0076         ELSEIF(IABS(K(I,2)).GT.80.AND.IABS(K(I,2)).LE.100) THEN
0077         ELSEIF(IABS(K(I,2)).GT.100.AND.MOD(IABS(K(I,2))/10,10).NE.0)
0078      &  THEN
0079         ELSE
0080           IM=I
0081   140     IM=K(IM,3)
0082           IF(IM.LE.0.OR.IM.GT.N) THEN
0083             NP=NP+1
0084           ELSEIF(K(IM,1).LE.0.OR.K(IM,1).GT.20) THEN
0085             NP=NP+1
0086           ELSEIF(IABS(K(IM,2)).GT.80.AND.IABS(K(IM,2)).LE.100) THEN
0087           ELSEIF(IABS(K(IM,2)).GT.100.AND.MOD(IABS(K(IM,2))/10,10).NE.0)
0088      &    THEN
0089           ELSE
0090             GOTO 140
0091           ENDIF
0092         ENDIF
0093   150   CONTINUE
0094         NPCO=MAX(NP,1)
0095         IF(NP.GE.6) NPCO=6
0096         IF(NP.GE.8) NPCO=7
0097         IF(NP.GE.11) NPCO=8
0098         IF(NP.GE.16) NPCO=9
0099         IF(NP.GE.26) NPCO=10
0100         NPIS(IKFIS,NPCO)=NPIS(IKFIS,NPCO)+1
0101         MSTU(62)=NP
0102
0103 C...Write statistics on initial parton state.
0104       ELSEIF(MTABU.EQ.12) THEN
0105         FAC=1./MAX(1,NEVIS)
0106         WRITE(MSTU(11),1000) NEVIS
0107         DO 160 I=1,NKFIS
0108         KFMN=KFIS(I,1)
0109         IF(KFMN.EQ.0) KFMN=KFIS(I,2)
0110         KFM1=(KFMN+1)/2
0111         IF(2*KFM1.EQ.KFMN) KFM1=-KFM1
0112         CALL LUNAME(KFM1,CHAU)
0113         CHIS(1)=CHAU(1:12)
0114         IF(CHAU(13:13).NE.' ') CHIS(1)(12:12)='?'
0115         KFMX=KFIS(I,2)
0116         IF(KFIS(I,1).EQ.0) KFMX=0
0117         KFM2=(KFMX+1)/2
0118         IF(2*KFM2.EQ.KFMX) KFM2=-KFM2
0119         CALL LUNAME(KFM2,CHAU)
0120         CHIS(2)=CHAU(1:12)
0121         IF(CHAU(13:13).NE.' ') CHIS(2)(12:12)='?'
0122   160   WRITE(MSTU(11),1100) CHIS(1),CHIS(2),FAC*NPIS(I,0),
0123      &  (NPIS(I,J)/FLOAT(NPIS(I,0)),J=1,10)
0124
0125 C...Copy statistics on initial parton state into /LUJETS/.
0126       ELSEIF(MTABU.EQ.13) THEN
0127         FAC=1./MAX(1,NEVIS)
0128         DO 170 I=1,NKFIS
0129         KFMN=KFIS(I,1)
0130         IF(KFMN.EQ.0) KFMN=KFIS(I,2)
0131         KFM1=(KFMN+1)/2
0132         IF(2*KFM1.EQ.KFMN) KFM1=-KFM1
0133         KFMX=KFIS(I,2)
0134         IF(KFIS(I,1).EQ.0) KFMX=0
0135         KFM2=(KFMX+1)/2
0136         IF(2*KFM2.EQ.KFMX) KFM2=-KFM2
0137         K(I,1)=32
0138         K(I,2)=99
0139         K(I,3)=KFM1
0140         K(I,4)=KFM2
0141         K(I,5)=NPIS(I,0)
0142         DO 170 J=1,5
0143         P(I,J)=FAC*NPIS(I,J)
0144   170   V(I,J)=FAC*NPIS(I,J+5)
0145         N=NKFIS
0146         DO 180 J=1,5
0147         K(N+1,J)=0
0148         P(N+1,J)=0.
0149   180   V(N+1,J)=0.
0150         K(N+1,1)=32
0151         K(N+1,2)=99
0152         K(N+1,5)=NEVIS
0153         MSTU(3)=1
0154
0155 C...Reset statistics on number of particles/partons.
0156       ELSEIF(MTABU.EQ.20) THEN
0157         NEVFS=0
0158         NPRFS=0
0159         NFIFS=0
0160         NCHFS=0
0161         NKFFS=0
0162
0163 C...Identify whether particle/parton is primary or not.
0164       ELSEIF(MTABU.EQ.21) THEN
0165         NEVFS=NEVFS+1
0166         MSTU(62)=0
0167         DO 230 I=1,N
0168         IF(K(I,1).LE.0.OR.K(I,1).GT.20.OR.K(I,1).EQ.13) GOTO 230
0169         MSTU(62)=MSTU(62)+1
0170         KC=LUCOMP(K(I,2))
0171         MPRI=0
0172         IF(K(I,3).LE.0.OR.K(I,3).GT.N) THEN
0173           MPRI=1
0174         ELSEIF(K(K(I,3),1).LE.0.OR.K(K(I,3),1).GT.20) THEN
0175           MPRI=1
0176         ELSEIF(KC.EQ.0) THEN
0177         ELSEIF(K(K(I,3),1).EQ.13) THEN
0178           IM=K(K(I,3),3)
0179           IF(IM.LE.0.OR.IM.GT.N) THEN
0180             MPRI=1
0181           ELSEIF(K(IM,1).LE.0.OR.K(IM,1).GT.20) THEN
0182             MPRI=1
0183           ENDIF
0184         ELSEIF(KCHG(KC,2).EQ.0) THEN
0185           KCM=LUCOMP(K(K(I,3),2))
0186           IF(KCM.NE.0) THEN
0187             IF(KCHG(KCM,2).NE.0) MPRI=1
0188           ENDIF
0189         ENDIF
0190         IF(KC.NE.0.AND.MPRI.EQ.1) THEN
0191           IF(KCHG(KC,2).EQ.0) NPRFS=NPRFS+1
0192         ENDIF
0193         IF(K(I,1).LE.10) THEN
0194           NFIFS=NFIFS+1
0195           IF(LUCHGE(K(I,2)).NE.0) NCHFS=NCHFS+1
0196         ENDIF
0197
0198 C...Fill statistics on number of particles/partons in event.
0199         KFA=IABS(K(I,2))
0200         KFS=3-ISIGN(1,K(I,2))-MPRI
0201         DO 190 IP=1,NKFFS
0202         IF(KFA.EQ.KFFS(IP)) THEN
0203           IKFFS=-IP
0204           GOTO 200
0205         ELSEIF(KFA.LT.KFFS(IP)) THEN
0206           IKFFS=IP
0207           GOTO 200
0208         ENDIF
0209   190   CONTINUE
0210         IKFFS=NKFFS+1
0211   200   IF(IKFFS.LT.0) THEN
0212           IKFFS=-IKFFS
0213         ELSE
0214           IF(NKFFS.GE.400) RETURN
0215           DO 210 IP=NKFFS,IKFFS,-1
0216           KFFS(IP+1)=KFFS(IP)
0217           DO 210 J=1,4
0218   210     NPFS(IP+1,J)=NPFS(IP,J)
0219           NKFFS=NKFFS+1
0220           KFFS(IKFFS)=KFA
0221           DO 220 J=1,4
0222   220     NPFS(IKFFS,J)=0
0223         ENDIF
0224         NPFS(IKFFS,KFS)=NPFS(IKFFS,KFS)+1
0225   230   CONTINUE
0226
0227 C...Write statistics on particle/parton composition of events.
0228       ELSEIF(MTABU.EQ.22) THEN
0229         FAC=1./MAX(1,NEVFS)
0230         WRITE(MSTU(11),1200) NEVFS,FAC*NPRFS,FAC*NFIFS,FAC*NCHFS
0231         DO 240 I=1,NKFFS
0232         CALL LUNAME(KFFS(I),CHAU)
0233         KC=LUCOMP(KFFS(I))
0234         MDCYF=0
0235         IF(KC.NE.0) MDCYF=MDCY(KC,1)
0236   240   WRITE(MSTU(11),1300) KFFS(I),CHAU,MDCYF,(FAC*NPFS(I,J),J=1,4),
0237      &  FAC*(NPFS(I,1)+NPFS(I,2)+NPFS(I,3)+NPFS(I,4))
0238
0239 C...Copy particle/parton composition information into /LUJETS/.
0240       ELSEIF(MTABU.EQ.23) THEN
0241         FAC=1./MAX(1,NEVFS)
0242         DO 260 I=1,NKFFS
0243         K(I,1)=32
0244         K(I,2)=99
0245         K(I,3)=KFFS(I)
0246         K(I,4)=0
0247         K(I,5)=NPFS(I,1)+NPFS(I,2)+NPFS(I,3)+NPFS(I,4)
0248         DO 250 J=1,4
0249         P(I,J)=FAC*NPFS(I,J)
0250   250   V(I,J)=0.
0251         P(I,5)=FAC*K(I,5)
0252   260   V(I,5)=0.
0253         N=NKFFS
0254         DO 270 J=1,5
0255         K(N+1,J)=0
0256         P(N+1,J)=0.
0257   270   V(N+1,J)=0.
0258         K(N+1,1)=32
0259         K(N+1,2)=99
0260         K(N+1,5)=NEVFS
0261         P(N+1,1)=FAC*NPRFS
0262         P(N+1,2)=FAC*NFIFS
0263         P(N+1,3)=FAC*NCHFS
0264         MSTU(3)=1
0265
0266 C...Reset factorial moments statistics.
0267       ELSEIF(MTABU.EQ.30) THEN
0268         NEVFM=0
0269         NMUFM=0
0270         DO 280 IM=1,3
0271         DO 280 IB=1,10
0272         DO 280 IP=1,4
0273         FM1FM(IM,IB,IP)=0.
0274   280   FM2FM(IM,IB,IP)=0.
0275
0276 C...Find particles to include, with (pion,pseudo)rapidity and azimuth.
0277       ELSEIF(MTABU.EQ.31) THEN
0278         NEVFM=NEVFM+1
0279         NLOW=N+MSTU(3)
0280         NUPP=NLOW
0281         DO 360 I=1,N
0282         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 360
0283         IF(MSTU(41).GE.2) THEN
0284           KC=LUCOMP(K(I,2))
0285           IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0286      &    KC.EQ.18) GOTO 360
0287           IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0)
0288      &    GOTO 360
0289         ENDIF
0290         PMR=0.
0291         IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) PMR=ULMASS(211)
0292         IF(MSTU(42).GE.2) PMR=P(I,5)
0293         PR=MAX(1E-20,PMR**2+P(I,1)**2+P(I,2)**2)
0294         YETA=SIGN(LOG(MIN((SQRT(PR+P(I,3)**2)+ABS(P(I,3)))/SQRT(PR),
0295      &  1E20)),P(I,3))
0296         IF(ABS(YETA).GT.PARU(57)) GOTO 360
0297         PHI=ULANGL(P(I,1),P(I,2))
0298         IYETA=512.*(YETA+PARU(57))/(2.*PARU(57))
0299         IYETA=MAX(0,MIN(511,IYETA))
0300         IPHI=512.*(PHI+PARU(1))/PARU(2)
0301         IPHI=MAX(0,MIN(511,IPHI))
0302         IYEP=0
0303         DO 290 IB=0,9
0304   290   IYEP=IYEP+4**IB*(2*MOD(IYETA/2**IB,2)+MOD(IPHI/2**IB,2))
0305
0306 C...Order particles in (pseudo)rapidity and/or azimuth.
0307         IF(NUPP.GT.MSTU(4)-5-MSTU(32)) THEN
0308           CALL LUERRM(11,'(LUTABU:) no more memory left in LUJETS')
0309           RETURN
0310         ENDIF
0311         NUPP=NUPP+1
0312         IF(NUPP.EQ.NLOW+1) THEN
0313           K(NUPP,1)=IYETA
0314           K(NUPP,2)=IPHI
0315           K(NUPP,3)=IYEP
0316         ELSE
0317           DO 300 I1=NUPP-1,NLOW+1,-1
0318           IF(IYETA.GE.K(I1,1)) GOTO 310
0319   300     K(I1+1,1)=K(I1,1)
0320   310     K(I1+1,1)=IYETA
0321           DO 320 I1=NUPP-1,NLOW+1,-1
0322           IF(IPHI.GE.K(I1,2)) GOTO 330
0323   320     K(I1+1,2)=K(I1,2)
0324   330     K(I1+1,2)=IPHI
0325           DO 340 I1=NUPP-1,NLOW+1,-1
0326           IF(IYEP.GE.K(I1,3)) GOTO 350
0327   340     K(I1+1,3)=K(I1,3)
0328   350     K(I1+1,3)=IYEP
0329         ENDIF
0330   360   CONTINUE
0331         K(NUPP+1,1)=2**10
0332         K(NUPP+1,2)=2**10
0333         K(NUPP+1,3)=4**10
0334
0335 C...Calculate sum of factorial moments in event.
0336         DO 400 IM=1,3
0337         DO 370 IB=1,10
0338         DO 370 IP=1,4
0339   370   FEVFM(IB,IP)=0.
0340         DO 380 IB=1,10
0341         IF(IM.LE.2) IBIN=2**(10-IB)
0342         IF(IM.EQ.3) IBIN=4**(10-IB)
0343         IAGR=K(NLOW+1,IM)/IBIN
0344         NAGR=1
0345         DO 380 I=NLOW+2,NUPP+1
0346         ICUT=K(I,IM)/IBIN
0347         IF(ICUT.EQ.IAGR) THEN
0348           NAGR=NAGR+1
0349         ELSE
0350           IF(NAGR.EQ.1) THEN
0351           ELSEIF(NAGR.EQ.2) THEN
0352             FEVFM(IB,1)=FEVFM(IB,1)+2.
0353           ELSEIF(NAGR.EQ.3) THEN
0354             FEVFM(IB,1)=FEVFM(IB,1)+6.
0355             FEVFM(IB,2)=FEVFM(IB,2)+6.
0356           ELSEIF(NAGR.EQ.4) THEN
0357             FEVFM(IB,1)=FEVFM(IB,1)+12.
0358             FEVFM(IB,2)=FEVFM(IB,2)+24.
0359             FEVFM(IB,3)=FEVFM(IB,3)+24.
0360           ELSE
0361             FEVFM(IB,1)=FEVFM(IB,1)+NAGR*(NAGR-1.)
0362             FEVFM(IB,2)=FEVFM(IB,2)+NAGR*(NAGR-1.)*(NAGR-2.)
0363             FEVFM(IB,3)=FEVFM(IB,3)+NAGR*(NAGR-1.)*(NAGR-2.)*(NAGR-3.)
0364             FEVFM(IB,4)=FEVFM(IB,4)+NAGR*(NAGR-1.)*(NAGR-2.)*(NAGR-3.)*
0365      &      (NAGR-4.)
0366           ENDIF
0367           IAGR=ICUT
0368           NAGR=1
0369         ENDIF
0370   380   CONTINUE
0371
0372 C...Add results to total statistics.
0373         DO 390 IB=10,1,-1
0374         DO 390 IP=1,4
0375         IF(FEVFM(1,IP).LT.0.5) THEN
0376           FEVFM(IB,IP)=0.
0377         ELSEIF(IM.LE.2) THEN
0378           FEVFM(IB,IP)=2**((IB-1)*IP)*FEVFM(IB,IP)/FEVFM(1,IP)
0379         ELSE
0380           FEVFM(IB,IP)=4**((IB-1)*IP)*FEVFM(IB,IP)/FEVFM(1,IP)
0381         ENDIF
0382         FM1FM(IM,IB,IP)=FM1FM(IM,IB,IP)+FEVFM(IB,IP)
0383   390   FM2FM(IM,IB,IP)=FM2FM(IM,IB,IP)+FEVFM(IB,IP)**2
0384   400   CONTINUE
0385         NMUFM=NMUFM+(NUPP-NLOW)
0386         MSTU(62)=NUPP-NLOW
0387
0388 C...Write accumulated statistics on factorial moments.
0389       ELSEIF(MTABU.EQ.32) THEN
0390         FAC=1./MAX(1,NEVFM)
0391         IF(MSTU(42).LE.0) WRITE(MSTU(11),1400) NEVFM,'eta'
0392         IF(MSTU(42).EQ.1) WRITE(MSTU(11),1400) NEVFM,'ypi'
0393         IF(MSTU(42).GE.2) WRITE(MSTU(11),1400) NEVFM,'y  '
0394         DO 420 IM=1,3
0395         WRITE(MSTU(11),1500)
0396         DO 420 IB=1,10
0397         BYETA=2.*PARU(57)
0398         IF(IM.NE.2) BYETA=BYETA/2**(IB-1)
0399         BPHI=PARU(2)
0400         IF(IM.NE.1) BPHI=BPHI/2**(IB-1)
0401         IF(IM.LE.2) BNAVE=FAC*NMUFM/FLOAT(2**(IB-1))
0402         IF(IM.EQ.3) BNAVE=FAC*NMUFM/FLOAT(4**(IB-1))
0403         DO 410 IP=1,4
0404         FMOMA(IP)=FAC*FM1FM(IM,IB,IP)
0405   410   FMOMS(IP)=SQRT(MAX(0.,FAC*(FAC*FM2FM(IM,IB,IP)-FMOMA(IP)**2)))
0406   420   WRITE(MSTU(11),1600) BYETA,BPHI,BNAVE,(FMOMA(IP),FMOMS(IP),
0407      &  IP=1,4)
0408
0409 C...Copy statistics on factorial moments into /LUJETS/.
0410       ELSEIF(MTABU.EQ.33) THEN
0411         FAC=1./MAX(1,NEVFM)
0412         DO 430 IM=1,3
0413         DO 430 IB=1,10
0414         I=10*(IM-1)+IB
0415         K(I,1)=32
0416         K(I,2)=99
0417         K(I,3)=1
0418         IF(IM.NE.2) K(I,3)=2**(IB-1)
0419         K(I,4)=1
0420         IF(IM.NE.1) K(I,4)=2**(IB-1)
0421         K(I,5)=0
0422         P(I,1)=2.*PARU(57)/K(I,3)
0423         V(I,1)=PARU(2)/K(I,4)
0424         DO 430 IP=1,4
0425         P(I,IP+1)=FAC*FM1FM(IM,IB,IP)
0426   430   V(I,IP+1)=SQRT(MAX(0.,FAC*(FAC*FM2FM(IM,IB,IP)-P(I,IP+1)**2)))
0427         N=30
0428         DO 440 J=1,5
0429         K(N+1,J)=0
0430         P(N+1,J)=0.
0431   440   V(N+1,J)=0.
0432         K(N+1,1)=32
0433         K(N+1,2)=99
0434         K(N+1,5)=NEVFM
0435         MSTU(3)=1
0436
0437 C...Reset statistics on Energy-Energy Correlation.
0438       ELSEIF(MTABU.EQ.40) THEN
0439         NEVEE=0
0440         DO 450 J=1,25
0441         FE1EC(J)=0.
0442         FE2EC(J)=0.
0443         FE1EC(51-J)=0.
0444         FE2EC(51-J)=0.
0445         FE1EA(J)=0.
0446   450   FE2EA(J)=0.
0447
0448 C...Find particles to include, with proper assumed mass.
0449       ELSEIF(MTABU.EQ.41) THEN
0450         NEVEE=NEVEE+1
0451         NLOW=N+MSTU(3)
0452         NUPP=NLOW
0453         ECM=0.
0454         DO 460 I=1,N
0455         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 460
0456         IF(MSTU(41).GE.2) THEN
0457           KC=LUCOMP(K(I,2))
0458           IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0459      &    KC.EQ.18) GOTO 460
0460           IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0)
0461      &    GOTO 460
0462         ENDIF
0463         PMR=0.
0464         IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) PMR=ULMASS(211)
0465         IF(MSTU(42).GE.2) PMR=P(I,5)
0466         IF(NUPP.GT.MSTU(4)-5-MSTU(32)) THEN
0467           CALL LUERRM(11,'(LUTABU:) no more memory left in LUJETS')
0468           RETURN
0469         ENDIF
0470         NUPP=NUPP+1
0471         P(NUPP,1)=P(I,1)
0472         P(NUPP,2)=P(I,2)
0473         P(NUPP,3)=P(I,3)
0474         P(NUPP,4)=SQRT(PMR**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
0475         P(NUPP,5)=MAX(1E-10,SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2))
0476         ECM=ECM+P(NUPP,4)
0477   460   CONTINUE
0478         IF(NUPP.EQ.NLOW) RETURN
0479
0480 C...Analyze Energy-Energy Correlation in event.
0481         FAC=(2./ECM**2)*50./PARU(1)
0482         DO 470 J=1,50
0483   470   FEVEE(J)=0.
0484         DO 480 I1=NLOW+2,NUPP
0485         DO 480 I2=NLOW+1,I1-1
0486         CTHE=(P(I1,1)*P(I2,1)+P(I1,2)*P(I2,2)+P(I1,3)*P(I2,3))/
0487      &  (P(I1,5)*P(I2,5))
0488         THE=ACOS(MAX(-1.,MIN(1.,CTHE)))
0489         ITHE=MAX(1,MIN(50,1+INT(50.*THE/PARU(1))))
0490   480   FEVEE(ITHE)=FEVEE(ITHE)+FAC*P(I1,4)*P(I2,4)
0491         DO 490 J=1,25
0492         FE1EC(J)=FE1EC(J)+FEVEE(J)
0493         FE2EC(J)=FE2EC(J)+FEVEE(J)**2
0494         FE1EC(51-J)=FE1EC(51-J)+FEVEE(51-J)
0495         FE2EC(51-J)=FE2EC(51-J)+FEVEE(51-J)**2
0496         FE1EA(J)=FE1EA(J)+(FEVEE(51-J)-FEVEE(J))
0497   490   FE2EA(J)=FE2EA(J)+(FEVEE(51-J)-FEVEE(J))**2
0498         MSTU(62)=NUPP-NLOW
0499
0500 C...Write statistics on Energy-Energy Correlation.
0501       ELSEIF(MTABU.EQ.42) THEN
0502         FAC=1./MAX(1,NEVEE)
0503         WRITE(MSTU(11),1700) NEVEE
0504         DO 500 J=1,25
0505         FEEC1=FAC*FE1EC(J)
0506         FEES1=SQRT(MAX(0.,FAC*(FAC*FE2EC(J)-FEEC1**2)))
0507         FEEC2=FAC*FE1EC(51-J)
0508         FEES2=SQRT(MAX(0.,FAC*(FAC*FE2EC(51-J)-FEEC2**2)))
0509         FEECA=FAC*FE1EA(J)
0510         FEESA=SQRT(MAX(0.,FAC*(FAC*FE2EA(J)-FEECA**2)))
0511   500   WRITE(MSTU(11),1800) 3.6*(J-1),3.6*J,FEEC1,FEES1,FEEC2,FEES2,
0512      &  FEECA,FEESA
0513
0514 C...Copy statistics on Energy-Energy Correlation into /LUJETS/.
0515       ELSEIF(MTABU.EQ.43) THEN
0516         FAC=1./MAX(1,NEVEE)
0517         DO 510 I=1,25
0518         K(I,1)=32
0519         K(I,2)=99
0520         K(I,3)=0
0521         K(I,4)=0
0522         K(I,5)=0
0523         P(I,1)=FAC*FE1EC(I)
0524         V(I,1)=SQRT(MAX(0.,FAC*(FAC*FE2EC(I)-P(I,1)**2)))
0525         P(I,2)=FAC*FE1EC(51-I)
0526         V(I,2)=SQRT(MAX(0.,FAC*(FAC*FE2EC(51-I)-P(I,2)**2)))
0527         P(I,3)=FAC*FE1EA(I)
0528         V(I,3)=SQRT(MAX(0.,FAC*(FAC*FE2EA(I)-P(I,3)**2)))
0529         P(I,4)=PARU(1)*(I-1)/50.
0530         P(I,5)=PARU(1)*I/50.
0531         V(I,4)=3.6*(I-1)
0532   510   V(I,5)=3.6*I
0533         N=25
0534         DO 520 J=1,5
0535         K(N+1,J)=0
0536         P(N+1,J)=0.
0537   520   V(N+1,J)=0.
0538         K(N+1,1)=32
0539         K(N+1,2)=99
0540         K(N+1,5)=NEVEE
0541         MSTU(3)=1
0542
0543 C...Reset statistics on decay channels.
0544       ELSEIF(MTABU.EQ.50) THEN
0545         NEVDC=0
0546         NKFDC=0
0547         NREDC=0
0548
0549 C...Identify and order flavour content of final state.
0550       ELSEIF(MTABU.EQ.51) THEN
0551         NEVDC=NEVDC+1
0552         NDS=0
0553         DO 550 I=1,N
0554         IF(K(I,1).LE.0.OR.K(I,1).GE.6) GOTO 550
0555         NDS=NDS+1
0556         IF(NDS.GT.8) THEN
0557           NREDC=NREDC+1
0558           RETURN
0559         ENDIF
0560         KFM=2*IABS(K(I,2))
0561         IF(K(I,2).LT.0) KFM=KFM-1
0562         DO 530 IDS=NDS-1,1,-1
0563         IIN=IDS+1
0564         IF(KFM.LT.KFDM(IDS)) GOTO 540
0565   530   KFDM(IDS+1)=KFDM(IDS)
0566         IIN=1
0567   540   KFDM(IIN)=KFM
0568   550   CONTINUE
0569
0570 C...Find whether old or new final state.
0571         DO 570 IDC=1,NKFDC
0572         IF(NDS.LT.KFDC(IDC,0)) THEN
0573           IKFDC=IDC
0574           GOTO 580
0575         ELSEIF(NDS.EQ.KFDC(IDC,0)) THEN
0576           DO 560 I=1,NDS
0577           IF(KFDM(I).LT.KFDC(IDC,I)) THEN
0578             IKFDC=IDC
0579             GOTO 580
0580           ELSEIF(KFDM(I).GT.KFDC(IDC,I)) THEN
0581             GOTO 570
0582           ENDIF
0583   560     CONTINUE
0584           IKFDC=-IDC
0585           GOTO 580
0586         ENDIF
0587   570   CONTINUE
0588         IKFDC=NKFDC+1
0589   580   IF(IKFDC.LT.0) THEN
0590           IKFDC=-IKFDC
0591         ELSEIF(NKFDC.GE.200) THEN
0592           NREDC=NREDC+1
0593           RETURN
0594         ELSE
0595           DO 590 IDC=NKFDC,IKFDC,-1
0596           NPDC(IDC+1)=NPDC(IDC)
0597           DO 590 I=0,8
0598   590     KFDC(IDC+1,I)=KFDC(IDC,I)
0599           NKFDC=NKFDC+1
0600           KFDC(IKFDC,0)=NDS
0601           DO 600 I=1,NDS
0602   600     KFDC(IKFDC,I)=KFDM(I)
0603           NPDC(IKFDC)=0
0604         ENDIF
0605         NPDC(IKFDC)=NPDC(IKFDC)+1
0606
0607 C...Write statistics on decay channels.
0608       ELSEIF(MTABU.EQ.52) THEN
0609         FAC=1./MAX(1,NEVDC)
0610         WRITE(MSTU(11),1900) NEVDC
0611         DO 620 IDC=1,NKFDC
0612         DO 610 I=1,KFDC(IDC,0)
0613         KFM=KFDC(IDC,I)
0614         KF=(KFM+1)/2
0615         IF(2*KF.NE.KFM) KF=-KF
0616         CALL LUNAME(KF,CHAU)
0617         CHDC(I)=CHAU(1:12)
0618   610   IF(CHAU(13:13).NE.' ') CHDC(I)(12:12)='?'
0619   620   WRITE(MSTU(11),2000) FAC*NPDC(IDC),(CHDC(I),I=1,KFDC(IDC,0))
0620         IF(NREDC.NE.0) WRITE(MSTU(11),2100) FAC*NREDC
0621
0622 C...Copy statistics on decay channels into /LUJETS/.
0623       ELSEIF(MTABU.EQ.53) THEN
0624         FAC=1./MAX(1,NEVDC)
0625         DO 650 IDC=1,NKFDC
0626         K(IDC,1)=32
0627         K(IDC,2)=99
0628         K(IDC,3)=0
0629         K(IDC,4)=0
0630         K(IDC,5)=KFDC(IDC,0)
0631         DO 630 J=1,5
0632         P(IDC,J)=0.
0633   630   V(IDC,J)=0.
0634         DO 640 I=1,KFDC(IDC,0)
0635         KFM=KFDC(IDC,I)
0636         KF=(KFM+1)/2
0637         IF(2*KF.NE.KFM) KF=-KF
0638         IF(I.LE.5) P(IDC,I)=KF
0639   640   IF(I.GE.6) V(IDC,I-5)=KF
0640   650   V(IDC,5)=FAC*NPDC(IDC)
0641         N=NKFDC
0642         DO 660 J=1,5
0643         K(N+1,J)=0
0644         P(N+1,J)=0.
0645   660   V(N+1,J)=0.
0646         K(N+1,1)=32
0647         K(N+1,2)=99
0648         K(N+1,5)=NEVDC
0649         V(N+1,5)=FAC*NREDC
0650         MSTU(3)=1
0651       ENDIF
0652
0653 C...Format statements for output on unit MSTU(11) (default 6).
0654  1000 FORMAT(///20X,'Event statistics - initial state'/
0655      &20X,'based on an analysis of ',I6,' events'//
0656      &3X,'Main flavours after',8X,'Fraction',4X,'Subfractions ',
0657      &'according to fragmenting system multiplicity'/
0658      &4X,'hard interaction',24X,'1',7X,'2',7X,'3',7X,'4',7X,'5',
0659      &6X,'6-7',5X,'8-10',3X,'11-15',3X,'16-25',4X,'>25'/)
0660  1100 FORMAT(3X,A12,1X,A12,F10.5,1X,10F8.4)
0661  1200 FORMAT(///20X,'Event statistics - final state'/
0662      &20X,'based on an analysis of ',I6,' events'//
0663      &5X,'Mean primary multiplicity =',F8.3/
0664      &5X,'Mean final   multiplicity =',F8.3/
0665      &5X,'Mean charged multiplicity =',F8.3//
0666      &5X,'Number of particles produced per event (directly and via ',
0667      &'decays/branchings)'/
0668      &5X,'KF    Particle/jet  MDCY',8X,'Particles',9X,'Antiparticles',
0669      &5X,'Total'/34X,'prim      seco      prim      seco'/)
0670  1300 FORMAT(1X,I6,4X,A16,I2,5(1X,F9.4))
0671  1400 FORMAT(///20X,'Factorial moments analysis of multiplicity'/
0672      &20X,'based on an analysis of ',I6,' events'//
0673      &3X,'delta-',A3,' delta-phi     <n>/bin',10X,'<F2>',18X,'<F3>',
0674      &18X,'<F4>',18X,'<F5>'/35X,4('     value     error  '))
0675  1500 FORMAT(10X)
0676  1600 FORMAT(2X,2F10.4,F12.4,4(F12.4,F10.4))
0677  1700 FORMAT(///20X,'Energy-Energy Correlation and Asymmetry'/
0678      &20X,'based on an analysis of ',I6,' events'//
0679      &2X,'theta range',8X,'EEC(theta)',8X,'EEC(180-theta)',7X,
0680      &'EECA(theta)'/2X,'in degrees ',3('      value    error')/)
0681  1800 FORMAT(2X,F4.1,' - ',F4.1,3(F11.4,F9.4))
0682  1900 FORMAT(///20X,'Decay channel analysis - final state'/
0683      &20X,'based on an analysis of ',I6,' events'//
0684      &2X,'Probability',10X,'Complete final state'/)
0685  2000 FORMAT(2X,F9.5,5X,8(A12,1X))
0686  2100 FORMAT(2X,F9.5,5X,'into other channels (more than 8 particles ',
0687      &'or table overflow)')
0688
0689       RETURN
0690       END