lutabu.f

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

    
C*********************************************************************  
    
      SUBROUTINE lutabu(MTABU)  
    
C...Purpose: to evaluate various properties of an event, with   
C...statistics accumulated during the course of the run and 
C...printed at the end. 
      common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
      SAVE /lujets/ 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)    
      SAVE /ludat2/ 
      common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)    
      SAVE /ludat3/ 
      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*0./,fm2fm/120*0./,   
     &nevee/0/,fe1ec/50*0./,fe2ec/50*0./,fe1ea/25*0./,fe2ea/25*0./, 
     &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 120 i=nkfis,ikfis,-1   
          kfis(i+1,1)=kfis(i,1) 
          kfis(i+1,2)=kfis(i,2) 
          DO 120 j=0,10 
  120     npis(i+1,j)=npis(i,j) 
          nkfis=nkfis+1 
          kfis(ikfis,1)=kfmn    
          kfis(ikfis,2)=kfmx    
          DO 130 j=0,10 
  130     npis(ikfis,j)=0   
        ENDIF   
        npis(ikfis,0)=npis(ikfis,0)+1   
    
C...Count number of partons in initial state.   
        np=0    
        DO 150 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  
  140     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 140    
          ENDIF 
        ENDIF   
  150   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=1./max(1,nevis) 
        WRITE(mstu(11),1000) nevis  
        DO 160 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 luname(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 luname(kfm2,chau)  
        chis(2)=chau(1:12)  
        IF(chau(13:13).NE.' ') chis(2)(12:12)='?'   
  160   WRITE(mstu(11),1100) chis(1),chis(2),fac*npis(i,0), 
     &  (npis(i,j)/float(npis(i,0)),j=1,10) 
    
C...Copy statistics on initial parton state into /LUJETS/.  
      ELSEIF(mtabu.EQ.13) THEN  
        fac=1./max(1,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   
        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 170 j=1,5    
        p(i,j)=fac*npis(i,j)    
  170   v(i,j)=fac*npis(i,j+5)  
        n=nkfis 
        DO 180 j=1,5    
        k(n+1,j)=0  
        p(n+1,j)=0. 
  180   v(n+1,j)=0. 
        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 230 i=1,n    
        IF(k(i,1).LE.0.OR.k(i,1).GT.20.OR.k(i,1).EQ.13) goto 230    
        mstu(62)=mstu(62)+1 
        kc=lucomp(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(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=lucomp(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(luchge(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 190 ip=1,nkffs   
        IF(kfa.EQ.kffs(ip)) THEN    
          ikffs=-ip 
          goto 200  
        ELSEIF(kfa.LT.kffs(ip)) THEN    
          ikffs=ip  
          goto 200  
        ENDIF   
  190   CONTINUE    
        ikffs=nkffs+1   
  200   IF(ikffs.LT.0) THEN 
          ikffs=-ikffs  
        ELSE    
          IF(nkffs.GE.400) RETURN   
          DO 210 ip=nkffs,ikffs,-1  
          kffs(ip+1)=kffs(ip)   
          DO 210 j=1,4  
  210     npfs(ip+1,j)=npfs(ip,j)   
          nkffs=nkffs+1 
          kffs(ikffs)=kfa   
          DO 220 j=1,4  
  220     npfs(ikffs,j)=0   
        ENDIF   
        npfs(ikffs,kfs)=npfs(ikffs,kfs)+1   
  230   CONTINUE    
    
C...Write statistics on particle/parton composition of events.  
      ELSEIF(mtabu.EQ.22) THEN  
        fac=1./max(1,nevfs) 
        WRITE(mstu(11),1200) nevfs,fac*nprfs,fac*nfifs,fac*nchfs    
        DO 240 i=1,nkffs    
        CALL luname(kffs(i),chau)   
        kc=lucomp(kffs(i))  
        mdcyf=0 
        IF(kc.NE.0) mdcyf=mdcy(kc,1)    
  240   WRITE(mstu(11),1300) kffs(i),chau,mdcyf,(fac*npfs(i,j),j=1,4),  
     &  fac*(npfs(i,1)+npfs(i,2)+npfs(i,3)+npfs(i,4))   
    
C...Copy particle/parton composition information into /LUJETS/. 
      ELSEIF(mtabu.EQ.23) THEN  
        fac=1./max(1,nevfs) 
        DO 260 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 250 j=1,4    
        p(i,j)=fac*npfs(i,j)    
  250   v(i,j)=0.   
        p(i,5)=fac*k(i,5)   
  260   v(i,5)=0.   
        n=nkffs 
        DO 270 j=1,5    
        k(n+1,j)=0  
        p(n+1,j)=0. 
  270   v(n+1,j)=0. 
        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 280 im=1,3   
        DO 280 ib=1,10  
        DO 280 ip=1,4   
        fm1fm(im,ib,ip)=0.  
  280   fm2fm(im,ib,ip)=0.  
    
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 360 i=1,n    
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 360    
        IF(mstu(41).GE.2) THEN  
          kc=lucomp(k(i,2)) 
          IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.    
     &    kc.EQ.18) goto 360    
          IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0) 
     &    goto 360  
        ENDIF   
        pmr=0.  
        IF(mstu(42).EQ.1.AND.k(i,2).NE.22) pmr=ulmass(211)  
        IF(mstu(42).GE.2) pmr=p(i,5)    
        pr=max(1e-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),    
     &  1e20)),p(i,3))  
        IF(abs(yeta).GT.paru(57)) goto 360  
        phi=ulangl(p(i,1),p(i,2))   
        iyeta=512.*(yeta+paru(57))/(2.*paru(57))    
        iyeta=max(0,min(511,iyeta)) 
        iphi=512.*(phi+paru(1))/paru(2) 
        iphi=max(0,min(511,iphi))   
        iyep=0  
        DO 290 ib=0,9   
  290   iyep=iyep+4**ib*(2*mod(iyeta/2**ib,2)+mod(iphi/2**ib,2))    
    
C...Order particles in (pseudo)rapidity and/or azimuth. 
        IF(nupp.GT.mstu(4)-5-mstu(32)) THEN 
          CALL luerrm(11,'(LUTABU:) no more memory left in LUJETS') 
          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 300 i1=nupp-1,nlow+1,-1    
          IF(iyeta.GE.k(i1,1)) goto 310 
  300     k(i1+1,1)=k(i1,1) 
  310     k(i1+1,1)=iyeta   
          DO 320 i1=nupp-1,nlow+1,-1    
          IF(iphi.GE.k(i1,2)) goto 330  
  320     k(i1+1,2)=k(i1,2) 
  330     k(i1+1,2)=iphi    
          DO 340 i1=nupp-1,nlow+1,-1    
          IF(iyep.GE.k(i1,3)) goto 350  
  340     k(i1+1,3)=k(i1,3) 
  350     k(i1+1,3)=iyep    
        ENDIF   
  360   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 400 im=1,3   
        DO 370 ib=1,10  
        DO 370 ip=1,4   
  370   fevfm(ib,ip)=0. 
        DO 380 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 380 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)+2.  
          ELSEIF(nagr.EQ.3) THEN    
            fevfm(ib,1)=fevfm(ib,1)+6.  
            fevfm(ib,2)=fevfm(ib,2)+6.  
          ELSEIF(nagr.EQ.4) THEN    
            fevfm(ib,1)=fevfm(ib,1)+12. 
            fevfm(ib,2)=fevfm(ib,2)+24. 
            fevfm(ib,3)=fevfm(ib,3)+24. 
          ELSE  
            fevfm(ib,1)=fevfm(ib,1)+nagr*(nagr-1.)  
            fevfm(ib,2)=fevfm(ib,2)+nagr*(nagr-1.)*(nagr-2.)    
            fevfm(ib,3)=fevfm(ib,3)+nagr*(nagr-1.)*(nagr-2.)*(nagr-3.)  
            fevfm(ib,4)=fevfm(ib,4)+nagr*(nagr-1.)*(nagr-2.)*(nagr-3.)* 
     &      (nagr-4.)   
          ENDIF 
          iagr=icut 
          nagr=1    
        ENDIF   
  380   CONTINUE    
    
C...Add results to total statistics.    
        DO 390 ib=10,1,-1   
        DO 390 ip=1,4   
        IF(fevfm(1,ip).LT.0.5) THEN 
          fevfm(ib,ip)=0.   
        ELSEIF(im.LE.2) THEN    
          fevfm(ib,ip)=2**((ib-1)*ip)*fevfm(ib,ip)/fevfm(1,ip)  
        ELSE    
          fevfm(ib,ip)=4**((ib-1)*ip)*fevfm(ib,ip)/fevfm(1,ip)  
        ENDIF   
        fm1fm(im,ib,ip)=fm1fm(im,ib,ip)+fevfm(ib,ip)    
  390   fm2fm(im,ib,ip)=fm2fm(im,ib,ip)+fevfm(ib,ip)**2 
  400   CONTINUE    
        nmufm=nmufm+(nupp-nlow) 
        mstu(62)=nupp-nlow  
    
C...Write accumulated statistics on factorial moments.  
      ELSEIF(mtabu.EQ.32) THEN  
        fac=1./max(1,nevfm) 
        IF(mstu(42).LE.0) WRITE(mstu(11),1400) nevfm,'eta'  
        IF(mstu(42).EQ.1) WRITE(mstu(11),1400) nevfm,'ypi'  
        IF(mstu(42).GE.2) WRITE(mstu(11),1400) nevfm,'y  '  
        DO 420 im=1,3   
        WRITE(mstu(11),1500)    
        DO 420 ib=1,10  
        byeta=2.*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/float(2**(ib-1))    
        IF(im.EQ.3) bnave=fac*nmufm/float(4**(ib-1))    
        DO 410 ip=1,4   
        fmoma(ip)=fac*fm1fm(im,ib,ip)   
  410   fmoms(ip)=sqrt(max(0.,fac*(fac*fm2fm(im,ib,ip)-fmoma(ip)**2)))  
  420   WRITE(mstu(11),1600) byeta,bphi,bnave,(fmoma(ip),fmoms(ip), 
     &  ip=1,4) 
    
C...Copy statistics on factorial moments into /LUJETS/. 
      ELSEIF(mtabu.EQ.33) THEN  
        fac=1./max(1,nevfm) 
        DO 430 im=1,3   
        DO 430 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)=2.*paru(57)/k(i,3)   
        v(i,1)=paru(2)/k(i,4)   
        DO 430 ip=1,4   
        p(i,ip+1)=fac*fm1fm(im,ib,ip)   
  430   v(i,ip+1)=sqrt(max(0.,fac*(fac*fm2fm(im,ib,ip)-p(i,ip+1)**2)))  
        n=30    
        DO 440 j=1,5    
        k(n+1,j)=0  
        p(n+1,j)=0. 
  440   v(n+1,j)=0. 
        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 450 j=1,25   
        fe1ec(j)=0. 
        fe2ec(j)=0. 
        fe1ec(51-j)=0.  
        fe2ec(51-j)=0.  
        fe1ea(j)=0. 
  450   fe2ea(j)=0. 
    
C...Find particles to include, with proper assumed mass.    
      ELSEIF(mtabu.EQ.41) THEN  
        nevee=nevee+1   
        nlow=n+mstu(3)  
        nupp=nlow   
        ecm=0.  
        DO 460 i=1,n    
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 460    
        IF(mstu(41).GE.2) THEN  
          kc=lucomp(k(i,2)) 
          IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.    
     &    kc.EQ.18) goto 460    
          IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0) 
     &    goto 460  
        ENDIF   
        pmr=0.  
        IF(mstu(42).EQ.1.AND.k(i,2).NE.22) pmr=ulmass(211)  
        IF(mstu(42).GE.2) pmr=p(i,5)    
        IF(nupp.GT.mstu(4)-5-mstu(32)) THEN 
          CALL luerrm(11,'(LUTABU:) no more memory left in LUJETS') 
          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(1e-10,sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2))    
        ecm=ecm+p(nupp,4)   
  460   CONTINUE    
        IF(nupp.EQ.nlow) RETURN 
    
C...Analyze Energy-Energy Correlation in event. 
        fac=(2./ecm**2)*50./paru(1) 
        DO 470 j=1,50   
  470   fevee(j)=0. 
        DO 480 i1=nlow+2,nupp   
        DO 480 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(-1.,min(1.,cthe))) 
        ithe=max(1,min(50,1+int(50.*the/paru(1))))  
  480   fevee(ithe)=fevee(ithe)+fac*p(i1,4)*p(i2,4) 
        DO 490 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))    
  490   fe2ea(j)=fe2ea(j)+(fevee(51-j)-fevee(j))**2 
        mstu(62)=nupp-nlow  
    
C...Write statistics on Energy-Energy Correlation.  
      ELSEIF(mtabu.EQ.42) THEN  
        fac=1./max(1,nevee) 
        WRITE(mstu(11),1700) nevee  
        DO 500 j=1,25   
        feec1=fac*fe1ec(j)  
        fees1=sqrt(max(0.,fac*(fac*fe2ec(j)-feec1**2))) 
        feec2=fac*fe1ec(51-j)   
        fees2=sqrt(max(0.,fac*(fac*fe2ec(51-j)-feec2**2)))  
        feeca=fac*fe1ea(j)  
        feesa=sqrt(max(0.,fac*(fac*fe2ea(j)-feeca**2))) 
  500   WRITE(mstu(11),1800) 3.6*(j-1),3.6*j,feec1,fees1,feec2,fees2,   
     &  feeca,feesa 
    
C...Copy statistics on Energy-Energy Correlation into /LUJETS/. 
      ELSEIF(mtabu.EQ.43) THEN  
        fac=1./max(1,nevee) 
        DO 510 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(0.,fac*(fac*fe2ec(i)-p(i,1)**2)))   
        p(i,2)=fac*fe1ec(51-i)  
        v(i,2)=sqrt(max(0.,fac*(fac*fe2ec(51-i)-p(i,2)**2)))    
        p(i,3)=fac*fe1ea(i) 
        v(i,3)=sqrt(max(0.,fac*(fac*fe2ea(i)-p(i,3)**2)))   
        p(i,4)=paru(1)*(i-1)/50.    
        p(i,5)=paru(1)*i/50.    
        v(i,4)=3.6*(i-1)    
  510   v(i,5)=3.6*i    
        n=25    
        DO 520 j=1,5    
        k(n+1,j)=0  
        p(n+1,j)=0. 
  520   v(n+1,j)=0. 
        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 550 i=1,n    
        IF(k(i,1).LE.0.OR.k(i,1).GE.6) goto 550 
        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 530 ids=nds-1,1,-1   
        iin=ids+1   
        IF(kfm.LT.kfdm(ids)) goto 540   
  530   kfdm(ids+1)=kfdm(ids)   
        iin=1   
  540   kfdm(iin)=kfm   
  550   CONTINUE    
    
C...Find whether old or new final state.    
        DO 570 idc=1,nkfdc  
        IF(nds.LT.kfdc(idc,0)) THEN 
          ikfdc=idc 
          goto 580  
        ELSEIF(nds.EQ.kfdc(idc,0)) THEN 
          DO 560 i=1,nds    
          IF(kfdm(i).LT.kfdc(idc,i)) THEN   
            ikfdc=idc   
            goto 580    
          ELSEIF(kfdm(i).GT.kfdc(idc,i)) THEN   
            goto 570    
          ENDIF 
  560     CONTINUE  
          ikfdc=-idc    
          goto 580  
        ENDIF   
  570   CONTINUE    
        ikfdc=nkfdc+1   
  580   IF(ikfdc.LT.0) THEN 
          ikfdc=-ikfdc  
        ELSEIF(nkfdc.GE.200) THEN   
          nredc=nredc+1 
          RETURN    
        ELSE    
          DO 590 idc=nkfdc,ikfdc,-1 
          npdc(idc+1)=npdc(idc) 
          DO 590 i=0,8  
  590     kfdc(idc+1,i)=kfdc(idc,i) 
          nkfdc=nkfdc+1 
          kfdc(ikfdc,0)=nds 
          DO 600 i=1,nds    
  600     kfdc(ikfdc,i)=kfdm(i) 
          npdc(ikfdc)=0 
        ENDIF   
        npdc(ikfdc)=npdc(ikfdc)+1   
    
C...Write statistics on decay channels. 
      ELSEIF(mtabu.EQ.52) THEN  
        fac=1./max(1,nevdc) 
        WRITE(mstu(11),1900) nevdc  
        DO 620 idc=1,nkfdc  
        DO 610 i=1,kfdc(idc,0)  
        kfm=kfdc(idc,i) 
        kf=(kfm+1)/2    
        IF(2*kf.NE.kfm) kf=-kf  
        CALL luname(kf,chau)    
        chdc(i)=chau(1:12)  
  610   IF(chau(13:13).NE.' ') chdc(i)(12:12)='?'   
  620   WRITE(mstu(11),2000) fac*npdc(idc),(chdc(i),i=1,kfdc(idc,0))    
        IF(nredc.NE.0) WRITE(mstu(11),2100) fac*nredc   
    
C...Copy statistics on decay channels into /LUJETS/.    
      ELSEIF(mtabu.EQ.53) THEN  
        fac=1./max(1,nevdc) 
        DO 650 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 630 j=1,5    
        p(idc,j)=0. 
  630   v(idc,j)=0. 
        DO 640 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  
  640   IF(i.GE.6) v(idc,i-5)=kf    
  650   v(idc,5)=fac*npdc(idc)  
        n=nkfdc 
        DO 660 j=1,5    
        k(n+1,j)=0  
        p(n+1,j)=0. 
  660   v(n+1,j)=0. 
        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).  
 1000 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'/)   
 1100 FORMAT(3x,a12,1x,a12,f10.5,1x,10f8.4) 
 1200 FORMAT(///20x,'Event statistics - final state'/   
     &20x,'based on an analysis of ',i6,' events'// 
     &5x,'Mean primary multiplicity =',f8.3/    
     &5x,'Mean final   multiplicity =',f8.3/    
     &5x,'Mean charged multiplicity =',f8.3//   
     &5x,'Number of particles produced per event (directly and via ',   
     &'decays/branchings)'/ 
     &5x,'KF    Particle/jet  MDCY',8x,'Particles',9x,'Antiparticles',  
     &5x,'Total'/34x,'prim      seco      prim      seco'/) 
 1300 FORMAT(1x,i6,4x,a16,i2,5(1x,f9.4))    
 1400 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  '))    
 1500 FORMAT(10x)   
 1600 FORMAT(2x,2f10.4,f12.4,4(f12.4,f10.4))    
 1700 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')/)    
 1800 FORMAT(2x,f4.1,' - ',f4.1,3(f11.4,f9.4))  
 1900 FORMAT(///20x,'Decay channel analysis - final state'/ 
     &20x,'based on an analysis of ',i6,' events'// 
     &2x,'Probability',10x,'Complete final state'/) 
 2000 FORMAT(2x,f9.5,5x,8(a12,1x))  
 2100 FORMAT(2x,f9.5,5x,'into other channels (more than 8 particles ',  
     &'or table overflow)') 
    
      RETURN    
      END