PYTHIA/jetset/lutabu.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE LUTABU(MTABU)
   5
   6 C...Purpose: to evaluate various properties of an event, with
   7 C...statistics accumulated during the course of the run and
   8 C...printed at the end.
   9       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
  10       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  11       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  12       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
  13       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/
  14       DIMENSION KFIS(100,2),NPIS(100,0:10),KFFS(400),NPFS(400,4),
  15      &FEVFM(10,4),FM1FM(3,10,4),FM2FM(3,10,4),FMOMA(4),FMOMS(4),
  16      &FEVEE(50),FE1EC(50),FE2EC(50),FE1EA(25),FE2EA(25),
  17      &KFDM(8),KFDC(200,0:8),NPDC(200)
  18       SAVE NEVIS,NKFIS,KFIS,NPIS,NEVFS,NPRFS,NFIFS,NCHFS,NKFFS,
  19      &KFFS,NPFS,NEVFM,NMUFM,FM1FM,FM2FM,NEVEE,FE1EC,FE2EC,FE1EA,
  20      &FE2EA,NEVDC,NKFDC,NREDC,KFDC,NPDC
  21       CHARACTER CHAU*16,CHIS(2)*12,CHDC(8)*12
  22       DATA NEVIS/0/,NKFIS/0/,NEVFS/0/,NPRFS/0/,NFIFS/0/,NCHFS/0/,
  23      &NKFFS/0/,NEVFM/0/,NMUFM/0/,FM1FM/120*0./,FM2FM/120*0./,
  24      &NEVEE/0/,FE1EC/50*0./,FE2EC/50*0./,FE1EA/25*0./,FE2EA/25*0./,
  25      &NEVDC/0/,NKFDC/0/,NREDC/0/
  26
  27 C...Reset statistics on initial parton state.
  28       IF(MTABU.EQ.10) THEN
  29         NEVIS=0
  30         NKFIS=0
  31
  32 C...Identify and order flavour content of initial state.
  33       ELSEIF(MTABU.EQ.11) THEN
  34         NEVIS=NEVIS+1
  35         KFM1=2*IABS(MSTU(161))
  36         IF(MSTU(161).GT.0) KFM1=KFM1-1
  37         KFM2=2*IABS(MSTU(162))
  38         IF(MSTU(162).GT.0) KFM2=KFM2-1
  39         KFMN=MIN(KFM1,KFM2)
  40         KFMX=MAX(KFM1,KFM2)
  41         DO 100 I=1,NKFIS
  42         IF(KFMN.EQ.KFIS(I,1).AND.KFMX.EQ.KFIS(I,2)) THEN
  43           IKFIS=-I
  44           GOTO 110
  45         ELSEIF(KFMN.LT.KFIS(I,1).OR.(KFMN.EQ.KFIS(I,1).AND.
  46      &  KFMX.LT.KFIS(I,2))) THEN
  47           IKFIS=I
  48           GOTO 110
  49         ENDIF
  50   100   CONTINUE
  51         IKFIS=NKFIS+1
  52   110   IF(IKFIS.LT.0) THEN
  53           IKFIS=-IKFIS
  54         ELSE
  55           IF(NKFIS.GE.100) RETURN
  56           DO 130 I=NKFIS,IKFIS,-1
  57           KFIS(I+1,1)=KFIS(I,1)
  58           KFIS(I+1,2)=KFIS(I,2)
  59           DO 120 J=0,10
  60           NPIS(I+1,J)=NPIS(I,J)
  61   120     CONTINUE
  62   130   CONTINUE
  63           NKFIS=NKFIS+1
  64           KFIS(IKFIS,1)=KFMN
  65           KFIS(IKFIS,2)=KFMX
  66           DO 140 J=0,10
  67           NPIS(IKFIS,J)=0
  68   140     CONTINUE
  69         ENDIF
  70         NPIS(IKFIS,0)=NPIS(IKFIS,0)+1
  71
  72 C...Count number of partons in initial state.
  73         NP=0
  74         DO 160 I=1,N
  75         IF(K(I,1).LE.0.OR.K(I,1).GT.12) THEN
  76         ELSEIF(IABS(K(I,2)).GT.80.AND.IABS(K(I,2)).LE.100) THEN
  77         ELSEIF(IABS(K(I,2)).GT.100.AND.MOD(IABS(K(I,2))/10,10).NE.0)
  78      &  THEN
  79         ELSE
  80           IM=I
  81   150     IM=K(IM,3)
  82           IF(IM.LE.0.OR.IM.GT.N) THEN
  83             NP=NP+1
  84           ELSEIF(K(IM,1).LE.0.OR.K(IM,1).GT.20) THEN
  85             NP=NP+1
  86           ELSEIF(IABS(K(IM,2)).GT.80.AND.IABS(K(IM,2)).LE.100) THEN
  87           ELSEIF(IABS(K(IM,2)).GT.100.AND.MOD(IABS(K(IM,2))/10,10).NE.0)
  88      &    THEN
  89           ELSE
  90             GOTO 150
  91           ENDIF
  92         ENDIF
  93   160   CONTINUE
  94         NPCO=MAX(NP,1)
  95         IF(NP.GE.6) NPCO=6
  96         IF(NP.GE.8) NPCO=7
  97         IF(NP.GE.11) NPCO=8
  98         IF(NP.GE.16) NPCO=9
  99         IF(NP.GE.26) NPCO=10
 100         NPIS(IKFIS,NPCO)=NPIS(IKFIS,NPCO)+1
 101         MSTU(62)=NP
 102
 103 C...Write statistics on initial parton state.
 104       ELSEIF(MTABU.EQ.12) THEN
 105         FAC=1./MAX(1,NEVIS)
 106         WRITE(MSTU(11),5000) NEVIS
 107         DO 170 I=1,NKFIS
 108         KFMN=KFIS(I,1)
 109         IF(KFMN.EQ.0) KFMN=KFIS(I,2)
 110         KFM1=(KFMN+1)/2
 111         IF(2*KFM1.EQ.KFMN) KFM1=-KFM1
 112         CALL LUNAME(KFM1,CHAU)
 113         CHIS(1)=CHAU(1:12)
 114         IF(CHAU(13:13).NE.' ') CHIS(1)(12:12)='?'
 115         KFMX=KFIS(I,2)
 116         IF(KFIS(I,1).EQ.0) KFMX=0
 117         KFM2=(KFMX+1)/2
 118         IF(2*KFM2.EQ.KFMX) KFM2=-KFM2
 119         CALL LUNAME(KFM2,CHAU)
 120         CHIS(2)=CHAU(1:12)
 121         IF(CHAU(13:13).NE.' ') CHIS(2)(12:12)='?'
 122         WRITE(MSTU(11),5100) CHIS(1),CHIS(2),FAC*NPIS(I,0),
 123      &  (NPIS(I,J)/FLOAT(NPIS(I,0)),J=1,10)
 124   170   CONTINUE
 125
 126 C...Copy statistics on initial parton state into /LUJETS/.
 127       ELSEIF(MTABU.EQ.13) THEN
 128         FAC=1./MAX(1,NEVIS)
 129         DO 190 I=1,NKFIS
 130         KFMN=KFIS(I,1)
 131         IF(KFMN.EQ.0) KFMN=KFIS(I,2)
 132         KFM1=(KFMN+1)/2
 133         IF(2*KFM1.EQ.KFMN) KFM1=-KFM1
 134         KFMX=KFIS(I,2)
 135         IF(KFIS(I,1).EQ.0) KFMX=0
 136         KFM2=(KFMX+1)/2
 137         IF(2*KFM2.EQ.KFMX) KFM2=-KFM2
 138         K(I,1)=32
 139         K(I,2)=99
 140         K(I,3)=KFM1
 141         K(I,4)=KFM2
 142         K(I,5)=NPIS(I,0)
 143         DO 180 J=1,5
 144         P(I,J)=FAC*NPIS(I,J)
 145         V(I,J)=FAC*NPIS(I,J+5)
 146   180   CONTINUE
 147   190   CONTINUE
 148         N=NKFIS
 149         DO 200 J=1,5
 150         K(N+1,J)=0
 151         P(N+1,J)=0.
 152         V(N+1,J)=0.
 153   200   CONTINUE
 154         K(N+1,1)=32
 155         K(N+1,2)=99
 156         K(N+1,5)=NEVIS
 157         MSTU(3)=1
 158
 159 C...Reset statistics on number of particles/partons.
 160       ELSEIF(MTABU.EQ.20) THEN
 161         NEVFS=0
 162         NPRFS=0
 163         NFIFS=0
 164         NCHFS=0
 165         NKFFS=0
 166
 167 C...Identify whether particle/parton is primary or not.
 168       ELSEIF(MTABU.EQ.21) THEN
 169         NEVFS=NEVFS+1
 170         MSTU(62)=0
 171         DO 260 I=1,N
 172         IF(K(I,1).LE.0.OR.K(I,1).GT.20.OR.K(I,1).EQ.13) GOTO 260
 173         MSTU(62)=MSTU(62)+1
 174         KC=LUCOMP(K(I,2))
 175         MPRI=0
 176         IF(K(I,3).LE.0.OR.K(I,3).GT.N) THEN
 177           MPRI=1
 178         ELSEIF(K(K(I,3),1).LE.0.OR.K(K(I,3),1).GT.20) THEN
 179           MPRI=1
 180         ELSEIF(K(K(I,3),2).GE.91.AND.K(K(I,3),2).LE.93) THEN
 181           MPRI=1
 182         ELSEIF(KC.EQ.0) THEN
 183         ELSEIF(K(K(I,3),1).EQ.13) THEN
 184           IM=K(K(I,3),3)
 185           IF(IM.LE.0.OR.IM.GT.N) THEN
 186             MPRI=1
 187           ELSEIF(K(IM,1).LE.0.OR.K(IM,1).GT.20) THEN
 188             MPRI=1
 189           ENDIF
 190         ELSEIF(KCHG(KC,2).EQ.0) THEN
 191           KCM=LUCOMP(K(K(I,3),2))
 192           IF(KCM.NE.0) THEN
 193             IF(KCHG(KCM,2).NE.0) MPRI=1
 194           ENDIF
 195         ENDIF
 196         IF(KC.NE.0.AND.MPRI.EQ.1) THEN
 197           IF(KCHG(KC,2).EQ.0) NPRFS=NPRFS+1
 198         ENDIF
 199         IF(K(I,1).LE.10) THEN
 200           NFIFS=NFIFS+1
 201           IF(LUCHGE(K(I,2)).NE.0) NCHFS=NCHFS+1
 202         ENDIF
 203
 204 C...Fill statistics on number of particles/partons in event.
 205         KFA=IABS(K(I,2))
 206         KFS=3-ISIGN(1,K(I,2))-MPRI
 207         DO 210 IP=1,NKFFS
 208         IF(KFA.EQ.KFFS(IP)) THEN
 209           IKFFS=-IP
 210           GOTO 220
 211         ELSEIF(KFA.LT.KFFS(IP)) THEN
 212           IKFFS=IP
 213           GOTO 220
 214         ENDIF
 215   210   CONTINUE
 216         IKFFS=NKFFS+1
 217   220   IF(IKFFS.LT.0) THEN
 218           IKFFS=-IKFFS
 219         ELSE
 220           IF(NKFFS.GE.400) RETURN
 221           DO 240 IP=NKFFS,IKFFS,-1
 222           KFFS(IP+1)=KFFS(IP)
 223           DO 230 J=1,4
 224           NPFS(IP+1,J)=NPFS(IP,J)
 225   230     CONTINUE
 226   240   CONTINUE
 227           NKFFS=NKFFS+1
 228           KFFS(IKFFS)=KFA
 229           DO 250 J=1,4
 230           NPFS(IKFFS,J)=0
 231   250     CONTINUE
 232         ENDIF
 233         NPFS(IKFFS,KFS)=NPFS(IKFFS,KFS)+1
 234   260   CONTINUE
 235
 236 C...Write statistics on particle/parton composition of events.
 237       ELSEIF(MTABU.EQ.22) THEN
 238         FAC=1./MAX(1,NEVFS)
 239         WRITE(MSTU(11),5200) NEVFS,FAC*NPRFS,FAC*NFIFS,FAC*NCHFS
 240         DO 270 I=1,NKFFS
 241         CALL LUNAME(KFFS(I),CHAU)
 242         KC=LUCOMP(KFFS(I))
 243         MDCYF=0
 244         IF(KC.NE.0) MDCYF=MDCY(KC,1)
 245         WRITE(MSTU(11),5300) KFFS(I),CHAU,MDCYF,(FAC*NPFS(I,J),J=1,4),
 246      &  FAC*(NPFS(I,1)+NPFS(I,2)+NPFS(I,3)+NPFS(I,4))
 247   270   CONTINUE
 248
 249 C...Copy particle/parton composition information into /LUJETS/.
 250       ELSEIF(MTABU.EQ.23) THEN
 251         FAC=1./MAX(1,NEVFS)
 252         DO 290 I=1,NKFFS
 253         K(I,1)=32
 254         K(I,2)=99
 255         K(I,3)=KFFS(I)
 256         K(I,4)=0
 257         K(I,5)=NPFS(I,1)+NPFS(I,2)+NPFS(I,3)+NPFS(I,4)
 258         DO 280 J=1,4
 259         P(I,J)=FAC*NPFS(I,J)
 260         V(I,J)=0.
 261   280   CONTINUE
 262         P(I,5)=FAC*K(I,5)
 263         V(I,5)=0.
 264   290   CONTINUE
 265         N=NKFFS
 266         DO 300 J=1,5
 267         K(N+1,J)=0
 268         P(N+1,J)=0.
 269         V(N+1,J)=0.
 270   300   CONTINUE
 271         K(N+1,1)=32
 272         K(N+1,2)=99
 273         K(N+1,5)=NEVFS
 274         P(N+1,1)=FAC*NPRFS
 275         P(N+1,2)=FAC*NFIFS
 276         P(N+1,3)=FAC*NCHFS
 277         MSTU(3)=1
 278
 279 C...Reset factorial moments statistics.
 280       ELSEIF(MTABU.EQ.30) THEN
 281         NEVFM=0
 282         NMUFM=0
 283         DO 330 IM=1,3
 284         DO 320 IB=1,10
 285         DO 310 IP=1,4
 286         FM1FM(IM,IB,IP)=0.
 287         FM2FM(IM,IB,IP)=0.
 288   310   CONTINUE
 289   320   CONTINUE
 290   330   CONTINUE
 291
 292 C...Find particles to include, with (pion,pseudo)rapidity and azimuth.
 293       ELSEIF(MTABU.EQ.31) THEN
 294         NEVFM=NEVFM+1
 295         NLOW=N+MSTU(3)
 296         NUPP=NLOW
 297         DO 410 I=1,N
 298         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 410
 299         IF(MSTU(41).GE.2) THEN
 300           KC=LUCOMP(K(I,2))
 301           IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
 302      &    KC.EQ.18) GOTO 410
 303           IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0)
 304      &    GOTO 410
 305         ENDIF
 306         PMR=0.
 307         IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) PMR=ULMASS(211)
 308         IF(MSTU(42).GE.2) PMR=P(I,5)
 309         PR=MAX(1E-20,PMR**2+P(I,1)**2+P(I,2)**2)
 310         YETA=SIGN(LOG(MIN((SQRT(PR+P(I,3)**2)+ABS(P(I,3)))/SQRT(PR),
 311      &  1E20)),P(I,3))
 312         IF(ABS(YETA).GT.PARU(57)) GOTO 410
 313         PHI=ULANGL(P(I,1),P(I,2))
 314         IYETA=512.*(YETA+PARU(57))/(2.*PARU(57))
 315         IYETA=MAX(0,MIN(511,IYETA))
 316         IPHI=512.*(PHI+PARU(1))/PARU(2)
 317         IPHI=MAX(0,MIN(511,IPHI))
 318         IYEP=0
 319         DO 340 IB=0,9
 320         IYEP=IYEP+4**IB*(2*MOD(IYETA/2**IB,2)+MOD(IPHI/2**IB,2))
 321   340   CONTINUE
 322
 323 C...Order particles in (pseudo)rapidity and/or azimuth.
 324         IF(NUPP.GT.MSTU(4)-5-MSTU(32)) THEN
 325           CALL LUERRM(11,'(LUTABU:) no more memory left in LUJETS')
 326           RETURN
 327         ENDIF
 328         NUPP=NUPP+1
 329         IF(NUPP.EQ.NLOW+1) THEN
 330           K(NUPP,1)=IYETA
 331           K(NUPP,2)=IPHI
 332           K(NUPP,3)=IYEP
 333         ELSE
 334           DO 350 I1=NUPP-1,NLOW+1,-1
 335           IF(IYETA.GE.K(I1,1)) GOTO 360
 336           K(I1+1,1)=K(I1,1)
 337   350     CONTINUE
 338   360     K(I1+1,1)=IYETA
 339           DO 370 I1=NUPP-1,NLOW+1,-1
 340           IF(IPHI.GE.K(I1,2)) GOTO 380
 341           K(I1+1,2)=K(I1,2)
 342   370     CONTINUE
 343   380     K(I1+1,2)=IPHI
 344           DO 390 I1=NUPP-1,NLOW+1,-1
 345           IF(IYEP.GE.K(I1,3)) GOTO 400
 346           K(I1+1,3)=K(I1,3)
 347   390     CONTINUE
 348   400     K(I1+1,3)=IYEP
 349         ENDIF
 350   410   CONTINUE
 351         K(NUPP+1,1)=2**10
 352         K(NUPP+1,2)=2**10
 353         K(NUPP+1,3)=4**10
 354
 355 C...Calculate sum of factorial moments in event.
 356         DO 480 IM=1,3
 357         DO 430 IB=1,10
 358         DO 420 IP=1,4
 359         FEVFM(IB,IP)=0.
 360   420   CONTINUE
 361   430   CONTINUE
 362         DO 450 IB=1,10
 363         IF(IM.LE.2) IBIN=2**(10-IB)
 364         IF(IM.EQ.3) IBIN=4**(10-IB)
 365         IAGR=K(NLOW+1,IM)/IBIN
 366         NAGR=1
 367         DO 440 I=NLOW+2,NUPP+1
 368         ICUT=K(I,IM)/IBIN
 369         IF(ICUT.EQ.IAGR) THEN
 370           NAGR=NAGR+1
 371         ELSE
 372           IF(NAGR.EQ.1) THEN
 373           ELSEIF(NAGR.EQ.2) THEN
 374             FEVFM(IB,1)=FEVFM(IB,1)+2.
 375           ELSEIF(NAGR.EQ.3) THEN
 376             FEVFM(IB,1)=FEVFM(IB,1)+6.
 377             FEVFM(IB,2)=FEVFM(IB,2)+6.
 378           ELSEIF(NAGR.EQ.4) THEN
 379             FEVFM(IB,1)=FEVFM(IB,1)+12.
 380             FEVFM(IB,2)=FEVFM(IB,2)+24.
 381             FEVFM(IB,3)=FEVFM(IB,3)+24.
 382           ELSE
 383             FEVFM(IB,1)=FEVFM(IB,1)+NAGR*(NAGR-1.)
 384             FEVFM(IB,2)=FEVFM(IB,2)+NAGR*(NAGR-1.)*(NAGR-2.)
 385             FEVFM(IB,3)=FEVFM(IB,3)+NAGR*(NAGR-1.)*(NAGR-2.)*(NAGR-3.)
 386             FEVFM(IB,4)=FEVFM(IB,4)+NAGR*(NAGR-1.)*(NAGR-2.)*(NAGR-3.)*
 387      &      (NAGR-4.)
 388           ENDIF
 389           IAGR=ICUT
 390           NAGR=1
 391         ENDIF
 392   440   CONTINUE
 393   450   CONTINUE
 394
 395 C...Add results to total statistics.
 396         DO 470 IB=10,1,-1
 397         DO 460 IP=1,4
 398         IF(FEVFM(1,IP).LT.0.5) THEN
 399           FEVFM(IB,IP)=0.
 400         ELSEIF(IM.LE.2) THEN
 401           FEVFM(IB,IP)=2.**((IB-1)*IP)*FEVFM(IB,IP)/FEVFM(1,IP)
 402         ELSE
 403           FEVFM(IB,IP)=4.**((IB-1)*IP)*FEVFM(IB,IP)/FEVFM(1,IP)
 404         ENDIF
 405         FM1FM(IM,IB,IP)=FM1FM(IM,IB,IP)+FEVFM(IB,IP)
 406         FM2FM(IM,IB,IP)=FM2FM(IM,IB,IP)+FEVFM(IB,IP)**2
 407   460   CONTINUE
 408   470   CONTINUE
 409   480   CONTINUE
 410         NMUFM=NMUFM+(NUPP-NLOW)
 411         MSTU(62)=NUPP-NLOW
 412
 413 C...Write accumulated statistics on factorial moments.
 414       ELSEIF(MTABU.EQ.32) THEN
 415         FAC=1./MAX(1,NEVFM)
 416         IF(MSTU(42).LE.0) WRITE(MSTU(11),5400) NEVFM,'eta'
 417         IF(MSTU(42).EQ.1) WRITE(MSTU(11),5400) NEVFM,'ypi'
 418         IF(MSTU(42).GE.2) WRITE(MSTU(11),5400) NEVFM,'y  '
 419         DO 510 IM=1,3
 420         WRITE(MSTU(11),5500)
 421         DO 500 IB=1,10
 422         BYETA=2.*PARU(57)
 423         IF(IM.NE.2) BYETA=BYETA/2**(IB-1)
 424         BPHI=PARU(2)
 425         IF(IM.NE.1) BPHI=BPHI/2**(IB-1)
 426         IF(IM.LE.2) BNAVE=FAC*NMUFM/FLOAT(2**(IB-1))
 427         IF(IM.EQ.3) BNAVE=FAC*NMUFM/FLOAT(4**(IB-1))
 428         DO 490 IP=1,4
 429         FMOMA(IP)=FAC*FM1FM(IM,IB,IP)
 430         FMOMS(IP)=SQRT(MAX(0.,FAC*(FAC*FM2FM(IM,IB,IP)-FMOMA(IP)**2)))
 431   490   CONTINUE
 432         WRITE(MSTU(11),5600) BYETA,BPHI,BNAVE,(FMOMA(IP),FMOMS(IP),
 433      &  IP=1,4)
 434   500   CONTINUE
 435   510   CONTINUE
 436
 437 C...Copy statistics on factorial moments into /LUJETS/.
 438       ELSEIF(MTABU.EQ.33) THEN
 439         FAC=1./MAX(1,NEVFM)
 440         DO 540 IM=1,3
 441         DO 530 IB=1,10
 442         I=10*(IM-1)+IB
 443         K(I,1)=32
 444         K(I,2)=99
 445         K(I,3)=1
 446         IF(IM.NE.2) K(I,3)=2**(IB-1)
 447         K(I,4)=1
 448         IF(IM.NE.1) K(I,4)=2**(IB-1)
 449         K(I,5)=0
 450         P(I,1)=2.*PARU(57)/K(I,3)
 451         V(I,1)=PARU(2)/K(I,4)
 452         DO 520 IP=1,4
 453         P(I,IP+1)=FAC*FM1FM(IM,IB,IP)
 454         V(I,IP+1)=SQRT(MAX(0.,FAC*(FAC*FM2FM(IM,IB,IP)-P(I,IP+1)**2)))
 455   520   CONTINUE
 456   530   CONTINUE
 457   540   CONTINUE
 458         N=30
 459         DO 550 J=1,5
 460         K(N+1,J)=0
 461         P(N+1,J)=0.
 462         V(N+1,J)=0.
 463   550   CONTINUE
 464         K(N+1,1)=32
 465         K(N+1,2)=99
 466         K(N+1,5)=NEVFM
 467         MSTU(3)=1
 468
 469 C...Reset statistics on Energy-Energy Correlation.
 470       ELSEIF(MTABU.EQ.40) THEN
 471         NEVEE=0
 472         DO 560 J=1,25
 473         FE1EC(J)=0.
 474         FE2EC(J)=0.
 475         FE1EC(51-J)=0.
 476         FE2EC(51-J)=0.
 477         FE1EA(J)=0.
 478         FE2EA(J)=0.
 479   560   CONTINUE
 480
 481 C...Find particles to include, with proper assumed mass.
 482       ELSEIF(MTABU.EQ.41) THEN
 483         NEVEE=NEVEE+1
 484         NLOW=N+MSTU(3)
 485         NUPP=NLOW
 486         ECM=0.
 487         DO 570 I=1,N
 488         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 570
 489         IF(MSTU(41).GE.2) THEN
 490           KC=LUCOMP(K(I,2))
 491           IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
 492      &    KC.EQ.18) GOTO 570
 493           IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0)
 494      &    GOTO 570
 495         ENDIF
 496         PMR=0.
 497         IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) PMR=ULMASS(211)
 498         IF(MSTU(42).GE.2) PMR=P(I,5)
 499         IF(NUPP.GT.MSTU(4)-5-MSTU(32)) THEN
 500           CALL LUERRM(11,'(LUTABU:) no more memory left in LUJETS')
 501           RETURN
 502         ENDIF
 503         NUPP=NUPP+1
 504         P(NUPP,1)=P(I,1)
 505         P(NUPP,2)=P(I,2)
 506         P(NUPP,3)=P(I,3)
 507         P(NUPP,4)=SQRT(PMR**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
 508         P(NUPP,5)=MAX(1E-10,SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2))
 509         ECM=ECM+P(NUPP,4)
 510   570   CONTINUE
 511         IF(NUPP.EQ.NLOW) RETURN
 512
 513 C...Analyze Energy-Energy Correlation in event.
 514         FAC=(2./ECM**2)*50./PARU(1)
 515         DO 580 J=1,50
 516         FEVEE(J)=0.
 517   580   CONTINUE
 518         DO 600 I1=NLOW+2,NUPP
 519         DO 590 I2=NLOW+1,I1-1
 520         CTHE=(P(I1,1)*P(I2,1)+P(I1,2)*P(I2,2)+P(I1,3)*P(I2,3))/
 521      &  (P(I1,5)*P(I2,5))
 522         THE=ACOS(MAX(-1.,MIN(1.,CTHE)))
 523         ITHE=MAX(1,MIN(50,1+INT(50.*THE/PARU(1))))
 524         FEVEE(ITHE)=FEVEE(ITHE)+FAC*P(I1,4)*P(I2,4)
 525   590   CONTINUE
 526   600   CONTINUE
 527         DO 610 J=1,25
 528         FE1EC(J)=FE1EC(J)+FEVEE(J)
 529         FE2EC(J)=FE2EC(J)+FEVEE(J)**2
 530         FE1EC(51-J)=FE1EC(51-J)+FEVEE(51-J)
 531         FE2EC(51-J)=FE2EC(51-J)+FEVEE(51-J)**2
 532         FE1EA(J)=FE1EA(J)+(FEVEE(51-J)-FEVEE(J))
 533         FE2EA(J)=FE2EA(J)+(FEVEE(51-J)-FEVEE(J))**2
 534   610   CONTINUE
 535         MSTU(62)=NUPP-NLOW
 536
 537 C...Write statistics on Energy-Energy Correlation.
 538       ELSEIF(MTABU.EQ.42) THEN
 539         FAC=1./MAX(1,NEVEE)
 540         WRITE(MSTU(11),5700) NEVEE
 541         DO 620 J=1,25
 542         FEEC1=FAC*FE1EC(J)
 543         FEES1=SQRT(MAX(0.,FAC*(FAC*FE2EC(J)-FEEC1**2)))
 544         FEEC2=FAC*FE1EC(51-J)
 545         FEES2=SQRT(MAX(0.,FAC*(FAC*FE2EC(51-J)-FEEC2**2)))
 546         FEECA=FAC*FE1EA(J)
 547         FEESA=SQRT(MAX(0.,FAC*(FAC*FE2EA(J)-FEECA**2)))
 548         WRITE(MSTU(11),5800) 3.6*(J-1),3.6*J,FEEC1,FEES1,FEEC2,FEES2,
 549      &  FEECA,FEESA
 550   620   CONTINUE
 551
 552 C...Copy statistics on Energy-Energy Correlation into /LUJETS/.
 553       ELSEIF(MTABU.EQ.43) THEN
 554         FAC=1./MAX(1,NEVEE)
 555         DO 630 I=1,25
 556         K(I,1)=32
 557         K(I,2)=99
 558         K(I,3)=0
 559         K(I,4)=0
 560         K(I,5)=0
 561         P(I,1)=FAC*FE1EC(I)
 562         V(I,1)=SQRT(MAX(0.,FAC*(FAC*FE2EC(I)-P(I,1)**2)))
 563         P(I,2)=FAC*FE1EC(51-I)
 564         V(I,2)=SQRT(MAX(0.,FAC*(FAC*FE2EC(51-I)-P(I,2)**2)))
 565         P(I,3)=FAC*FE1EA(I)
 566         V(I,3)=SQRT(MAX(0.,FAC*(FAC*FE2EA(I)-P(I,3)**2)))
 567         P(I,4)=PARU(1)*(I-1)/50.
 568         P(I,5)=PARU(1)*I/50.
 569         V(I,4)=3.6*(I-1)
 570         V(I,5)=3.6*I
 571   630   CONTINUE
 572         N=25
 573         DO 640 J=1,5
 574         K(N+1,J)=0
 575         P(N+1,J)=0.
 576         V(N+1,J)=0.
 577   640   CONTINUE
 578         K(N+1,1)=32
 579         K(N+1,2)=99
 580         K(N+1,5)=NEVEE
 581         MSTU(3)=1
 582
 583 C...Reset statistics on decay channels.
 584       ELSEIF(MTABU.EQ.50) THEN
 585         NEVDC=0
 586         NKFDC=0
 587         NREDC=0
 588
 589 C...Identify and order flavour content of final state.
 590       ELSEIF(MTABU.EQ.51) THEN
 591         NEVDC=NEVDC+1
 592         NDS=0
 593         DO 670 I=1,N
 594         IF(K(I,1).LE.0.OR.K(I,1).GE.6) GOTO 670
 595         NDS=NDS+1
 596         IF(NDS.GT.8) THEN
 597           NREDC=NREDC+1
 598           RETURN
 599         ENDIF
 600         KFM=2*IABS(K(I,2))
 601         IF(K(I,2).LT.0) KFM=KFM-1
 602         DO 650 IDS=NDS-1,1,-1
 603         IIN=IDS+1
 604         IF(KFM.LT.KFDM(IDS)) GOTO 660
 605         KFDM(IDS+1)=KFDM(IDS)
 606   650   CONTINUE
 607         IIN=1
 608   660   KFDM(IIN)=KFM
 609   670   CONTINUE
 610
 611 C...Find whether old or new final state.
 612         DO 690 IDC=1,NKFDC
 613         IF(NDS.LT.KFDC(IDC,0)) THEN
 614           IKFDC=IDC
 615           GOTO 700
 616         ELSEIF(NDS.EQ.KFDC(IDC,0)) THEN
 617           DO 680 I=1,NDS
 618           IF(KFDM(I).LT.KFDC(IDC,I)) THEN
 619             IKFDC=IDC
 620             GOTO 700
 621           ELSEIF(KFDM(I).GT.KFDC(IDC,I)) THEN
 622             GOTO 690
 623           ENDIF
 624   680     CONTINUE
 625           IKFDC=-IDC
 626           GOTO 700
 627         ENDIF
 628   690   CONTINUE
 629         IKFDC=NKFDC+1
 630   700   IF(IKFDC.LT.0) THEN
 631           IKFDC=-IKFDC
 632         ELSEIF(NKFDC.GE.200) THEN
 633           NREDC=NREDC+1
 634           RETURN
 635         ELSE
 636           DO 720 IDC=NKFDC,IKFDC,-1
 637           NPDC(IDC+1)=NPDC(IDC)
 638           DO 710 I=0,8
 639           KFDC(IDC+1,I)=KFDC(IDC,I)
 640   710     CONTINUE
 641   720     CONTINUE
 642           NKFDC=NKFDC+1
 643           KFDC(IKFDC,0)=NDS
 644           DO 730 I=1,NDS
 645           KFDC(IKFDC,I)=KFDM(I)
 646   730     CONTINUE
 647           NPDC(IKFDC)=0
 648         ENDIF
 649         NPDC(IKFDC)=NPDC(IKFDC)+1
 650
 651 C...Write statistics on decay channels.
 652       ELSEIF(MTABU.EQ.52) THEN
 653         FAC=1./MAX(1,NEVDC)
 654         WRITE(MSTU(11),5900) NEVDC
 655         DO 750 IDC=1,NKFDC
 656         DO 740 I=1,KFDC(IDC,0)
 657         KFM=KFDC(IDC,I)
 658         KF=(KFM+1)/2
 659         IF(2*KF.NE.KFM) KF=-KF
 660         CALL LUNAME(KF,CHAU)
 661         CHDC(I)=CHAU(1:12)
 662         IF(CHAU(13:13).NE.' ') CHDC(I)(12:12)='?'
 663   740   CONTINUE
 664         WRITE(MSTU(11),6000) FAC*NPDC(IDC),(CHDC(I),I=1,KFDC(IDC,0))
 665   750   CONTINUE
 666         IF(NREDC.NE.0) WRITE(MSTU(11),6100) FAC*NREDC
 667
 668 C...Copy statistics on decay channels into /LUJETS/.
 669       ELSEIF(MTABU.EQ.53) THEN
 670         FAC=1./MAX(1,NEVDC)
 671         DO 780 IDC=1,NKFDC
 672         K(IDC,1)=32
 673         K(IDC,2)=99
 674         K(IDC,3)=0
 675         K(IDC,4)=0
 676         K(IDC,5)=KFDC(IDC,0)
 677         DO 760 J=1,5
 678         P(IDC,J)=0.
 679         V(IDC,J)=0.
 680   760   CONTINUE
 681         DO 770 I=1,KFDC(IDC,0)
 682         KFM=KFDC(IDC,I)
 683         KF=(KFM+1)/2
 684         IF(2*KF.NE.KFM) KF=-KF
 685         IF(I.LE.5) P(IDC,I)=KF
 686         IF(I.GE.6) V(IDC,I-5)=KF
 687   770   CONTINUE
 688         V(IDC,5)=FAC*NPDC(IDC)
 689   780   CONTINUE
 690         N=NKFDC
 691         DO 790 J=1,5
 692         K(N+1,J)=0
 693         P(N+1,J)=0.
 694         V(N+1,J)=0.
 695   790   CONTINUE
 696         K(N+1,1)=32
 697         K(N+1,2)=99
 698         K(N+1,5)=NEVDC
 699         V(N+1,5)=FAC*NREDC
 700         MSTU(3)=1
 701       ENDIF
 702
 703 C...Format statements for output on unit MSTU(11) (default 6).
 704  5000 FORMAT(///20X,'Event statistics - initial state'/
 705      &20X,'based on an analysis of ',I6,' events'//
 706      &3X,'Main flavours after',8X,'Fraction',4X,'Subfractions ',
 707      &'according to fragmenting system multiplicity'/
 708      &4X,'hard interaction',24X,'1',7X,'2',7X,'3',7X,'4',7X,'5',
 709      &6X,'6-7',5X,'8-10',3X,'11-15',3X,'16-25',4X,'>25'/)
 710  5100 FORMAT(3X,A12,1X,A12,F10.5,1X,10F8.4)
 711  5200 FORMAT(///20X,'Event statistics - final state'/
 712      &20X,'based on an analysis of ',I7,' events'//
 713      &5X,'Mean primary multiplicity =',F10.4/
 714      &5X,'Mean final   multiplicity =',F10.4/
 715      &5X,'Mean charged multiplicity =',F10.4//
 716      &5X,'Number of particles produced per event (directly and via ',
 717      &'decays/branchings)'/
 718      &5X,'KF    Particle/jet  MDCY',10X,'Particles',13X,'Antiparticles',
 719      &8X,'Total'/35X,'prim        seco        prim        seco'/)
 720  5300 FORMAT(1X,I6,4X,A16,I2,5(1X,F11.6))
 721  5400 FORMAT(///20X,'Factorial moments analysis of multiplicity'/
 722      &20X,'based on an analysis of ',I6,' events'//
 723      &3X,'delta-',A3,' delta-phi     <n>/bin',10X,'<F2>',18X,'<F3>',
 724      &18X,'<F4>',18X,'<F5>'/35X,4('     value     error  '))
 725  5500 FORMAT(10X)
 726  5600 FORMAT(2X,2F10.4,F12.4,4(F12.4,F10.4))
 727  5700 FORMAT(///20X,'Energy-Energy Correlation and Asymmetry'/
 728      &20X,'based on an analysis of ',I6,' events'//
 729      &2X,'theta range',8X,'EEC(theta)',8X,'EEC(180-theta)',7X,
 730      &'EECA(theta)'/2X,'in degrees ',3('      value    error')/)
 731  5800 FORMAT(2X,F4.1,' - ',F4.1,3(F11.4,F9.4))
 732  5900 FORMAT(///20X,'Decay channel analysis - final state'/
 733      &20X,'based on an analysis of ',I6,' events'//
 734      &2X,'Probability',10X,'Complete final state'/)
 735  6000 FORMAT(2X,F9.5,5X,8(A12,1X))
 736  6100 FORMAT(2X,F9.5,5X,'into other channels (more than 8 particles ',
 737      &'or table overflow)')
 738
 739       RETURN
 740       END