pytabu.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pytabu.f

 
C*********************************************************************
 
C...PYTABU
C...Evaluates various properties of an event, with statistics
C...accumulated during the course of the run and
C...printed at the end.
 
      SUBROUTINE pytabu(MTABU)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Parameter statement to help give large particle numbers.
      parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
     &kexcit=4000000,kdimen=5000000)
C...Commonblocks.
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
      SAVE /pyjets/,/pydat1/,/pydat2/,/pydat3/
C...Local arrays, character variables, saved variables and data.
      dimension kfis(100,2),npis(100,0:10),kffs(400),npfs(400,4),
     &fevfm(10,4),fm1fm(3,10,4),fm2fm(3,10,4),fmoma(4),fmoms(4),
     &fevee(50),fe1ec(50),fe2ec(50),fe1ea(25),fe2ea(25),
     &kfdm(8),kfdc(200,0:8),npdc(200)
      SAVE nevis,nkfis,kfis,npis,nevfs,nprfs,nfifs,nchfs,nkffs,
     &kffs,npfs,nevfm,nmufm,fm1fm,fm2fm,nevee,fe1ec,fe2ec,fe1ea,
     &fe2ea,nevdc,nkfdc,nredc,kfdc,npdc
      CHARACTER chau*16,chis(2)*12,chdc(8)*12
      DATA nevis/0/,nkfis/0/,nevfs/0/,nprfs/0/,nfifs/0/,nchfs/0/,
     &nkffs/0/,nevfm/0/,nmufm/0/,fm1fm/120*0d0/,fm2fm/120*0d0/,
     &nevee/0/,fe1ec/50*0d0/,fe2ec/50*0d0/,fe1ea/25*0d0/,fe2ea/25*0d0/,
     &nevdc/0/,nkfdc/0/,nredc/0/
 
C...Reset statistics on initial parton state.
      IF(mtabu.EQ.10) THEN
        nevis=0
        nkfis=0
 
C...Identify and order flavour content of initial state.
      ELSEIF(mtabu.EQ.11) THEN
        nevis=nevis+1
        kfm1=2*iabs(mstu(161))
        IF(mstu(161).GT.0) kfm1=kfm1-1
        kfm2=2*iabs(mstu(162))
        IF(mstu(162).GT.0) kfm2=kfm2-1
        kfmn=min(kfm1,kfm2)
        kfmx=max(kfm1,kfm2)
        DO 100 i=1,nkfis
          IF(kfmn.EQ.kfis(i,1).AND.kfmx.EQ.kfis(i,2)) THEN
            ikfis=-i
            goto 110
          ELSEIF(kfmn.LT.kfis(i,1).OR.(kfmn.EQ.kfis(i,1).AND.
     &      kfmx.LT.kfis(i,2))) THEN
            ikfis=i
            goto 110
          ENDIF
  100   CONTINUE
        ikfis=nkfis+1
  110   IF(ikfis.LT.0) THEN
          ikfis=-ikfis
        ELSE
          IF(nkfis.GE.100) RETURN
          DO 130 i=nkfis,ikfis,-1
            kfis(i+1,1)=kfis(i,1)
            kfis(i+1,2)=kfis(i,2)
            DO 120 j=0,10
              npis(i+1,j)=npis(i,j)
  120       CONTINUE
  130     CONTINUE
          nkfis=nkfis+1
          kfis(ikfis,1)=kfmn
          kfis(ikfis,2)=kfmx
          DO 140 j=0,10
            npis(ikfis,j)=0
  140     CONTINUE
        ENDIF
        npis(ikfis,0)=npis(ikfis,0)+1
 
C...Count number of partons in initial state.
        np=0
        DO 160 i=1,n
          IF(k(i,1).LE.0.OR.k(i,1).GT.12) THEN
          ELSEIF(iabs(k(i,2)).GT.80.AND.iabs(k(i,2)).LE.100) THEN
          ELSEIF(iabs(k(i,2)).GT.100.AND.mod(iabs(k(i,2))/10,10).NE.0)
     &      THEN
          ELSE
            im=i
  150       im=k(im,3)
            IF(im.LE.0.OR.im.GT.n) THEN
              np=np+1
            ELSEIF(k(im,1).LE.0.OR.k(im,1).GT.20) THEN
              np=np+1
            ELSEIF(iabs(k(im,2)).GT.80.AND.iabs(k(im,2)).LE.100) THEN
            ELSEIF(iabs(k(im,2)).GT.100.AND.mod(iabs(k(im,2))/10,10)
     &        .NE.0) THEN
            ELSE
              goto 150
            ENDIF
          ENDIF
  160   CONTINUE
        npco=max(np,1)
        IF(np.GE.6) npco=6
        IF(np.GE.8) npco=7
        IF(np.GE.11) npco=8
        IF(np.GE.16) npco=9
        IF(np.GE.26) npco=10
        npis(ikfis,npco)=npis(ikfis,npco)+1
        mstu(62)=np
 
C...Write statistics on initial parton state.
      ELSEIF(mtabu.EQ.12) THEN
        fac=1d0/max(1,nevis)
        WRITE(mstu(11),5000) nevis
        DO 170 i=1,nkfis
          kfmn=kfis(i,1)
          IF(kfmn.EQ.0) kfmn=kfis(i,2)
          kfm1=(kfmn+1)/2
          IF(2*kfm1.EQ.kfmn) kfm1=-kfm1
          CALL pyname(kfm1,chau)
          chis(1)=chau(1:12)
          IF(chau(13:13).NE.' ') chis(1)(12:12)='?'
          kfmx=kfis(i,2)
          IF(kfis(i,1).EQ.0) kfmx=0
          kfm2=(kfmx+1)/2
          IF(2*kfm2.EQ.kfmx) kfm2=-kfm2
          CALL pyname(kfm2,chau)
          chis(2)=chau(1:12)
          IF(chau(13:13).NE.' ') chis(2)(12:12)='?'
          WRITE(mstu(11),5100) chis(1),chis(2),fac*npis(i,0),
     &    (npis(i,j)/dble(npis(i,0)),j=1,10)
  170   CONTINUE
 
C...Copy statistics on initial parton state into /PYJETS/.
      ELSEIF(mtabu.EQ.13) THEN
        fac=1d0/max(1,nevis)
        DO 190 i=1,nkfis
          kfmn=kfis(i,1)
          IF(kfmn.EQ.0) kfmn=kfis(i,2)
          kfm1=(kfmn+1)/2
          IF(2*kfm1.EQ.kfmn) kfm1=-kfm1
          kfmx=kfis(i,2)
          IF(kfis(i,1).EQ.0) kfmx=0
          kfm2=(kfmx+1)/2
          IF(2*kfm2.EQ.kfmx) kfm2=-kfm2
          k(i,1)=32
          k(i,2)=99
          k(i,3)=kfm1
          k(i,4)=kfm2
          k(i,5)=npis(i,0)
          DO 180 j=1,5
            p(i,j)=fac*npis(i,j)
            v(i,j)=fac*npis(i,j+5)
  180     CONTINUE
  190   CONTINUE
        n=nkfis
        DO 200 j=1,5
          k(n+1,j)=0
          p(n+1,j)=0d0
          v(n+1,j)=0d0
  200   CONTINUE
        k(n+1,1)=32
        k(n+1,2)=99
        k(n+1,5)=nevis
        mstu(3)=1
 
C...Reset statistics on number of particles/partons.
      ELSEIF(mtabu.EQ.20) THEN
        nevfs=0
        nprfs=0
        nfifs=0
        nchfs=0
        nkffs=0
 
C...Identify whether particle/parton is primary or not.
      ELSEIF(mtabu.EQ.21) THEN
        nevfs=nevfs+1
        mstu(62)=0
        DO 260 i=1,n
          IF(k(i,1).LE.0.OR.k(i,1).GT.20.OR.k(i,1).EQ.13) goto 260
          mstu(62)=mstu(62)+1
          kc=pycomp(k(i,2))
          mpri=0
          IF(k(i,3).LE.0.OR.k(i,3).GT.n) THEN
            mpri=1
          ELSEIF(k(k(i,3),1).LE.0.OR.k(k(i,3),1).GT.20) THEN
            mpri=1
          ELSEIF(k(k(i,3),2).GE.91.AND.k(k(i,3),2).LE.93) THEN
            mpri=1
          ELSEIF(kc.EQ.0) THEN
          ELSEIF(k(k(i,3),1).EQ.13) THEN
            im=k(k(i,3),3)
            IF(im.LE.0.OR.im.GT.n) THEN
              mpri=1
            ELSEIF(k(im,1).LE.0.OR.k(im,1).GT.20) THEN
              mpri=1
            ENDIF
          ELSEIF(kchg(kc,2).EQ.0) THEN
            kcm=pycomp(k(k(i,3),2))
            IF(kcm.NE.0) THEN
              IF(kchg(kcm,2).NE.0) mpri=1
            ENDIF
          ENDIF
          IF(kc.NE.0.AND.mpri.EQ.1) THEN
            IF(kchg(kc,2).EQ.0) nprfs=nprfs+1
          ENDIF
          IF(k(i,1).LE.10) THEN
            nfifs=nfifs+1
            IF(pychge(k(i,2)).NE.0) nchfs=nchfs+1
          ENDIF
 
C...Fill statistics on number of particles/partons in event.
          kfa=iabs(k(i,2))
          kfs=3-isign(1,k(i,2))-mpri
          DO 210 ip=1,nkffs
            IF(kfa.EQ.kffs(ip)) THEN
              ikffs=-ip
              goto 220
            ELSEIF(kfa.LT.kffs(ip)) THEN
              ikffs=ip
              goto 220
            ENDIF
  210     CONTINUE
          ikffs=nkffs+1
  220     IF(ikffs.LT.0) THEN
            ikffs=-ikffs
          ELSE
            IF(nkffs.GE.400) RETURN
            DO 240 ip=nkffs,ikffs,-1
              kffs(ip+1)=kffs(ip)
              DO 230 j=1,4
                npfs(ip+1,j)=npfs(ip,j)
  230         CONTINUE
  240       CONTINUE
            nkffs=nkffs+1
            kffs(ikffs)=kfa
            DO 250 j=1,4
              npfs(ikffs,j)=0
  250       CONTINUE
          ENDIF
          npfs(ikffs,kfs)=npfs(ikffs,kfs)+1
  260   CONTINUE
 
C...Write statistics on particle/parton composition of events.
      ELSEIF(mtabu.EQ.22) THEN
        fac=1d0/max(1,nevfs)
        WRITE(mstu(11),5200) nevfs,fac*nprfs,fac*nfifs,fac*nchfs
        DO 270 i=1,nkffs
          CALL pyname(kffs(i),chau)
          kc=pycomp(kffs(i))
          mdcyf=0
          IF(kc.NE.0) mdcyf=mdcy(kc,1)
          WRITE(mstu(11),5300) kffs(i),chau,mdcyf,(fac*npfs(i,j),j=1,4),
     &    fac*(npfs(i,1)+npfs(i,2)+npfs(i,3)+npfs(i,4))
  270   CONTINUE
 
C...Copy particle/parton composition information into /PYJETS/.
      ELSEIF(mtabu.EQ.23) THEN
        fac=1d0/max(1,nevfs)
        DO 290 i=1,nkffs
          k(i,1)=32
          k(i,2)=99
          k(i,3)=kffs(i)
          k(i,4)=0
          k(i,5)=npfs(i,1)+npfs(i,2)+npfs(i,3)+npfs(i,4)
          DO 280 j=1,4
            p(i,j)=fac*npfs(i,j)
            v(i,j)=0d0
  280     CONTINUE
          p(i,5)=fac*k(i,5)
          v(i,5)=0d0
  290   CONTINUE
        n=nkffs
        DO 300 j=1,5
          k(n+1,j)=0
          p(n+1,j)=0d0
          v(n+1,j)=0d0
  300   CONTINUE
        k(n+1,1)=32
        k(n+1,2)=99
        k(n+1,5)=nevfs
        p(n+1,1)=fac*nprfs
        p(n+1,2)=fac*nfifs
        p(n+1,3)=fac*nchfs
        mstu(3)=1
 
C...Reset factorial moments statistics.
      ELSEIF(mtabu.EQ.30) THEN
        nevfm=0
        nmufm=0
        DO 330 im=1,3
          DO 320 ib=1,10
            DO 310 ip=1,4
              fm1fm(im,ib,ip)=0d0
              fm2fm(im,ib,ip)=0d0
  310       CONTINUE
  320     CONTINUE
  330   CONTINUE
 
C...Find particles to include, with (pion,pseudo)rapidity and azimuth.
      ELSEIF(mtabu.EQ.31) THEN
        nevfm=nevfm+1
        nlow=n+mstu(3)
        nupp=nlow
        DO 410 i=1,n
          IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 410
          IF(mstu(41).GE.2) THEN
            kc=pycomp(k(i,2))
            IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
     &      kc.EQ.18.OR.k(i,2).EQ.ksusy1+22.OR.k(i,2).EQ.39.OR.
     &      k(i,2).EQ.ksusy1+39) goto 410
            IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.
     &      pychge(k(i,2)).EQ.0) goto 410
          ENDIF
          pmr=0d0
          IF(mstu(42).EQ.1.AND.k(i,2).NE.22) pmr=pymass(211)
          IF(mstu(42).GE.2) pmr=p(i,5)
          pr=max(1d-20,pmr**2+p(i,1)**2+p(i,2)**2)
          yeta=sign(log(min((sqrt(pr+p(i,3)**2)+abs(p(i,3)))/sqrt(pr),
     &    1d20)),p(i,3))
          IF(abs(yeta).GT.paru(57)) goto 410
          phi=pyangl(p(i,1),p(i,2))
          iyeta=512d0*(yeta+paru(57))/(2d0*paru(57))
          iyeta=max(0,min(511,iyeta))
          iphi=512d0*(phi+paru(1))/paru(2)
          iphi=max(0,min(511,iphi))
          iyep=0
          DO 340 ib=0,9
            iyep=iyep+4**ib*(2*mod(iyeta/2**ib,2)+mod(iphi/2**ib,2))
  340     CONTINUE
 
C...Order particles in (pseudo)rapidity and/or azimuth.
          IF(nupp.GT.mstu(4)-5-mstu(32)) THEN
            CALL pyerrm(11,'(PYTABU:) no more memory left in PYJETS')
            RETURN
          ENDIF
          nupp=nupp+1
          IF(nupp.EQ.nlow+1) THEN
            k(nupp,1)=iyeta
            k(nupp,2)=iphi
            k(nupp,3)=iyep
          ELSE
            DO 350 i1=nupp-1,nlow+1,-1
              IF(iyeta.GE.k(i1,1)) goto 360
              k(i1+1,1)=k(i1,1)
  350       CONTINUE
  360       k(i1+1,1)=iyeta
            DO 370 i1=nupp-1,nlow+1,-1
              IF(iphi.GE.k(i1,2)) goto 380
              k(i1+1,2)=k(i1,2)
  370       CONTINUE
  380       k(i1+1,2)=iphi
            DO 390 i1=nupp-1,nlow+1,-1
              IF(iyep.GE.k(i1,3)) goto 400
              k(i1+1,3)=k(i1,3)
  390       CONTINUE
  400       k(i1+1,3)=iyep
          ENDIF
  410   CONTINUE
        k(nupp+1,1)=2**10
        k(nupp+1,2)=2**10
        k(nupp+1,3)=4**10
 
C...Calculate sum of factorial moments in event.
        DO 480 im=1,3
          DO 430 ib=1,10
            DO 420 ip=1,4
              fevfm(ib,ip)=0d0
  420       CONTINUE
  430     CONTINUE
          DO 450 ib=1,10
            IF(im.LE.2) ibin=2**(10-ib)
            IF(im.EQ.3) ibin=4**(10-ib)
            iagr=k(nlow+1,im)/ibin
            nagr=1
            DO 440 i=nlow+2,nupp+1
              icut=k(i,im)/ibin
              IF(icut.EQ.iagr) THEN
                nagr=nagr+1
              ELSE
                IF(nagr.EQ.1) THEN
                ELSEIF(nagr.EQ.2) THEN
                  fevfm(ib,1)=fevfm(ib,1)+2d0
                ELSEIF(nagr.EQ.3) THEN
                  fevfm(ib,1)=fevfm(ib,1)+6d0
                  fevfm(ib,2)=fevfm(ib,2)+6d0
                ELSEIF(nagr.EQ.4) THEN
                  fevfm(ib,1)=fevfm(ib,1)+12d0
                  fevfm(ib,2)=fevfm(ib,2)+24d0
                  fevfm(ib,3)=fevfm(ib,3)+24d0
                ELSE
                  fevfm(ib,1)=fevfm(ib,1)+nagr*(nagr-1d0)
                  fevfm(ib,2)=fevfm(ib,2)+nagr*(nagr-1d0)*(nagr-2d0)
                  fevfm(ib,3)=fevfm(ib,3)+nagr*(nagr-1d0)*(nagr-2d0)*
     &            (nagr-3d0)
                  fevfm(ib,4)=fevfm(ib,4)+nagr*(nagr-1d0)*(nagr-2d0)*
     &            (nagr-3d0)*(nagr-4d0)
                ENDIF
                iagr=icut
                nagr=1
              ENDIF
  440       CONTINUE
  450     CONTINUE
 
C...Add results to total statistics.
          DO 470 ib=10,1,-1
            DO 460 ip=1,4
              IF(fevfm(1,ip).LT.0.5d0) THEN
                fevfm(ib,ip)=0d0
              ELSEIF(im.LE.2) THEN
                fevfm(ib,ip)=2d0**((ib-1)*ip)*fevfm(ib,ip)/fevfm(1,ip)
              ELSE
                fevfm(ib,ip)=4d0**((ib-1)*ip)*fevfm(ib,ip)/fevfm(1,ip)
              ENDIF
              fm1fm(im,ib,ip)=fm1fm(im,ib,ip)+fevfm(ib,ip)
              fm2fm(im,ib,ip)=fm2fm(im,ib,ip)+fevfm(ib,ip)**2
  460       CONTINUE
  470     CONTINUE
  480   CONTINUE
        nmufm=nmufm+(nupp-nlow)
        mstu(62)=nupp-nlow
 
C...Write accumulated statistics on factorial moments.
      ELSEIF(mtabu.EQ.32) THEN
        fac=1d0/max(1,nevfm)
        IF(mstu(42).LE.0) WRITE(mstu(11),5400) nevfm,'eta'
        IF(mstu(42).EQ.1) WRITE(mstu(11),5400) nevfm,'ypi'
        IF(mstu(42).GE.2) WRITE(mstu(11),5400) nevfm,'y  '
        DO 510 im=1,3
          WRITE(mstu(11),5500)
          DO 500 ib=1,10
            byeta=2d0*paru(57)
            IF(im.NE.2) byeta=byeta/2**(ib-1)
            bphi=paru(2)
            IF(im.NE.1) bphi=bphi/2**(ib-1)
            IF(im.LE.2) bnave=fac*nmufm/dble(2**(ib-1))
            IF(im.EQ.3) bnave=fac*nmufm/dble(4**(ib-1))
            DO 490 ip=1,4
              fmoma(ip)=fac*fm1fm(im,ib,ip)
              fmoms(ip)=sqrt(max(0d0,fac*(fac*fm2fm(im,ib,ip)-
     &        fmoma(ip)**2)))
  490       CONTINUE
            WRITE(mstu(11),5600) byeta,bphi,bnave,(fmoma(ip),fmoms(ip),
     &      ip=1,4)
  500     CONTINUE
  510   CONTINUE
 
C...Copy statistics on factorial moments into /PYJETS/.
      ELSEIF(mtabu.EQ.33) THEN
        fac=1d0/max(1,nevfm)
        DO 540 im=1,3
          DO 530 ib=1,10
            i=10*(im-1)+ib
            k(i,1)=32
            k(i,2)=99
            k(i,3)=1
            IF(im.NE.2) k(i,3)=2**(ib-1)
            k(i,4)=1
            IF(im.NE.1) k(i,4)=2**(ib-1)
            k(i,5)=0
            p(i,1)=2d0*paru(57)/k(i,3)
            v(i,1)=paru(2)/k(i,4)
            DO 520 ip=1,4
              p(i,ip+1)=fac*fm1fm(im,ib,ip)
              v(i,ip+1)=sqrt(max(0d0,fac*(fac*fm2fm(im,ib,ip)-
     &        p(i,ip+1)**2)))
  520       CONTINUE
  530     CONTINUE
  540   CONTINUE
        n=30
        DO 550 j=1,5
          k(n+1,j)=0
          p(n+1,j)=0d0
          v(n+1,j)=0d0
  550   CONTINUE
        k(n+1,1)=32
        k(n+1,2)=99
        k(n+1,5)=nevfm
        mstu(3)=1
 
C...Reset statistics on Energy-Energy Correlation.
      ELSEIF(mtabu.EQ.40) THEN
        nevee=0
        DO 560 j=1,25
          fe1ec(j)=0d0
          fe2ec(j)=0d0
          fe1ec(51-j)=0d0
          fe2ec(51-j)=0d0
          fe1ea(j)=0d0
          fe2ea(j)=0d0
  560   CONTINUE
 
C...Find particles to include, with proper assumed mass.
      ELSEIF(mtabu.EQ.41) THEN
        nevee=nevee+1
        nlow=n+mstu(3)
        nupp=nlow
        ecm=0d0
        DO 570 i=1,n
          IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 570
          IF(mstu(41).GE.2) THEN
            kc=pycomp(k(i,2))
            IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
     &      kc.EQ.18.OR.k(i,2).EQ.ksusy1+22.OR.k(i,2).EQ.39.OR.
     &      k(i,2).EQ.ksusy1+39) goto 570
            IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.
     &      pychge(k(i,2)).EQ.0) goto 570
          ENDIF
          pmr=0d0
          IF(mstu(42).EQ.1.AND.k(i,2).NE.22) pmr=pymass(211)
          IF(mstu(42).GE.2) pmr=p(i,5)
          IF(nupp.GT.mstu(4)-5-mstu(32)) THEN
            CALL pyerrm(11,'(PYTABU:) no more memory left in PYJETS')
            RETURN
          ENDIF
          nupp=nupp+1
          p(nupp,1)=p(i,1)
          p(nupp,2)=p(i,2)
          p(nupp,3)=p(i,3)
          p(nupp,4)=sqrt(pmr**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)
          p(nupp,5)=max(1d-10,sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2))
          ecm=ecm+p(nupp,4)
  570   CONTINUE
        IF(nupp.EQ.nlow) RETURN
 
C...Analyze Energy-Energy Correlation in event.
        fac=(2d0/ecm**2)*50d0/paru(1)
        DO 580 j=1,50
          fevee(j)=0d0
  580   CONTINUE
        DO 600 i1=nlow+2,nupp
          DO 590 i2=nlow+1,i1-1
            cthe=(p(i1,1)*p(i2,1)+p(i1,2)*p(i2,2)+p(i1,3)*p(i2,3))/
     &      (p(i1,5)*p(i2,5))
            the=acos(max(-1d0,min(1d0,cthe)))
            ithe=max(1,min(50,1+int(50d0*the/paru(1))))
            fevee(ithe)=fevee(ithe)+fac*p(i1,4)*p(i2,4)
  590     CONTINUE
  600   CONTINUE
        DO 610 j=1,25
          fe1ec(j)=fe1ec(j)+fevee(j)
          fe2ec(j)=fe2ec(j)+fevee(j)**2
          fe1ec(51-j)=fe1ec(51-j)+fevee(51-j)
          fe2ec(51-j)=fe2ec(51-j)+fevee(51-j)**2
          fe1ea(j)=fe1ea(j)+(fevee(51-j)-fevee(j))
          fe2ea(j)=fe2ea(j)+(fevee(51-j)-fevee(j))**2
  610   CONTINUE
        mstu(62)=nupp-nlow
 
C...Write statistics on Energy-Energy Correlation.
      ELSEIF(mtabu.EQ.42) THEN
        fac=1d0/max(1,nevee)
        WRITE(mstu(11),5700) nevee
        DO 620 j=1,25
          feec1=fac*fe1ec(j)
          fees1=sqrt(max(0d0,fac*(fac*fe2ec(j)-feec1**2)))
          feec2=fac*fe1ec(51-j)
          fees2=sqrt(max(0d0,fac*(fac*fe2ec(51-j)-feec2**2)))
          feeca=fac*fe1ea(j)
          feesa=sqrt(max(0d0,fac*(fac*fe2ea(j)-feeca**2)))
          WRITE(mstu(11),5800) 3.6d0*(j-1),3.6d0*j,feec1,fees1,
     &    feec2,fees2,feeca,feesa
  620   CONTINUE
 
C...Copy statistics on Energy-Energy Correlation into /PYJETS/.
      ELSEIF(mtabu.EQ.43) THEN
        fac=1d0/max(1,nevee)
        DO 630 i=1,25
          k(i,1)=32
          k(i,2)=99
          k(i,3)=0
          k(i,4)=0
          k(i,5)=0
          p(i,1)=fac*fe1ec(i)
          v(i,1)=sqrt(max(0d0,fac*(fac*fe2ec(i)-p(i,1)**2)))
          p(i,2)=fac*fe1ec(51-i)
          v(i,2)=sqrt(max(0d0,fac*(fac*fe2ec(51-i)-p(i,2)**2)))
          p(i,3)=fac*fe1ea(i)
          v(i,3)=sqrt(max(0d0,fac*(fac*fe2ea(i)-p(i,3)**2)))
          p(i,4)=paru(1)*(i-1)/50d0
          p(i,5)=paru(1)*i/50d0
          v(i,4)=3.6d0*(i-1)
          v(i,5)=3.6d0*i
  630   CONTINUE
        n=25
        DO 640 j=1,5
          k(n+1,j)=0
          p(n+1,j)=0d0
          v(n+1,j)=0d0
  640   CONTINUE
        k(n+1,1)=32
        k(n+1,2)=99
        k(n+1,5)=nevee
        mstu(3)=1
 
C...Reset statistics on decay channels.
      ELSEIF(mtabu.EQ.50) THEN
        nevdc=0
        nkfdc=0
        nredc=0
 
C...Identify and order flavour content of final state.
      ELSEIF(mtabu.EQ.51) THEN
        nevdc=nevdc+1
        nds=0
        DO 670 i=1,n
          IF(k(i,1).LE.0.OR.k(i,1).GE.6) goto 670
          nds=nds+1
          IF(nds.GT.8) THEN
            nredc=nredc+1
            RETURN
          ENDIF
          kfm=2*iabs(k(i,2))
          IF(k(i,2).LT.0) kfm=kfm-1
          DO 650 ids=nds-1,1,-1
            iin=ids+1
            IF(kfm.LT.kfdm(ids)) goto 660
            kfdm(ids+1)=kfdm(ids)
  650     CONTINUE
          iin=1
  660     kfdm(iin)=kfm
  670   CONTINUE
 
C...Find whether old or new final state.
        DO 690 idc=1,nkfdc
          IF(nds.LT.kfdc(idc,0)) THEN
            ikfdc=idc
            goto 700
          ELSEIF(nds.EQ.kfdc(idc,0)) THEN
            DO 680 i=1,nds
              IF(kfdm(i).LT.kfdc(idc,i)) THEN
                ikfdc=idc
                goto 700
              ELSEIF(kfdm(i).GT.kfdc(idc,i)) THEN
                goto 690
              ENDIF
  680       CONTINUE
            ikfdc=-idc
            goto 700
          ENDIF
  690   CONTINUE
        ikfdc=nkfdc+1
  700   IF(ikfdc.LT.0) THEN
          ikfdc=-ikfdc
        ELSEIF(nkfdc.GE.200) THEN
          nredc=nredc+1
          RETURN
        ELSE
          DO 720 idc=nkfdc,ikfdc,-1
            npdc(idc+1)=npdc(idc)
            DO 710 i=0,8
              kfdc(idc+1,i)=kfdc(idc,i)
  710       CONTINUE
  720     CONTINUE
          nkfdc=nkfdc+1
          kfdc(ikfdc,0)=nds
          DO 730 i=1,nds
            kfdc(ikfdc,i)=kfdm(i)
  730     CONTINUE
          npdc(ikfdc)=0
        ENDIF
        npdc(ikfdc)=npdc(ikfdc)+1
 
C...Write statistics on decay channels.
      ELSEIF(mtabu.EQ.52) THEN
        fac=1d0/max(1,nevdc)
        WRITE(mstu(11),5900) nevdc
        DO 750 idc=1,nkfdc
          DO 740 i=1,kfdc(idc,0)
            kfm=kfdc(idc,i)
            kf=(kfm+1)/2
            IF(2*kf.NE.kfm) kf=-kf
            CALL pyname(kf,chau)
            chdc(i)=chau(1:12)
            IF(chau(13:13).NE.' ') chdc(i)(12:12)='?'
  740     CONTINUE
          WRITE(mstu(11),6000) fac*npdc(idc),(chdc(i),i=1,kfdc(idc,0))
  750   CONTINUE
        IF(nredc.NE.0) WRITE(mstu(11),6100) fac*nredc
 
C...Copy statistics on decay channels into /PYJETS/.
      ELSEIF(mtabu.EQ.53) THEN
        fac=1d0/max(1,nevdc)
        DO 780 idc=1,nkfdc
          k(idc,1)=32
          k(idc,2)=99
          k(idc,3)=0
          k(idc,4)=0
          k(idc,5)=kfdc(idc,0)
          DO 760 j=1,5
            p(idc,j)=0d0
            v(idc,j)=0d0
  760     CONTINUE
          DO 770 i=1,kfdc(idc,0)
            kfm=kfdc(idc,i)
            kf=(kfm+1)/2
            IF(2*kf.NE.kfm) kf=-kf
            IF(i.LE.5) p(idc,i)=kf
            IF(i.GE.6) v(idc,i-5)=kf
  770     CONTINUE
          v(idc,5)=fac*npdc(idc)
  780   CONTINUE
        n=nkfdc
        DO 790 j=1,5
          k(n+1,j)=0
          p(n+1,j)=0d0
          v(n+1,j)=0d0
  790   CONTINUE
        k(n+1,1)=32
        k(n+1,2)=99
        k(n+1,5)=nevdc
        v(n+1,5)=fac*nredc
        mstu(3)=1
      ENDIF
 
C...Format statements for output on unit MSTU(11) (default 6).
 5000 FORMAT(///20x,'Event statistics - initial state'/
     &20x,'based on an analysis of ',i6,' events'//
     &3x,'Main flavours after',8x,'Fraction',4x,'Subfractions ',
     &'according to fragmenting system multiplicity'/
     &4x,'hard interaction',24x,'1',7x,'2',7x,'3',7x,'4',7x,'5',
     &6x,'6-7',5x,'8-10',3x,'11-15',3x,'16-25',4x,'>25'/)
 5100 FORMAT(3x,a12,1x,a12,f10.5,1x,10f8.4)
 5200 FORMAT(///20x,'Event statistics - final state'/
     &20x,'based on an analysis of ',i7,' events'//
     &5x,'Mean primary multiplicity =',f10.4/
     &5x,'Mean final   multiplicity =',f10.4/
     &5x,'Mean charged multiplicity =',f10.4//
     &5x,'Number of particles produced per event (directly and via ',
     &'decays/branchings)'/
     &8x,'KF    Particle/jet  MDCY',10x,'Particles',13x,'Antiparticles',
     &8x,'Total'/35x,'prim        seco        prim        seco'/)
 5300 FORMAT(1x,i9,4x,a16,i2,5(1x,f11.6))
 5400 FORMAT(///20x,'Factorial moments analysis of multiplicity'/
     &20x,'based on an analysis of ',i6,' events'//
     &3x,'delta-',a3,' delta-phi     <n>/bin',10x,'<F2>',18x,'<F3>',
     &18x,'<F4>',18x,'<F5>'/35x,4('     value     error  '))
 5500 FORMAT(10x)
 5600 FORMAT(2x,2f10.4,f12.4,4(f12.4,f10.4))
 5700 FORMAT(///20x,'Energy-Energy Correlation and Asymmetry'/
     &20x,'based on an analysis of ',i6,' events'//
     &2x,'theta range',8x,'EEC(theta)',8x,'EEC(180-theta)',7x,
     &'EECA(theta)'/2x,'in degrees ',3('      value    error')/)
 5800 FORMAT(2x,f4.1,' - ',f4.1,3(f11.4,f9.4))
 5900 FORMAT(///20x,'Decay channel analysis - final state'/
     &20x,'based on an analysis of ',i6,' events'//
     &2x,'Probability',10x,'Complete final state'/)
 6000 FORMAT(2x,f9.5,5x,8(a12,1x))
 6100 FORMAT(2x,f9.5,5x,'into other channels (more than 8 particles ',
     &'or table overflow)')
 
      RETURN
      END