      program pp
c* makes web page results for W mass (Aleph) . Adapted from G. Quast program
c* reads input file : wm_inputs.cards
c* to compile and linkit: 
c* f77 filename.f -L/cern/pro/lib -lmathlib -lpacklib -o filename.exe
c* output web page : mass-test.html

      call readcards

      call aver_lvqq            !lvqq per ecm
      call xpect_lvqq           !expected lvqq error per ecm
      call aver_0               !lvqq+4q per ecm (e+mu+tau+4q)
c      call aver_0a              !lvqq+4q per ecm (lv+4q)
      call aver_0e              !expected lvqq+4q error per ecm
      call aver_1               !each channel for all ecm
      call aver_1e              !expected error per channel for all ecm
      call aver_2               !all lvqq
      call aver_2e              !expected all lvqq
      call aver_all             !averages everything
c      call aver_alla             !averages everything(inputs from aver_1)
      call aver_1_19x           !each channel for 192 to 202 GeV data
      call aver_1_20x           !each channel for 202 to 209 GeV data
      call aver_all19x          !average 192 to 202 GeV data
      call aver_all20x          !average 202 to 207 GeV data
      call aver_219x            !all lvqq at 19x energies
      call aver_220x            !all lvqq at 20x energies
      call aver_1e19x           !expected error per channel for all ecm
      call aver_1e20x           !expected error per channel for all ecm
      call aver_2e19x           !expected 19x lvqq
      call aver_2e20x           !expected 20x lvqq

      call writweb
      end
c
c***********************************************************

      subroutine readcards
      implicit none
      include 'pp.inc'

      integer i,j,iecm,itypsys,iene
      integer nenesys,ienesys(necmx),jene,jene1,jene2,jecm
      real esystmp(nchan),rhotmp,emin
      character*80 tit
      character*10 titr
      lout    = 6                ! output
      lin     = 8                ! input

      call vzero(irange_sys,necmx*necmx)

      OPEN(Lin,FILE='wm_inputs.cards',STATUS='unknown')
      REWIND Lin

      iecm=0

 1    READ (Lin,100) TIT        
      WRITE(LOut,100) TIT
      if(tit(1:3).eq.'Ecm') then !Ecm lumi
        iecm=iecm+1
        necm=iecm
        READ (Lin,*) ecm(iecm), lumi(iecm)
        WRITE(LOut,*) ecm(iecm), lumi(iecm)
        read(lin,100) tit       ! n1 n2 dmexp m dmstat
        write(lout,100) tit
        do i=1,4                !loop e mu tau 4q
        read(lin,*) n1(i,iecm),n2(i,iecm),dmxpect(i,iecm)
     &            , xmw(i,iecm),estat(i,iecm)
        write(lout,*) n1(i,iecm),n2(i,iecm),dmxpect(i,iecm)
     &            , xmw(i,iecm),estat(i,iecm)
        enddo
      elseif(tit(1:4).eq.'Syst') then !type of sys error

c* check the list of Ecm for which the following syst errors are valid
        nenesys=0          
        do jene=1,necm
          jene1=6+(jene-1)*4
          jene2=jene1+2
          read(tit(jene1:jene2),*,err=998) iene
c          print *,'iene ',iene
          if(iene.gt.0) then
            emin=9999.
            do iecm=1,necm
              if(abs(ecm(iecm)-float(iene)).lt.emin)then
                emin=abs(ecm(iecm)-float(iene))
                jecm=iecm
              endif
            enddo
            nenesys=nenesys+1
            ienesys(nenesys)=jecm
c            print *,'nenesys ienesys ',nenesys,jecm
          endif
        enddo
 998    continue
c        print *,'nenesys ienesys ',nenesys,ienesys
        nrange_sys(ienesys(1))=nenesys
        do jene=1,nenesys
          irange_sys(jene,ienesys(1))=ienesys(jene)
        enddo
 2      read(lin,'(A10,5f5.1)') titr
     &        , (esystmp(i),i=1,4), rhotmp
         if(titr(1:6).eq.'ENDSYS') goto 1
         do i=1,ntypsys
            if(titr(1:10).eq.typsys(i))itypsys=i
         enddo
         rho(itypsys)=rhotmp
         do j=1,nenesys
           do i=1,4
             esys(itypsys,i,ienesys(j))=esystmp(i)/1000. !in GeV
           enddo
           write(lout,'(A,5f5.1)') typsys(itypsys)
     &          ,(esys(itypsys,i,ienesys(j))*1000.,i=1,4), rho(itypsys)
         enddo
         goto 2
      elseif(tit(1:7).eq.'CorrEcm') then
         do j=1,necmx
            read(lin,*) ecmc(j),decmc(j),(rhoecm(i,j),i=1,necmx)
            write(lout,'(2f6.1,8f5.2)')ecmc(j),decmc(j)
     &           ,(rhoecm(i,j),i=1,necmx)
            decmc(j)=decmc(j)/1000. !in GeV
         enddo
c         print *,'decmc',decmc
      elseif(tit(1:3).eq.'END') then

        write(lout,*)' ------------------------------------------'
        write(lout,*)'|       Data Cards Reading Finished        |'
        write(lout,*)' ------------------------------------------'
        write(lout,*)' '
        return  
      endif
      goto 1

 100  format(A80)
c 999  write(lout,*)'bad format in input file'
c      stop
      end
 
      subroutine aver_alla
c* average lvqq and 4q results for all c.m. energies (and m4q-mlvqq) 
c*from already combined lvqq and 4q results     
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=nchan,nerr=3)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene
      real emin

      integer ndfit,ndfitmx
      parameter (ndfitmx=36)
      real wmtofit(ndfitmx)
      integer ilq4q(ndfitmx)
      double precision dcovif(ndfitmx,ndfitmx)
      common/forfits/ndfit,ilq4q,wmtofit
      common/forfitd/dcovif

      print *,'|-------------------------------------------------------' 
      print *,'|                    AVERAGE ALL (2nd way)             |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,nd)
      nda=0
      xm_all=0.
      em_all=0.
      es_all1=0.
      es_all2=0.
      do j=1,4                  !loop over lvqq+4q results at this ecm
        if(xm_chan(j).gt.0.) then
          nda=nda+1
          iener(j)=1.
          wm(nda)=xm_chan(j)
          if(nda.le.ndfitmx) then
            wmtofit(nda)=xm_chan(j)
            ilq4q(nda)=1
            if(j.eq.4) ilq4q(nda)=2
          endif
          err(1,nda)=exm_chan(j)
          rhoerr(1)=0.
          err(2,nda)=sxm_chan1(j)
          rhoerr(2)=1.
          err(3,nda)=sxm_chan2(j)
          rhoerr(3)=0.
        endif                   !check if result is available
      enddo                     !loop on channels
      print *,' wm wfit'
      do i=1,nda
        print *,wm(i),wmtofit(i)
      enddo
      call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &     ,cov,covec)
      call covinv(nda,nd,cov,dcovi,w)
      print *,'nda nd weights ',nda,nd,w
      call aver(nda,nd,wm,cov,w,nerr,err,covec
     &     ,wmav,wmer,e2,0)
      print *,'wmav er ',wmav,wmer
      xm_all=wmav
      em_all=sqrt(e2(1))
      es_all1=sqrt(wmer**2-e2(1)-e2(3))
      es_all2=sqrt(e2(3))
      print *,'ALL m em es'
      print *,xm_all,em_all,es_all1,es_all2

        if(nda.ge.4) then
          ndfit=nda
          do j=1,nd
            do k=1,nd
              dcovif(j,k)=dcovi(j,k)
            enddo
          enddo
          call wfits(mwdiftot,emwdiftot)
        endif

      end
      subroutine aver_all
c* average lvqq and 4q results for all c.m. energies (and m4q-mlvqq)     
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=nchan*necmx,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene,ndaw
      real emin

      integer ndfit,ndfitmx
      parameter (ndfitmx=36)
      real wmtofit(ndfitmx)
      integer ilq4q(ndfitmx)
      double precision dcovif(ndfitmx,ndfitmx)
      common/forfits/ndfit,ilq4q,wmtofit
      common/forfitd/dcovif

      print *,'|-------------------------------------------------------' 
      print *,'|                    AVERAGE ALL                       |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,nd)
      nda=0
      xm_all=0.
      em_all=0.
      es_all1=0.
      es_all2=0.
      do i=1,necm
        emin=9999.
        do ii=1,necmx           !check index of c.m. energy
          if(abs(ecm(i)-ecmc(ii)).lt.emin) then
            iene=ii
            emin=abs(ecm(i)-ecmc(ii))
          endif
        enddo
c        print *,'i iene ',i,iene
        do j=1,nchan                !loop over lvqq+4q results at this ecm
c          print *,'j xmw ',j,xmw(j,i)
          if(xmw(j,i).gt.0.) then
            nda=nda+1
            iener(nda)=iene
            wm(nda)=xmw(j,i)
            if(nda.le.ndfitmx) then
              wmtofit(nda)=xmw(j,i)
              ilq4q(nda)=1
              if(j.eq.4) ilq4q(nda)=2
            endif
            err(1,nda)=estat(j,i)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
c              if(k.eq.10)print *,'averall iene decmc esys '
c     &             ,iene,decmc(iene),esys(k,j)
              err(k+1,nda)=esys(k,j,iene)
              rhoerr(k+1)=rho(k)
            enddo
          endif                 !check if result is available
        enddo                   !loop on channels
      enddo                     !loop on c.m. energies
      print *,'iene ecm wm wfit'
      do i=1,nda
        print *,iener(i),ecm(iener(i)),wm(i),wmtofit(i)
      enddo
      call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &     ,cov,covec)
      call covinv(nda,nd,cov,dcovi,w)
      call aver(nda,nd,wm,cov,w,nerr,err,covec
     &     ,wmav,wmer,e2,1)
      xm_all=wmav
      em_all=sqrt(e2(1))
      es_all1=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
      es_all2=sqrt(e2(14)+e2(15))
      chi2_all=0
      do i=1,nda
c        pull_all(i) = (wm(i)-wmav) * dcovi(i,i) * (wm(i)-wmav)
        pull_all(i) = (wm(i)-wmav) * sqrt(dcovi(i,i))
        do j=1,nda
          chi2_all=chi2_all + (wm(i)-wmav) * dcovi(i,j) * (wm(j)-wmav)
        enddo
      enddo
      write(6,'(A,f5.2,A,i2)') 'ALL chi2/ndf :',chi2_all,'/',nda-1
      print *,'pulls'
      write(70,'(1f7.2)') (pull_all(i),i=1,nda)
      print *,'ALL m em es'
      print *,xm_all,em_all,es_all1,es_all2
      call vzero(wei_mall,nchan)
      ndaw=0
      do i=1,necm
        do j=1,nchan
          if(xmw(j,i).gt.0.) then
            ndaw=ndaw+1
            wei_mall(j)=w(ndaw)+wei_mall(j)
          endif
        enddo
      enddo
      
        if(nda.ge.4) then
          ndfit=nda
c* switch off FSR error for mass difference fit
          do j=1,nd
             err(14,j)=0.
             err(15,j)=0.
          enddo
          call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &         ,cov,covec)
          call covinv(nda,nd,cov,dcovi,w)
          do j=1,nd
            do k=1,nd
              dcovif(j,k)=dcovi(j,k)
            enddo
          enddo
          call wfits(mwdiftot,emwdiftot)
        endif

      end
      subroutine aver_all19x
c* average lvqq and 4q results for 192-196 GeV energies (and m4q-mlvqq)     
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=nchan*necmx,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene,ndaw
      real emin

      integer ndfit,ndfitmx
      parameter (ndfitmx=36)
      real wmtofit(ndfitmx)
      integer ilq4q(ndfitmx)
      double precision dcovif(ndfitmx,ndfitmx)
      common/forfits/ndfit,ilq4q,wmtofit
      common/forfitd/dcovif

      print *,'|-------------------------------------------------------' 
      print *,'|                    AVERAGE 19X                       |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,nd)
      nda=0
      xm_all19x=0.
      em_all19x=0.
      es_all119x=0.
      es_all219x=0.
      do i=1,necm
        if(ecm(i).gt.190..and.ecm(i).lt.202.) then
        emin=9999.
        do ii=1,necmx           !check index of c.m. energy
          if(abs(ecm(i)-ecmc(ii)).lt.emin) then
            iene=ii
            emin=abs(ecm(i)-ecmc(ii))
          endif
        enddo
        do j=1,nchan                !loop over lvqq+4q results at this ecm
          if(xmw(j,i).gt.0.) then
            nda=nda+1
            iener(nda)=iene
            wm(nda)=xmw(j,i)
            if(nda.le.ndfitmx) then
              wmtofit(nda)=xmw(j,i)
              ilq4q(nda)=1
              if(j.eq.4) ilq4q(nda)=2
            endif
            err(1,nda)=estat(j,i)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,j,iene)
              rhoerr(k+1)=rho(k)
            enddo
          endif                 !check if result is available
        enddo                   !loop on channels
      endif                     !check 19x cm energy range
      enddo                     !loop on c.m. energies
      print *,'iene ecm wm wfit'
      do i=1,nda
        print *,iener(i),ecm(iener(i)),wm(i),wmtofit(i)
      enddo
      call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &     ,cov,covec)
      call covinv(nda,nd,cov,dcovi,w)
      call aver(nda,nd,wm,cov,w,nerr,err,covec
     &     ,wmav,wmer,e2,0)
      xm_all19x=wmav
      em_all19x=sqrt(e2(1))
      es_all119x=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
      es_all219x=sqrt(e2(14)+e2(15))
      print *,'ALL m em es'
      print *,xm_all19x,em_all19x,es_all119x,es_all219x
      call vzero(wei_mall19x,nchan)
      ndaw=0
      do i=1,necm
        if(ecm(i).gt.190..and.ecm(i).lt.202.) then
        do j=1,nchan
          if(xmw(j,i).gt.0.) then
            ndaw=ndaw+1
            wei_mall19x(j)=w(ndaw)+wei_mall19x(j)
          endif
        enddo
        endif
      enddo
      
        if(nda.ge.4) then
          ndfit=nda
c* switch off FSR eror for mass difference fit
          do j=1,nd
             err(14,j)=0.
             err(15,j)=0.
          enddo
          call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &         ,cov,covec)
          call covinv(nda,nd,cov,dcovi,w)
          do j=1,nd
            do k=1,nd
              dcovif(j,k)=dcovi(j,k)
            enddo
          enddo
          call wfits(mwdif19x,emwdif19x)
        endif

      end
      subroutine aver_all20x
c* average lvqq and 4q results for 192-196 GeV energies (and m4q-mlvqq)     
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=nchan*necmx,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene,ndaw
      real emin

      integer ndfit,ndfitmx
      parameter (ndfitmx=36)
      real wmtofit(ndfitmx)
      integer ilq4q(ndfitmx)
      double precision dcovif(ndfitmx,ndfitmx)
      common/forfits/ndfit,ilq4q,wmtofit
      common/forfitd/dcovif

      print *,'|-------------------------------------------------------' 
      print *,'|                    AVERAGE 20X                       |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,nd)
      nda=0
      xm_all20x=0.
      em_all20x=0.
      es_all120x=0.
      es_all220x=0.
      do i=1,necm
        if(ecm(i).gt.202..and.ecm(i).lt.209.) then
        emin=9999.
        do ii=1,necmx           !check index of c.m. energy
          if(abs(ecm(i)-ecmc(ii)).lt.emin) then
            iene=ii
            emin=abs(ecm(i)-ecmc(ii))
          endif
        enddo
        do j=1,nchan                !loop over lvqq+4q results at this ecm
          if(xmw(j,i).gt.0.) then
            nda=nda+1
            iener(nda)=iene
            wm(nda)=xmw(j,i)
            if(nda.le.ndfitmx) then
              wmtofit(nda)=xmw(j,i)
              ilq4q(nda)=1
              if(j.eq.4) ilq4q(nda)=2
            endif
            err(1,nda)=estat(j,i)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,j,iene)
              rhoerr(k+1)=rho(k)
            enddo
          endif                 !check if result is available
        enddo                   !loop on channels
      endif                     !check 20x cm energy range
      enddo                     !loop on c.m. energies
      print *,'iene ecm wm wfit'
      do i=1,nda
        print *,iener(i),ecm(iener(i)),wm(i),wmtofit(i)
      enddo
      call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &     ,cov,covec)
      call covinv(nda,nd,cov,dcovi,w)
      call aver(nda,nd,wm,cov,w,nerr,err,covec
     &     ,wmav,wmer,e2,0)
      xm_all20x=wmav
      em_all20x=sqrt(e2(1))
      es_all120x=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
      es_all220x=sqrt(e2(14)+e2(15))
      print *,'ALL m em es'
      print *,xm_all20x,em_all20x,es_all120x,es_all220x
      call vzero(wei_mall20x,nchan)
      ndaw=0
      do i=1,necm
        if(ecm(i).gt.202..and.ecm(i).lt.209.) then
        do j=1,nchan
          if(xmw(j,i).gt.0.) then
            ndaw=ndaw+1
            wei_mall20x(j)=w(ndaw)+wei_mall20x(j)
          endif
        enddo
        endif
      enddo
      
        if(nda.ge.4) then
          ndfit=nda
c* switch off FSR eror for mass difference fit
          do j=1,nd
             err(14,j)=0.
             err(15,j)=0.
          enddo
          call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &         ,cov,covec)
          call covinv(nda,nd,cov,dcovi,w)
          do j=1,nd
            do k=1,nd
              dcovif(j,k)=dcovi(j,k)
            enddo
          enddo
          call wfits(mwdif20x,emwdif20x)
        endif

      end

      subroutine aver_2
c* average lvqq for all c.m. energies     
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=3*necmx,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|                    AVERAGE ALL lvqq                  |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,nd)
      nda=0
      xm_lvqq=0.
      em_lvqq=0.
      es_lvqq1=0.
      es_lvqq2=0.
      do i=1,necm
        emin=9999.
        do ii=1,necmx           !check index of c.m. energy
          if(abs(ecm(i)-ecmc(ii)).lt.emin) then
            iene=ii
            emin=abs(ecm(i)-ecmc(ii))
          endif
        enddo
        do j=1,3                !loop over lvqq+4q results at this ecm
          if(xmw(j,i).gt.0.) then
            nda=nda+1
            iener(nda)=iene
            wm(nda)=xmw(j,i)
            err(1,nda)=estat(j,i)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,j,iene)
              rhoerr(k+1)=rho(k)
            enddo
          endif                 !check if result is available
        enddo
      enddo                     !loop on c.m. energies
      call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &     ,cov,covec)
      call covinv(nda,nd,cov,dcovi,w)
      call aver(nda,nd,wm,cov,w,nerr,err,covec
     &     ,wmav,wmer,e2,0)
      xm_lvqq=wmav
      em_lvqq=sqrt(e2(1))
      es_lvqq1=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
      es_lvqq2=sqrt(e2(14)+e2(15))
      print *,'ALLlvqq  m em es'
      print *,xm_lvqq,em_lvqq,es_lvqq1,es_lvqq2
      end

      subroutine aver_2e
c* expected average lvqq for all c.m. energies     
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=3*necmx,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene,ndaw
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|          expected  AVERAGE ALL lvqq                  |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,nd)
      ndaw=0
      do i=1,necm
        do j=1,3
          if(xmw(j,i).gt.0.) ndaw=ndaw+1
        enddo
      enddo
      nda=0
      xm_lvqq_e=0.
      em_lvqq_e=0.
      es_lvqq1_e=0.
      es_lvqq2_e=0.
      do i=1,necm
        emin=9999.
        do ii=1,necmx           !check index of c.m. energy
          if(abs(ecm(i)-ecmc(ii)).lt.emin) then
            iene=ii
            emin=abs(ecm(i)-ecmc(ii))
          endif
        enddo
        do j=1,3                !loop over lvqq+4q results at this ecm
          if(dmxpect(j,i).gt.0.) then
            nda=nda+1
            iener(nda)=iene
            wm(nda)=80.35
            err(1,nda)=dmxpect(j,i)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,j,iene)
              rhoerr(k+1)=rho(k)
            enddo
          endif                 !check if result is available
        enddo
      enddo                     !loop on c.m. energies
      if(nda.gt.0.and.ndaw.eq.nda) then
        call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &       ,cov,covec)
        call covinv(nda,nd,cov,dcovi,w)
        call aver(nda,nd,wm,cov,w,nerr,err,covec
     &       ,wmav,wmer,e2,0)
        xm_lvqq_e=wmav
        em_lvqq_e=sqrt(e2(1))
        es_lvqq1_e=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
        es_lvqq2_e=sqrt(e2(14)+e2(15))
        print *,'ALLlvqq  m em es'
        print *,xm_lvqq_e,em_lvqq_e,es_lvqq1_e,es_lvqq2_e
      endif
      end

      subroutine aver_2e19x
c* expected average lvqq for all c.m. energies     
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=3*necmx,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene,ndaw
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|          expected  AVERAGE ALL lvqq 19X              |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,nd)
      ndaw=0
      do i=1,necm
        if(ecm(i).gt.190..and.ecm(i).lt.202.) then
          do j=1,3
            if(xmw(j,i).gt.0.) ndaw=ndaw+1
          enddo
        endif
      enddo
      nda=0
      xm_lvqq_e19x=0.
      em_lvqq_e19x=0.
      es_lvqq1_e19x=0.
      es_lvqq2_e19x=0.
      do i=1,necm
        if(ecm(i).gt.190..and.ecm(i).lt.202.) then
        emin=9999.
        do ii=1,necmx           !check index of c.m. energy
          if(abs(ecm(i)-ecmc(ii)).lt.emin) then
            iene=ii
            emin=abs(ecm(i)-ecmc(ii))
          endif
        enddo
        do j=1,3                !loop over lvqq+4q results at this ecm
          if(dmxpect(j,i).gt.0.) then
            nda=nda+1
            iener(nda)=iene
            wm(nda)=80.35
            err(1,nda)=dmxpect(j,i)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,j,iene)
              rhoerr(k+1)=rho(k)
            enddo
          endif                 !check if result is available
        enddo
      endif                     !check energy range 19x
      enddo                     !loop on c.m. energies
      if(nda.gt.0.and.ndaw.eq.nda) then
        call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &       ,cov,covec)
        call covinv(nda,nd,cov,dcovi,w)
        call aver(nda,nd,wm,cov,w,nerr,err,covec
     &       ,wmav,wmer,e2,0)
        xm_lvqq_e19x=wmav
        em_lvqq_e19x=sqrt(e2(1))
        es_lvqq1_e19x=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
        es_lvqq2_e19x=sqrt(e2(14)+e2(15))
        print *,'ALLlvqq  m em es'
        print *,xm_lvqq_e19x,em_lvqq_e19x,es_lvqq1_e19x,es_lvqq2_e19x
      endif
      end

      subroutine aver_2e20x
c* expected average lvqq for all c.m. energies     
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=3*necmx,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene,ndaw
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|          expected  AVERAGE ALL lvqq 20X              |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,nd)
      ndaw=0
      do i=1,necm
        if(ecm(i).gt.202..and.ecm(i).lt.209.) then
          do j=1,3
            if(xmw(j,i).gt.0.) ndaw=ndaw+1
          enddo
        endif
      enddo
      nda=0
      xm_lvqq_e20x=0.
      em_lvqq_e20x=0.
      es_lvqq1_e20x=0.
      es_lvqq2_e20x=0.
      do i=1,necm
        if(ecm(i).gt.202..and.ecm(i).lt.209.) then
        emin=9999.
        do ii=1,necmx           !check index of c.m. energy
          if(abs(ecm(i)-ecmc(ii)).lt.emin) then
            iene=ii
            emin=abs(ecm(i)-ecmc(ii))
          endif
        enddo
        do j=1,3                !loop over lvqq+4q results at this ecm
          if(dmxpect(j,i).gt.0.) then
            nda=nda+1
            iener(nda)=iene
            wm(nda)=80.35
            err(1,nda)=dmxpect(j,i)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,j,iene)
              rhoerr(k+1)=rho(k)
            enddo
          endif                 !check if result is available
        enddo
      endif                     !check energy range 20x
      enddo                     !loop on c.m. energies
      if(nda.gt.0.and.ndaw.eq.nda) then
        call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &       ,cov,covec)
        call covinv(nda,nd,cov,dcovi,w)
        call aver(nda,nd,wm,cov,w,nerr,err,covec
     &       ,wmav,wmer,e2,0)
        xm_lvqq_e20x=wmav
        em_lvqq_e20x=sqrt(e2(1))
        es_lvqq1_e20x=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
        es_lvqq2_e20x=sqrt(e2(14)+e2(15))
        print *,'ALLlvqq  m em es'
        print *,xm_lvqq_e20x,em_lvqq_e20x,es_lvqq1_e20x,es_lvqq2_e20x
      endif
      end

      subroutine aver_220X
c* average lvqq for all c.m. energies     
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=3*necmx,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|                    AVERAGE 20X lvqq                  |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,nd)
      nda=0
      xm_lvqq20x=0.
      em_lvqq20x=0.
      es_lvqq120x=0.
      es_lvqq220x=0.
      do i=1,necm
        if(ecm(i).gt.202..and.ecm(i).lt.209.) then
        emin=9999.
        do ii=1,necmx           !check index of c.m. energy
          if(abs(ecm(i)-ecmc(ii)).lt.emin) then
            iene=ii
            emin=abs(ecm(i)-ecmc(ii))
          endif
        enddo
        do j=1,3                !loop over lvqq+4q results at this ecm
          if(xmw(j,i).gt.0.) then
            nda=nda+1
            iener(nda)=iene
            wm(nda)=xmw(j,i)
            err(1,nda)=estat(j,i)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,j,iene)
              rhoerr(k+1)=rho(k)
            enddo
          endif                 !check if result is available
        enddo
      endif                     !check 19x ecm range
      enddo                     !loop on c.m. energies
      call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &     ,cov,covec)
      call covinv(nda,nd,cov,dcovi,w)
      call aver(nda,nd,wm,cov,w,nerr,err,covec
     &     ,wmav,wmer,e2,0)
      xm_lvqq20x=wmav
      em_lvqq20x=sqrt(e2(1))
      es_lvqq120x=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
      es_lvqq220x=sqrt(e2(14)+e2(15))
      print *,'ALLlvqq  m em es'
      print *,xm_lvqq20x,em_lvqq20x,es_lvqq120x,es_lvqq220x
      end


      subroutine aver_219X
c* average lvqq for all c.m. energies     
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=3*necmx,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|                    AVERAGE 19X lvqq                  |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,nd)
      nda=0
      xm_lvqq19x=0.
      em_lvqq19x=0.
      es_lvqq119x=0.
      es_lvqq219x=0.
      do i=1,necm
        if(ecm(i).gt.190..and.ecm(i).lt.202.) then
        emin=9999.
        do ii=1,necmx           !check index of c.m. energy
          if(abs(ecm(i)-ecmc(ii)).lt.emin) then
            iene=ii
            emin=abs(ecm(i)-ecmc(ii))
          endif
        enddo
        do j=1,3                !loop over lvqq+4q results at this ecm
          if(xmw(j,i).gt.0.) then
            nda=nda+1
            iener(nda)=iene
            wm(nda)=xmw(j,i)
            err(1,nda)=estat(j,i)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,j,iene)
              rhoerr(k+1)=rho(k)
            enddo
          endif                 !check if result is available
        enddo
      endif                     !check 19x ecm range
      enddo                     !loop on c.m. energies
      call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &     ,cov,covec)
      call covinv(nda,nd,cov,dcovi,w)
      call aver(nda,nd,wm,cov,w,nerr,err,covec
     &     ,wmav,wmer,e2,0)
      xm_lvqq19x=wmav
      em_lvqq19x=sqrt(e2(1))
      es_lvqq119x=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
      es_lvqq219x=sqrt(e2(14)+e2(15))
      print *,'ALLlvqq  m em es'
      print *,xm_lvqq19x,em_lvqq19x,es_lvqq119x,es_lvqq219x
      end

      subroutine aver_1_19x
c* average each channel results for 192 to 196 GeV c.m. energies     
      implicit none
      include 'pp.inc'
      integer nda,nerr
      parameter (nerr=ntypsys+1)
      integer iener(necmx)
      real wm(necmx),err(nerr,necmx),rhoerr(nerr)
      real cov(necmx,necmx), w(necmx)
     &     ,covec(necmx,necmx,nerr)
      double precision dcovi(necmx,necmx)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|               AVERAGE each channel 19X               |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,necmx)
      do i=1,nchan
c        print *,'19x:ichan=',i
        nda=0
        xm_chan19x(i)=0.
        exm_chan19x(i)=0.
        sxm_chan119x(i)=0.
        sxm_chan219x(i)=0.
        do j=1,necm         !loop over different ecm results for this channel
          if(ecm(j).gt.190..and.ecm(j).lt.202.) then
          emin=9999.
          do ii=1,necmx         !check index of c.m. energy
            if(abs(ecm(j)-ecmc(ii)).lt.emin) then
              iene=ii
              emin=abs(ecm(j)-ecmc(ii))
            endif
          enddo
          if(xmw(i,j).gt.0.) then
            nda=nda+1
            wm(nda)=xmw(i,j)
            err(1,nda)=estat(i,j)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,i,iene)=decmc(iene)*2./ecmc(iene)*80.35
               err(k+1,nda)=esys(k,i,iene)
               rhoerr(k+1)=rho(k)
c               print *,'19x: k err rho ',k,err(k+1,nda),rhoerr(k+1)
            enddo
            emin=9999.
            do ii=1,necmx       !check index of c.m. energy
               if(abs(ecm(j)-ecmc(ii)).lt.emin) then
                  iener(nda)=ii
                  emin=abs(ecm(j)-ecmc(ii))
               endif
            enddo
          endif
        endif                   !check 19x ecm range
      enddo                     !loop on c.m. energies
         call mkcov(necmx,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &             ,cov,covec)
         call covinv(nda,necmx,cov,dcovi,w)
         call aver(nda,necmx,wm,cov,w,nerr,err,covec
     &            ,wmav,wmer,e2,0)
         xm_chan19x(i)=wmav
         exm_chan19x(i)=sqrt(e2(1))
         sxm_chan119x(i)=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
         sxm_chan219x(i)=sqrt(e2(14)+e2(15))
         print *,'i m em es channel'
         print *,i,xm_chan19x(i),exm_chan19x(i)
     &        ,sxm_chan119x(i),sxm_chan219x(i)
      enddo                     !loop on channels
      end

      subroutine aver_1_20x
c* average each channel results for 192 to 196 GeV c.m. energies     
      implicit none
      include 'pp.inc'
      integer nda,nerr
      parameter (nerr=ntypsys+1)
      integer iener(necmx)
      real wm(necmx),err(nerr,necmx),rhoerr(nerr)
      real cov(necmx,necmx), w(necmx)
     &     ,covec(necmx,necmx,nerr)
      double precision dcovi(necmx,necmx)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|               AVERAGE each channel 20X               |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,necmx)
      do i=1,nchan
        nda=0
        xm_chan20x(i)=0.
        exm_chan20x(i)=0.
        sxm_chan120x(i)=0.
        sxm_chan220x(i)=0.
        do j=1,necm         !loop over different ecm results for this channel
          if(ecm(j).gt.202..and.ecm(j).lt.209.) then
          emin=9999.
          do ii=1,necmx         !check index of c.m. energy
            if(abs(ecm(j)-ecmc(ii)).lt.emin) then
              iene=ii
              emin=abs(ecm(j)-ecmc(ii))
            endif
          enddo
          if(xmw(i,j).gt.0.) then
            nda=nda+1
            wm(nda)=xmw(i,j)
            err(1,nda)=estat(i,j)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,i,iene)=decmc(iene)*2./ecmc(iene)*80.35
               err(k+1,nda)=esys(k,i,iene)
               rhoerr(k+1)=rho(k)
            enddo
            emin=9999.
            do ii=1,necmx       !check index of c.m. energy
               if(abs(ecm(j)-ecmc(ii)).lt.emin) then
                  iener(nda)=ii
                  emin=abs(ecm(j)-ecmc(ii))
               endif
            enddo
          endif
        endif                   !check 19x ecm range
      enddo                     !loop on c.m. energies
         call mkcov(necmx,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &             ,cov,covec)
         call covinv(nda,necmx,cov,dcovi,w)
         call aver(nda,necmx,wm,cov,w,nerr,err,covec
     &            ,wmav,wmer,e2,0)
         xm_chan20x(i)=wmav
         exm_chan20x(i)=sqrt(e2(1))
         sxm_chan120x(i)=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
         sxm_chan220x(i)=sqrt(e2(14)+e2(15))
         print *,'i m em es channel'
         print *,i,xm_chan20x(i),exm_chan20x(i)
     &        ,sxm_chan120x(i),sxm_chan220x(i)
      enddo                     !loop on channels
      end

      subroutine aver_1
c* average each channel results for all c.m. energies     
      implicit none
      include 'pp.inc'
      integer nda,nerr
      parameter (nerr=ntypsys+1)
      integer iener(necmx)
      real wm(necmx),err(nerr,necmx),rhoerr(nerr)
      real cov(necmx,necmx), w(necmx)
     &     ,covec(necmx,necmx,nerr)
      double precision dcovi(necmx,necmx)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|               AVERAGE each channel                   |' 
      print *,'|-------------------------------------------------------' 
      call vzero(iener,necmx)
      do i=1,nchan
        nda=0
        xm_chan(i)=0.
        exm_chan(i)=0.
        sxm_chan1(i)=0.
        sxm_chan2(i)=0.
        do j=1,necm             !loop over different ecm results for this channel
          emin=9999.
          do ii=1,necmx         !check index of c.m. energy
            if(abs(ecm(j)-ecmc(ii)).lt.emin) then
              iene=ii
              emin=abs(ecm(j)-ecmc(ii))
            endif
          enddo
          if(xmw(i,j).gt.0.) then
            nda=nda+1
            wm(nda)=xmw(i,j)
            err(1,nda)=estat(i,j)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,i,iene)=decmc(iene)*2./ecmc(iene)*80.35
               err(k+1,nda)=esys(k,i,iene)
               rhoerr(k+1)=rho(k)
            enddo
            emin=9999.
            do ii=1,necmx       !check index of c.m. energy
               if(abs(ecm(j)-ecmc(ii)).lt.emin) then
                  iener(nda)=ii
                  emin=abs(ecm(j)-ecmc(ii))
               endif
            enddo
          endif
         enddo                  !loop on c.m. energies
         call mkcov(necmx,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &             ,cov,covec)
         call covinv(nda,necmx,cov,dcovi,w)
         call aver(nda,necmx,wm,cov,w,nerr,err,covec
     &            ,wmav,wmer,e2,0)
         xm_chan(i)=wmav
         exm_chan(i)=sqrt(e2(1))
         sxm_chan1(i)=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
         sxm_chan2(i)=sqrt(e2(14)+e2(15))
         print *,'i m em es channel'
         print *,i,xm_chan(i),exm_chan(i),sxm_chan1(i),sxm_chan2(i)
      enddo                     !loop on channels
      end

      subroutine aver_1e
c* expected error on average of each channel results for all c.m. energies     
      implicit none
      include 'pp.inc'
      integer nda,nerr
      parameter (nerr=ntypsys+1)
      integer iener(necmx)
      real wm(necmx),err(nerr,necmx),rhoerr(nerr)
      real cov(necmx,necmx), w(necmx)
     &     ,covec(necmx,necmx,nerr)
      double precision dcovi(necmx,necmx)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,ndaw,iene
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|            expected AVERAGE per channel              |' 
      print *,'|-------------------------------------------------------' 

      call vzero(iener,necmx)
      do i=1,nchan
        nda=0
        xm_chan_e(i)=0.
        exm_chan_e(i)=0.
        ndaw=0
        do j=1,necm
          if(xmw(i,j).gt.0.)ndaw=ndaw+1
        enddo
        do j=1,necm          !loop over different ecm results for this channel
          emin=9999.
          do ii=1,necmx         !check index of c.m. energy
            if(abs(ecm(j)-ecmc(ii)).lt.emin) then
              iene=ii
              emin=abs(ecm(j)-ecmc(ii))
            endif
          enddo
          if(dmxpect(i,j).gt.0..and.xmw(i,j).gt.0.) then
            nda=nda+1
            wm(nda)=80.35
            err(1,nda)=dmxpect(i,j)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,i,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,i,iene)
              rhoerr(k+1)=rho(k)
            enddo
            emin=9999.
            do ii=1,necmx       !check index of c.m. energy
               if(abs(ecm(j)-ecmc(ii)).lt.emin) then
                  iener(nda)=ii
                  emin=abs(ecm(j)-ecmc(ii))
               endif
            enddo
          endif
         enddo                  !loop on c.m. energies
         if(nda.ge.1.and.nda.eq.ndaw) then
           call mkcov(necmx,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &          ,cov,covec)
           call covinv(nda,necmx,cov,dcovi,w)
           call aver(nda,necmx,wm,cov,w,nerr,err,covec
     &          ,wmav,wmer,e2,0)
           xm_chan_e(i)=wmav
           exm_chan_e(i)=sqrt(e2(1))
           print *,'i m em es channel expected'
           print *,i,xm_chan_e(i),exm_chan_e(i)
         endif
      enddo                     !loop on channels
      end

      subroutine aver_1e19x
c* expected error on average of each channel results for all c.m. energies     
      implicit none
      include 'pp.inc'
      integer nda,nerr
      parameter (nerr=ntypsys+1)
      integer iener(necmx)
      real wm(necmx),err(nerr,necmx),rhoerr(nerr)
      real cov(necmx,necmx), w(necmx)
     &     ,covec(necmx,necmx,nerr)
      double precision dcovi(necmx,necmx)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,ndaw,iene
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|            expected AVERAGE per channel  19x         |' 
      print *,'|-------------------------------------------------------' 

      call vzero(iener,necmx)
      do i=1,nchan
        nda=0
        xm_chan_e19x(i)=0.
        exm_chan_e19x(i)=0.
        ndaw=0
        do j=1,necm
          if(ecm(j).gt.190..and.ecm(j).le.202.) then
            if(xmw(i,j).gt.0.)ndaw=ndaw+1
          endif
        enddo
        do j=1,necm          !loop over different ecm results for this channel
          if(ecm(j).gt.190..and.ecm(j).le.202.) then
          emin=9999.
          do ii=1,necmx         !check index of c.m. energy
            if(abs(ecm(j)-ecmc(ii)).lt.emin) then
              iene=ii
              emin=abs(ecm(j)-ecmc(ii))
            endif
          enddo
          if(dmxpect(i,j).gt.0..and.xmw(i,j).gt.0.) then
            nda=nda+1
            wm(nda)=80.35
            err(1,nda)=dmxpect(i,j)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,i,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,i,iene)
              rhoerr(k+1)=rho(k)
            enddo
            emin=9999.
            do ii=1,necmx       !check index of c.m. energy
               if(abs(ecm(j)-ecmc(ii)).lt.emin) then
                  iener(nda)=ii
                  emin=abs(ecm(j)-ecmc(ii))
               endif
            enddo
          endif
        endif                   !check 19x range
         enddo                  !loop on c.m. energies
         if(nda.ge.1.and.nda.eq.ndaw) then
           call mkcov(necmx,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &          ,cov,covec)
           call covinv(nda,necmx,cov,dcovi,w)
           call aver(nda,necmx,wm,cov,w,nerr,err,covec
     &          ,wmav,wmer,e2,0)
           xm_chan_e19x(i)=wmav
           exm_chan_e19x(i)=sqrt(e2(1))
           print *,'i m em es channel expected'
           print *,i,xm_chan_e19x(i),exm_chan_e19x(i)
         endif
      enddo                     !loop on channels
      end

      subroutine aver_1e20x
c* expected error on average of each channel results for all c.m. energies     
      implicit none
      include 'pp.inc'
      integer nda,nerr
      parameter (nerr=ntypsys+1)
      integer iener(necmx)
      real wm(necmx),err(nerr,necmx),rhoerr(nerr)
      real cov(necmx,necmx), w(necmx)
     &     ,covec(necmx,necmx,nerr)
      double precision dcovi(necmx,necmx)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,ndaw,iene
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|            expected AVERAGE per channel  20x         |' 
      print *,'|-------------------------------------------------------' 

      call vzero(iener,necmx)
      do i=1,nchan
        nda=0
        xm_chan_e20x(i)=0.
        exm_chan_e20x(i)=0.
        ndaw=0
        do j=1,necm
          if(ecm(j).gt.202..and.ecm(j).le.209.) then
            if(xmw(i,j).gt.0.)ndaw=ndaw+1
          endif
        enddo
        do j=1,necm          !loop over different ecm results for this channel
          if(ecm(j).gt.202..and.ecm(j).le.209.) then
          emin=9999.
          do ii=1,necmx         !check index of c.m. energy
            if(abs(ecm(j)-ecmc(ii)).lt.emin) then
              iene=ii
              emin=abs(ecm(j)-ecmc(ii))
            endif
          enddo
          if(dmxpect(i,j).gt.0..and.xmw(i,j).gt.0.) then
            nda=nda+1
            wm(nda)=80.35
            err(1,nda)=dmxpect(i,j)
            rhoerr(1)=0.
            do k=1,ntypsys
              if(k.eq.10) esys(k,i,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,i,iene)
              rhoerr(k+1)=rho(k)
            enddo
            emin=9999.
            do ii=1,necmx       !check index of c.m. energy
               if(abs(ecm(j)-ecmc(ii)).lt.emin) then
                  iener(nda)=ii
                  emin=abs(ecm(j)-ecmc(ii))
               endif
            enddo
          endif
        endif                   !check 19x range
         enddo                  !loop on c.m. energies
         if(nda.ge.1.and.nda.eq.ndaw) then
           call mkcov(necmx,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &          ,cov,covec)
           call covinv(nda,necmx,cov,dcovi,w)
           call aver(nda,necmx,wm,cov,w,nerr,err,covec
     &          ,wmav,wmer,e2,0)
           xm_chan_e20x(i)=wmav
           exm_chan_e20x(i)=sqrt(e2(1))
           print *,'i m em es channel expected'
           print *,i,xm_chan_e20x(i),exm_chan_e20x(i)
         endif
      enddo                     !loop on channels
      end

      subroutine aver_0
c* average lvqq and 4q results for each c.m. energy (and m4q-mlvqq)     
      implicit none
      include 'pp.inc'

      integer ndfit,ndfitmx
      parameter (ndfitmx=36)
      real wmtofit(ndfitmx)
      integer ilq4q(ndfitmx)
      double precision dcovif(ndfitmx,ndfitmx)
      common/forfits/ndfit,ilq4q,wmtofit
      common/forfitd/dcovif

      integer nd,nerr,nda
      parameter (nd=4,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|         AVERAGE mW and m4q-ml at each energy         |' 
      print *,'|-------------------------------------------------------' 

      call vzero(iener,nd)
      do i=1,necm
        emin=9999.
        do ii=1,necmx           !check index of c.m. energy
          if(abs(ecm(i)-ecmc(ii)).lt.emin) then
            iene=ii
            emin=abs(ecm(i)-ecmc(ii))
          endif
        enddo
        nda=0
        xmtot(i)=0.
        emtot(i)=0.
        estot1(i)=0.
        estot2(i)=0.
        do j=1,nd               !loop over lvqq+4q results at this ecm
           es1(j,i)=0.
           es2(j,i)=0.
          if(xmw(j,i).gt.0.) then
            nda=nda+1
            iener(nda)=iene
            wm(nda)=xmw(j,i)
            if(nda.le.ndfitmx) wmtofit(nda)=xmw(j,i)
            err(1,nda)=estat(j,i)
            rhoerr(1)=0.
            do k=1,ntypsys
              rhoerr(k+1)=0.
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,j,iene)
              rhoerr(k+1)=rho(k)
              if(k.le.12) then
                 es1(j,i)=es1(j,i)+esys(k,j,iene)**2
              else
                 es2(j,i)=es2(j,i)+esys(k,j,iene)**2
              endif
            enddo
            es1(j,i)=sqrt(es1(j,i))
            es2(j,i)=sqrt(es2(j,i))
          endif                 !check result available
        enddo
        call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &       ,cov,covec)
        call covinv(nda,nd,cov,dcovi,w)
        call aver(nda,nd,wm,cov,w,nerr,err,covec
     &       ,wmav,wmer,e2,0)
        xmtot(i)=wmav
        emtot(i)=sqrt(e2(1))
        estot1(i)=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
        estot2(i)=sqrt(e2(14)+e2(15))
        print *,'i ecm m em es'
        print *,i,ecm(i),xmtot(i),emtot(i),estot1(i),estot2(i)
        do j=1,nd
          if(xmw(j,i).gt.0.) wei_mtot(j,i)=w(j)
        enddo

        if(nda.eq.4) then
          ndfit=nda
c* switch off FSR eror for mass difference fit
          do j=1,nd
             err(14,j)=0.
             err(15,j)=0.
          enddo
          call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &         ,cov,covec)
          call covinv(nda,nd,cov,dcovi,w)
          do j=1,ndfit
            ilq4q(j)=1
            if(j.eq.4) ilq4q(j)=2
            do k=1,ndfit
              dcovif(j,k)=dcovi(j,k)
            enddo
          enddo
          call wfits(mwdiff(i),emwdiff(i))
        endif

      enddo                     !loop on c.m. energies
      end
      subroutine aver_0a
c* average lvqq and 4q results (first lvqq, then lvqq+4q)     
      implicit none
      include 'pp.inc'

      integer ndfit,ndfitmx
      parameter (ndfitmx=36)
      real wmtofit(ndfitmx)
      integer ilq4q(ndfitmx)
      double precision dcovif(ndfitmx,ndfitmx)
      common/forfits/ndfit,ilq4q,wmtofit
      common/forfitd/dcovif

      integer nd,nerr,nda
      parameter (nd=2,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene
      real emin

      print *,'|-------------------------------------------------------' 
      print *,'|         AVERAGE mW and m4q-ml at each energy (alt)   |' 
      print *,'|-------------------------------------------------------' 

      call vzero(iener,nd)
      do i=1,necm
        emin=9999.
        do ii=1,necmx           !check index of c.m. energy
          if(abs(ecm(i)-ecmc(ii)).lt.emin) then
            iene=ii
            emin=abs(ecm(i)-ecmc(ii))
          endif
        enddo
        nda=0
        xmtot(i)=0.
        emtot(i)=0.
        estot1(i)=0.
        estot2(i)=0.
        do j=1,2               !loop over lvqq+4q results at this ecm
          es1(j,i)=0.
          es2(j,i)=0.
        enddo
        if(xmlvqq(i).gt.0.) then
          nda=nda+1
          iener(nda)=iene
          wm(nda)=xmlvqq(i)
          if(nda.le.ndfitmx) wmtofit(nda)=xmlvqq(i)
          err(1,nda)=emlvqq(i)
          do k=1,ntypsys
            err(k+1,nda)=es_lq(k,iene)
          enddo
        endif
        if(xmw(4,i).gt.0.) then
          nda=nda+1
          iener(nda)=iene
          wm(nda)=xmw(4,i)
          if(nda.le.ndfitmx) wmtofit(nda)=xmw(4,i)
          err(1,nda)=estat(4,i)
          rhoerr(1)=0.
          do k=1,ntypsys
            if(k.eq.10) esys(k,4,iene)=decmc(iene)*2./ecmc(iene)*80.35
            err(k+1,nda)=esys(k,4,iene)
            rhoerr(k+1)=rho(k)
            if(k.le.12) then
              es1(4,i)=es1(4,i)+esys(k,4,iene)**2
            else
              es2(4,i)=es2(4,i)+esys(k,4,iene)**2
            endif
          enddo
          es1(4,i)=sqrt(es1(4,i))
          es2(4,i)=sqrt(es2(4,i))
        endif                   !check result available
        call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &       ,cov,covec)
        call covinv(nda,nd,cov,dcovi,w)
        call aver(nda,nd,wm,cov,w,nerr,err,covec
     &       ,wmav,wmer,e2,0)
        xmtot(i)=wmav
        emtot(i)=sqrt(e2(1))
        estot1(i)=sqrt(wmer**2-e2(1)-e2(14)-e2(15))
        estot2(i)=sqrt(e2(14)+e2(15))
        print *,'i ecm m em es'
        print *,i,ecm(i),xmtot(i),emtot(i),estot1(i),estot2(i)

        if(nda.eq.2) then
          ndfit=nda
c* switch off FSR eror for mass difference fit
          do j=1,nd
             err(14,j)=0.
             err(15,j)=0.
          enddo
          call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &         ,cov,covec)
          call covinv(nda,nd,cov,dcovi,w)
          do j=1,ndfit
            ilq4q(j)=1
            if(j.eq.2) ilq4q(j)=2
            do k=1,ndfit
              dcovif(j,k)=dcovi(j,k)
            enddo
          enddo
          call wfits(mwdiff(i),emwdiff(i))
        endif

      enddo                     !loop on c.m. energies
      end

      subroutine aver_0e
c* expected error on average of lvqq and 4q results for each c.m. energy     
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=4,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr)
      integer i,j,ii,k,iene
      real emin

      integer ndfit,ndfitmx
      parameter (ndfitmx=36)
      real wmtofit(ndfitmx)
      integer ilq4q(ndfitmx)
      double precision dcovif(ndfitmx,ndfitmx)
      common/forfits/ndfit,ilq4q,wmtofit
      common/forfitd/dcovif

      print *,'|-------------------------------------------------------' 
      print *,'|     expected AVERAGE mW and m4q-ml at each energy    |' 
      print *,'|-------------------------------------------------------' 

      call vzero(iener,nd)
      do i=1,necm
        emin=9999.
        do ii=1,necmx           !check index of c.m. energy
          if(abs(ecm(i)-ecmc(ii)).lt.emin) then
            iene=ii
            emin=abs(ecm(i)-ecmc(ii))
          endif
        enddo
        nda=0
        xmtot_e(i)=0.
        emtot_e(i)=0.
        do j=1,nd               !loop over lvqq+4q results at this ecm
          if(dmxpect(j,i).gt.0.) then
            nda=nda+1
            iener(nda)=iene
            wm(nda)=80.35
            if(nda.le.ndfitmx) wmtofit(nda)=80.35
            err(1,nda)=dmxpect(j,i)
            rhoerr(1)=0.
            do k=1,ntypsys
              rhoerr(k+1)=0.
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
              err(k+1,nda)=esys(k,j,iene)
              rhoerr(k+1)=rho(k)
            enddo
          endif                 !check result available
        enddo
        if(nda.eq.4) then
          call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &         ,cov,covec)
          call covinv(nda,nd,cov,dcovi,w)
         call aver(nda,nd,wm,cov,w,nerr,err,covec
     &         ,wmav,wmer,e2,0)
          xmtot_e(i)=wmav
          emtot_e(i)=sqrt(e2(1))
        endif

        if(nda.eq.4) then
          ndfit=nda
c* switch off FSR eror for mass difference fit
          do j=1,nd
            err(14,j)=0.
            err(15,j)=0.
          enddo
          call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &         ,cov,covec)
          call covinv(nda,nd,cov,dcovi,w)
          do j=1,ndfit
            ilq4q(j)=1
            if(j.eq.4) ilq4q(j)=2
            do k=1,ndfit
              dcovif(j,k)=dcovi(j,k)
            enddo
          enddo
         call aver(nda,nd,wm,cov,w,nerr,err,covec
     &         ,wmav,wmer,e2,0)
          call wfits(mwdiff_e(i),emwdiff_e(i))
        endif
      enddo                     !loop on c.m. energies
      end

      subroutine aver_lvqq
c* average evqq, muvqq and tauvqq result at each c.m. energy
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=3,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr),emin
      integer i,j,ii,k,iene,iprint

      print *,'-------------------------------------------------------' 
      print *,'|               AVERAGE mlvqq at each energy          |' 
      print *,'-------------------------------------------------------' 

      call vzero(iener,nd)
      do i=1,necm
         iprint=0
         emin=9999.
         do ii=1,necmx          !check index of c.m. energy
            if(abs(ecm(i)-ecmc(ii)).lt.emin) then
               iene=ii
               emin=abs(ecm(i)-ecmc(ii))
            endif
         enddo
         nda=0
         xmlvqq(i)=0.
         emlvqq(i)=0.
         eslvqq(i)=0.
         do j=1,nd              !loop over lvqq results at this ecm
           if(xmw(j,i).gt.0.) then
             nda=nda+1
             iener(nda)=iene
             wm(nda)=xmw(j,i)
             err(1,nda)=estat(j,i)
             rhoerr(1)=0.
             es1(j,i)=0.
             do k=1,ntypsys
               rhoerr(k+1)=0.
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
               err(k+1,nda)=esys(k,j,iene)
               rhoerr(k+1)=rho(k)
               es1(j,i)=es1(j,i)+esys(k,j,iene)**2
             enddo
             es1(j,i)=sqrt(es1(j,i)) !total systematic for this channel
           endif                !check result available
         enddo
         call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &             ,cov,covec)
         call covinv(nda,nd,cov,dcovi,w)
         call aver(nda,nd,wm,cov,w,nerr,err,covec
     &            ,wmav,wmer,e2,iprint)
         xmlvqq(i)=wmav
         emlvqq(i)=sqrt(e2(1))
         eslvqq(i)=sqrt(wmer**2-e2(1))
c store syst for overall lvqq 
         do j=1,ntypsys
            es_lq(j,iene)=sqrt(e2(j+1))
         enddo
         print *,'i ecm m em es lvqq'
         print *,i,ecm(i),xmlvqq(i),emlvqq(i),eslvqq(i)
      enddo                     !loop on c.m. energies
      end

      subroutine xpect_lvqq
c* expected error on average lvqq result at each c.m. energy
      implicit none
      include 'pp.inc'
      integer nd,nerr,nda
      parameter (nd=3,nerr=ntypsys+1)
      integer iener(nd)
      real wm(nd),err(nerr,nd),rhoerr(nerr)
      real cov(nd,nd), w(nd),covec(nd,nd,nerr)
      double precision dcovi(nd,nd)
      real wmav,wmer,e2(nerr),emin
      integer i,j,ii,k,iene,iprint

      print *,'-------------------------------------------------------' 
      print *,'|     expected  AVERAGE mlvqq at each energy          |' 
      print *,'-------------------------------------------------------' 

      call vzero(iener,nd)
      do i=1,necm
         iprint=0
         emin=9999.
         do ii=1,necmx          !check index of c.m. energy
            if(abs(ecm(i)-ecmc(ii)).lt.emin) then
               iene=ii
               emin=abs(ecm(i)-ecmc(ii))
            endif
         enddo
         nda=0
         xmlvqq_e(i)=0.
         emlvqq_e(i)=0.
         eslvqq_e(i)=0.
         do j=1,nd              !loop over lvqq results at this ecm
           if(dmxpect(j,i).gt.0.) then
             nda=nda+1
             iener(nda)=iene
             wm(nda)=80.35
             err(1,nda)=dmxpect(j,i)
             rhoerr(1)=0.
             do k=1,ntypsys
               rhoerr(k+1)=0.
              if(k.eq.10) esys(k,j,iene)=decmc(iene)*2./ecmc(iene)*80.35
               err(k+1,j)=esys(k,j,iene)
               rhoerr(k+1)=rho(k)
             enddo
           endif                !check result available
         enddo
         if(nda.eq.3) then
           call mkcov(nd,wm,nerr,err,rhoerr,necm,rhoecm,iener
     &          ,cov,covec)
           call covinv(nda,nd,cov,dcovi,w)
           call aver(nda,nd,wm,cov,w,nerr,err,covec
     &          ,wmav,wmer,e2,iprint)
           xmlvqq_e(i)=wmav
           emlvqq_e(i)=sqrt(e2(1))
           eslvqq_e(i)=sqrt(wmer**2-e2(1))
           print *,'i ecm m em es lvqq expected'
           print *,i,ecm(i),xmlvqq_e(i),emlvqq_e(i),eslvqq_e(i)
         endif
      enddo                     !loop on c.m. energies
      end

      SUBROUTINE mkcov(nd,wmass,nerr,err,rhoerr,necm,rhoecm,iener
     &                ,cov,covec)
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*... MaKe COVariance and error correlation matrices
*********************************************************************
      IMPLICIT none
      integer nd,nerr,necm
      integer iener(nd)
      real wmass(nd),err(nerr,nd),rhoerr(nerr)
      real rhoecm(necm,necm)
C
      real COV(Nd,nd),COVEC(nd,nd,nerr)
      INTEGER IErr,I,J
      REAL R
C
C -----------------------------------------------------------------
C
      Do i=1,nd
         Do j=1,nd
            cov(i,j)=0.
            if(iener(j).gt.0) then
            Do ierr=1,nerr
               R=1.
               if(i.ne.j) R=rhoerr(ierr)
               if(abs(r-3.).lt.0.001) r=rhoecm(iener(i),iener(j))
               covec(i,j,ierr)=R*err(ierr,i)*err(ierr,j)
               cov(i,j)=cov(i,j)+covec(i,j,ierr)
            ENDDO
            endif
         ENDDO                  ! second index
      ENDDO                     ! first index
C
      END

      SUBROUTINE COVINV(nd,ndimx,cov,dcovi,w)
      IMPLICIT none
      integer nd,ndimx
      REAL COV(ndimx,ndimx), COVI(ndimx,ndimx) ,W(ndimx)
C
      DOUBLE PRECISION DCOVI(ndimx,ndimx)
      DOUBLE PRECISION WORK(ndimx)
      INTEGER I,J,K,L,Ndim,Ierr
      REAL ANORM
C
C --------------------------------------------------------------------
C
      NDIM=ND
C
      Do I=1,NDIM
        Do J=1,NDIM
          DCOVI(I,J)=COV(I,J)
        ENDDO
      ENDDO
C
      CALL DINV(NDIM,DCOVI,Ndimx,WORK,IERR)
      IF(IERR.NE.0) STOP '*==* COVINV Inversion of COV matrix failed'
C
C  compute weights
      ANORM=0.
      Do K=1,NDIM
       Do L=1,NDIM
        ANORM=ANORM+DCOVI(K,L)
       ENDDO
      ENDDO
      DO I=1,NDIM
       W(I)=0.
       DO J=1,NDIM
        W(I)=W(I)+DCOVI(I,J)
       ENDDO
       W(I)=W(I)/Anorm
      ENDDO
C
      RETURN
      END


      subroutine aver(nd,ndmx,wm,cov,w,nerr,err,covec
     &               ,wmav,wmer,e2,iprint)
      implicit none
      integer nd,nerr,ndmx
      real wm(ndmx),cov(ndmx,ndmx),w(ndmx),err(nerr,ndmx)
      real covec(ndmx,ndmx,nerr)
      real wmav,wmer,e2(nerr)
      real etotsys
      integer i,j,ierr,iprint

      WMAv=0.
      WmEr=0.
      CALL VZERO(E2,Nerr)
C
      DO I=1,ND
C     calculate average
         WMAv=WMAv+W(I)*WM(I)
         DO J=1,ND
C     calculate total error ...
            WMEr=WMEr+W(I)*W(J)*COV(I,J)
C     ... and error break-down
            DO Ierr=1,Nerr
               E2(Ierr)=E2(Ierr)
     &              +W(I)*W(J)*COVEC(I,J,Ierr)
            ENDDO
         ENDDO
      ENDDO
      WMEr=SQRT(WMEr)
C     
      if(iprint.eq.1) then
        PRINT*,' '
        PRINT*,
     &     '*==* AveMW: +++++++++ analytical average ++++++++++++++++++'
C     
        PRINT('(A)'),
     &       '     weights of channels'
        PRINT('(4X,6X,5X,4(1X,F4.2))'),
     &       (W(I),I=1,Nd)
C     
        PRINT*,
     &       '     average W mass= ',WMAv,'+/-',WMEr
        etotsys=0.
        do i=2,nerr
          etotsys=etotsys+e2(i)
        enddo
        etotsys=sqrt(etotsys)
        PRINT*,
     &       '    -----------------------------------------------'
        PRINT*,'    error break-down'
        PRINT*,'        stat. ',SQRT(E2(1))
        print*,'      totsys. ',etotsys
        PRINT*,'        char. ',SQRT(E2(2))
        PRINT*,'        lang. ',SQRT(E2(3))
        PRINT*,'        lmom. ',SQRT(E2(4))
        PRINT*,'        jcor. ',SQRT(E2(5))
        PRINT*,'        jres. ',SQRT(E2(6))
        PRINT*,'        calo. ',SQRT(E2(7))
        PRINT*,'        jang. ',SQRT(E2(8))
        PRINT*,'       fragm. ',SQRT(E2(9))
        PRINT*,'         rad. ',SQRT(E2(10))
        PRINT*,'         lep. ',SQRT(E2(11))
        PRINT*,'          mc. ',SQRT(E2(12))
        PRINT*,'        back. ',SQRT(E2(13))
        PRINT*,'          CR. ',SQRT(E2(14))
        PRINT*,'          BE. ',SQRT(E2(15))
        PRINT*,'       BE&CR  ',SQRT(E2(14)+E2(15))
        PRINT*,'      Detector',
     &       SQRT(E2(2)+E2(3)+E2(4)+E2(5)+E2(6)+E2(7)+E2(8))
        PRINT*,' '
      endif
C
      RETURN
      END
      SUBROUTINE WFITS(mwdif,emwdif)
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*... WFITs    : Fit W mass(s) with MINUIT
*********************************************************************
      IMPLICIT none
      INTEGER LUNI,LUNO,LUNS
      REAL MMW,MMWlq,MMW4q,MClq4q,MMWDIF
     &   ,MEMW,MEMWlq,MEMW4q,MEMWDIF,c2MW2,c2MW
c      COMMON /MINFIT/MMW,MMWlq,MMW4q,MClq4q,MMWDIF
c     &   ,MEMW,MEMWlq,MEMW4q,MEMWDIF,c2MW2,c2MW
C
* variables for call to MNSTAT
      DOUBLE PRECISION FMIN,FEDM,ERRDEF
      INTEGER NPARI,NPARX,ISTAT
* variables for call to MNPOUT
      CHARACTER*10 CHNAM
      DOUBLE PRECISION VAL,ERR,B1,B2
      INTEGER ISTV
* variables for MNEXCM
      INTEGER IERR,NARG
      DOUBLE PRECISION  ARG(10)
      CHARACTER*10 CMD
* variable for MNSETI
      CHARACTER*50 CTITLE
* variables for MNPARS
      CHARACTER*80 CHpar
      INTEGER Icond
C
      DOUBLE PRECISION CMAT(2,2)
C
      real mwdif, emwdif
      EXTERNAL WFCN
C
C ----------------------------------------------------------------------
C
C initialise output variables
      MMW=0.
      MMWlq=0.
      MMW4q=0.
      MClq4q=0.
      MMWDIF=0.
      MEMW=0.
      MEMWlq=0.
      MEMW4q=0.
      MEMWDIF=0.
C
*  initialize MINUIT
      LUNI=5
      LUNO=6
      LUNS=10
      CALL MNINIT(LUNI,LUNO,LUNS)
C
C first is a fit for the W mass
      CTITLE='W mass fit'
      CALL MNSETI(CTITLE)
      WRITE(CHPAR,'(A)') '1 ''MW'' 80.4 0.050'
      CALL MNPARS(CHpar,Icond)
      IF(Icond.NE.0) THEN
       PRINT*, '*==* WFITS: error from MNPARS ',Icond
       STOP 'error in WFITS'
      ENDIF
C
* minimisation
      NARG=2
      ARG(1)=1000.D0
      ARG(2)=.1D0
      CMD='MINIMIZE'
      CALL MNEXCM(WFCN,CMD,ARG,NARG,IERR,0)
      IF(IERR.NE.0) STOP '!!! failed to MINIMIZE'
C get chi2 from Minuit
      CALL MNSTAT(FMIN,FEDM,ERRDEF,NPARI,NPARX,ISTAT)
      c2MW=FMIN
c get results
      CALL MNPOUT(1,CHNAM,VAL,ERR,B1,B2,ISTV)
      MMW=VAL
      MEMW=ERR
C
* now come the correlated MWql,MWqq masses
      CTITLE='MWlq, MWqq  fit'
      CALL MNSETI(CTITLE)
      WRITE(CHPAR,'(A)') '1 ''MWlq'' 80.4 0.050'
      CALL MNPARS(CHpar,Icond)
      IF(Icond.NE.0) THEN
       PRINT*, '*==* WFITS: error from MNPARS ',Icond
       STOP 'error in WFITS'
      ENDIF
      WRITE(CHPAR,'(A)') '2 ''MWqq'' 80.4 0.050'
      CALL MNPARS(CHpar,Icond)
      IF(Icond.NE.0) THEN
       PRINT*, '*==* WFITS: error from MNPARS ',Icond
       STOP 'error in WFITS'
      ENDIF
* minimisation
      NARG=2
      ARG(1)=1000.D0
      ARG(2)=.1D0
      CMD='MINIMIZE'
      CALL MNEXCM(WFCN,CMD,ARG,NARG,IERR,0)
      IF(IERR.NE.0) STOP '!!! failed to MINIMIZE'
C get chi2 from Minuit
      CALL MNSTAT(FMIN,FEDM,ERRDEF,NPARI,NPARX,ISTAT)
      c2MW2=FMIN
c get results
      CALL MNPOUT(1,CHNAM,VAL,ERR,B1,B2,ISTV)
      MMWlq=VAL
      MEMWlq=ERR
      CALL MNPOUT(2,CHNAM,VAL,ERR,B1,B2,ISTV)
      MMW4q=VAL
      MEMW4q=ERR
      CALL MNEMAT(CMAT,2)
      MClq4q=CMAT(1,2)/MEMWlq/MEMW4q
C
* MWql and mass difference
      CTITLE='MWlq, MWqq  fit'
      CALL MNSETI(CTITLE)
      WRITE(CHPAR,'(A)') '1 ''MWlq'' 80.4 0.050'
      CALL MNPARS(CHpar,Icond)
      IF(Icond.NE.0) THEN
       PRINT*, '*==* WFITS: error from MNPARS ',Icond
       STOP 'error in WFITS'
      ENDIF
      WRITE(CHPAR,'(A)') '2 ''MWdiff'' 0. 0.050'
      CALL MNPARS(CHpar,Icond)
      IF(Icond.NE.0) THEN
       PRINT*, '*==* WFITS: error from MNPARS ',Icond
       STOP 'error in WFITS'
      ENDIF
* minimisation
      NARG=2
      ARG(1)=1000.D0
      ARG(2)=.1D0
      CMD='MINIMIZE'
      CALL MNEXCM(WFCN,CMD,ARG,NARG,IERR,0)
      IF(IERR.NE.0) STOP '!!! failed to MINIMIZE'
      CALL MNPOUT(2,CHNAM,VAL,ERR,B1,B2,ISTV)
      MMWDIF=VAL
      MEMWDIF=ERR
      mwdif=mmwdif
      emwdif=memwdif
C
  999 CONTINUE
      RETURN
      END
      SUBROUTINE WFCN(NPAR,GRAD,FVAL,XVAL,IFLAG,FUTIL)
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*... WFCN     : MINUIT interface
*                DCOVI in common /DWCOVI/ from S/R COVINV
*********************************************************************
      IMPLICIT none
      INTEGER NPAR,IFLAG,FUTIL
      DOUBLE PRECISION GRAD(*),XVAL(*),FVAL
* variables for call to MNPOUT
      CHARACTER*10 CHNAM
      DOUBLE PRECISION VAL,ERR,B1,B2
      INTEGER ISTV
C
      INTEGER I,J
      DOUBLE PRECISION MW(2)

      integer ndfit,ndfitmx
      parameter (ndfitmx=36)
      real wmtofit(ndfitmx)
      integer ilq4q(ndfitmx)
      double precision dcovif(ndfitmx,ndfitmx)
      common/forfits/ndfit,ilq4q,wmtofit
      common/forfitd/dcovif
C
C ----------------------------------------------------------------------
C
      FVAL=0.D0
C
      IF(IFLAG.EQ.1) THEN
C INITIALIZATION...
C    ... nothing to do here
       RETURN
      ENDIF
*
      IF(IFLAG.EQ.2) THEN
       PRINT *,
     & '*==* ZFCN Computation of derivatives not implemented'
C OPTIONAL COMPUTATION OF DERIVATIVES
      ENDIF
*
      IF(IFLAG.EQ.3) THEN
C END
       RETURN
      ENDIF
C
* compute CHI2
      IF(NPAR.EQ.2) THEN
       MW(1)=XVAL(1)  ! qqlv mass in any case
       CALL MNPOUT(2,CHNAM,VAL,ERR,B1,B2,ISTV)
       IF(CHNAM.EQ.'MWqq') THEN
        MW(2)=XVAL(2) ! qqqq mass
       ELSE IF(CHNAM.EQ.'MWdiff') THEN
        MW(2)=MW(1)+XVAL(2) ! mass difference
       ELSE
        PRINT*,'*==* WFCN unknown variable to fit for '//CHNAM
        STOP 'error in WFCN'
       ENDIF
      ELSE IF(NPAR.EQ.1) THEN
       MW(1)=XVAL(1)
       MW(2)=XVAL(1)
      ELSE
       PRINT*, '*==* WFCN: unexpected # of parameters ',NPAR
       STOP 'error in WFCN'
      ENDIF

      DO I=1,Ndfit
       Do J=1,Ndfit
        FVAL=FVAL +
     &     (WMtofit(I)-MW(Ilq4q(I))) * DCOVIf(I,J)
     &                          * (WMtofit(J)-MW(Ilq4q(J)))
       ENDDO
      ENDDO
C
      RETURN
      END

      subroutine writweb
      implicit none
      include 'pp.inc'
      integer elearn(nchan,necmx)
      real cutlearn(nchan,necmx)
      data elearn / 4*0                !172 GeV
     &            , 4*189              !183
     &            , 4*189              !189 
     &            , 3*196, 189         !192
     &            , 4*196              !196
     &            , 3*200, 196         !200
     &            , 3*200, 196         !202
     &            , 4*207              !205
     &            , 4*207  /           !207
      data cutlearn / 0.6, 0.6, 0.4, 0.3
     &              , 0.6, 0.6, 0.4, 0.3
     &              , 0.6, 0.6, 0.4, 0.3
     &              , 0.6, 0.6, 0.4, 0.3
     &              , 0.6, 0.6, 0.4, 0.3
     &              , 0.6, 0.6, 0.4, 0.3
     &              , 0.6, 0.6, 0.3, 0.3
     &              , 0.6, 0.6, 0.3, 0.3
     &              , 0.6, 0.6, 0.4, 0.3 /
      character*80 chw(nchan), titr, tit_align, weblink
      data chw / 'e', '&micro;', '<font FACE="Symbol">t</font>', '4q'/ 
      character*8 cqdate, cqtime
      integer jdat, jtim
      DIMENSION JDAT(2), JTIM(2)
      integer lweb
      real sumlum,sumlum19x,sumlum20x
      integer n1tot(nchan),n2tot(nchan)
      integer n119x(nchan),n219x(nchan),n120x(nchan),n220x(nchan)
      integer i,j,iecm,i1,i2,k,nn1,nn2,ic,ii,necmtot
      real x,ex,sx1,sx2,dmx,dmx2
      real esl,est(5),poids
      character*40 chsys(ntypsys)
      data chsys/'Charged tracking'
     &          ,'e+mu angle bias/resolution'
     &          ,'e+mu momentum bias/resolution'  
     &          ,'Jet energy corrections'
     &          ,'Jet energy resolution'
     &          ,'Calorimeter simulation'
     &          ,'Jet angle bias/resolution'
     &          ,'Fragmentation'
     &          ,'Missing radiative corrections'
     &          ,'LEP energy'
     &          ,'Reference MC statistics'
     &          ,'Background contamination'
     &          ,'Colour reconnection'
     &          ,'Bose-Einstein effects'  /

      lweb=8
      OPEN(Lweb,FILE='mass-test.html',STATUS='unknown')

      CALL DATIMH (JDAT,JTIM)
      WRITE(CQDATE(1:4),'(A4)') JDAT(1)
      WRITE(CQDATE(5:8),'(A4)') JDAT(2)
      WRITE(CQTIME(1:4),'(A4)') JTIM(1)
      WRITE(CQTIME(5:8),'(A4)') JTIM(2)

      write(lweb,*)'<html>'
      write(lweb,*)
     &     '<body bgcolor="#F0FFF0" LINK="#0000ff" VLINK="#800080">'
      write(lweb,*)'<b><font SIZE="6" COLOR="#ff0000">'
      write(lweb,*)'<i><p>Preliminary </i></font>'
      write(lweb,*)'<font SIZE="6">Mass Results' 
      write(lweb,*)'(<i>Last update ',cqdate,' </i>)' 
      write(lweb,*)'</font></b></p>'
      write(lweb,*)' '
      write(lweb,*)'<hr>'
      write(lweb,*)' '
      write(lweb,*)'<p>&nbsp;</p>'
      write(lweb,*)'<p>this page is created from a'
      write(lweb,*)'<a href="wm_inputs.cards">card file</a>,'
      write(lweb,*)'a fortran <a href="pp2.f">program</a>'
      write(lweb,*)' and an <a href="pp.inc">include file</a><br>'
      write(lweb,*)
     &' (which you will want to improve if you look at it!)</p>'
      write(lweb,*)' '
      write(lweb,*)'<p><img SRC="http://alephwww.cern.ch/LEP2ANALYSIS/'
      write(lweb,*)
     &'WTF/Picts/Bullet2.gif" WIDTH="15" HEIGHT="12">&nbsp;'
      write(lweb,*)'<a HREF="#General">General</a></p>'
      write(lweb,*)' '
      write(lweb,*)'<p><img SRC="http://alephwww.cern.ch/LEP2ANALYSIS/'
      write(lweb,*)
     &'WTF/Picts/Bullet2.gif" WIDTH="15" HEIGHT="12">&nbsp;'
      write(lweb,*)'<a HREF="#Results">Results</a></p>'
      write(lweb,*)' '
      write(lweb,*)'<p><img SRC="http://alephwww.cern.ch/LEP2ANALYSIS/'
      write(lweb,*)
     &'WTF/Picts/Bullet2.gif" WIDTH="15" HEIGHT="12">&nbsp;'
      write(lweb,*)'<a HREF="#Systematics">Systematics</a></p>'
      write(lweb,*)' '
      write(lweb,*)
     &     '<h2>Results obtained with:<a NAME="General"></a></h2>'
      write(lweb,*)' '
      write(lweb,*)'<ul>'
      write(lweb,*)'  <li>GAL version xxx </li>'
      write(lweb,*)'  <li>JUL version xxx </li>'
      write(lweb,*)'  <li>MIN version 202 </li>'
      write(lweb,*)'  <li>DBA version 307 </li>'
      write(lweb,*)'  <li>ntuple maker version 1.2 </li>'
      write(lweb,*)'  <li>selections were made using </li>'
      write(lweb,*)'</ul>'
      write(lweb,*)' '

      write(lweb,*)'<table BORDER="1" CELLSPACING="1" WIDTH="618">'
      write(lweb,*)'  <tr>'
      do i=1,4
        write(lweb,*)
     &     '    <td WIDTH="11%" VALIGN="MIDDLE" COLSPAN="2">&nbsp;</td>'
        write(lweb,*)
     &       '    <td WIDTH="9%" VALIGN="MIDDLE"><i><p ALIGN="CENTER">'
        write(lweb,*)'NN trained at</i> </td>'
        write(lweb,*)
     &       '    <td WIDTH="5%" VALIGN="MIDDLE"><i><p ALIGN="CENTER">'
        write(lweb,*)'cut at</i> </td>'
      enddo
      write(lweb,*)'  </tr>'
      write(lweb,*)'  <tr>'

      do k=1,2
        i1=2
        i2=5
        if(k.eq.2) i1=6
        if(k.eq.2) i2=9
 1      do i=i1,i2
      write(lweb,*)'    <td WIDTH="8%" VALIGN="MIDDLE" ROWSPAN="4">'
      write(lweb,'(A,f7.3,A)')
     &   '    <font COLOR="#ff0000"><p ALIGN="CENTER">',ecm(i),'</font>'
      write(lweb,*)'  </td>'
      write(lweb,*)'  <td WIDTH="3%" VALIGN="MIDDLE">',chw(1)(1:1)
     &     ,'</td>'
      write(lweb,*)'  <td WIDTH="9%" VALIGN="MIDDLE">'
     &       ,elearn(1,i),'</td>'
      write(lweb,'(A,f4.2,A)')'    <td WIDTH="5%" VALIGN="MIDDLE">'
     &       ,cutlearn(1,i),'</td>'
      enddo

      write(lweb,*)'  </tr>'
      write(lweb,*)'  <tr>'

      do i=i1,i2
      write(lweb,*)'    <td WIDTH="3%" VALIGN="MIDDLE">',chw(2)(1:7)
     &       ,'</td>'
      write(lweb,*)'  <td WIDTH="9%" VALIGN="MIDDLE">'
     &       ,elearn(2,i),'</td>'
      write(lweb,'(A,f4.2,A)')'    <td WIDTH="5%" VALIGN="MIDDLE">'
     &       ,cutlearn(2,i),'</td>'
      enddo

      write(lweb,*)'  </tr>'
      write(lweb,*)'  <tr>'

      do i=i1,i2
      write(lweb,*)'    <td WIDTH="3%" VALIGN="MIDDLE">',chw(3)(1:28)
     &     ,'</td>'
      write(lweb,*)'  <td WIDTH="9%" VALIGN="MIDDLE">'
     &       ,elearn(3,i),'</td>'
      write(lweb,'(A,f4.2,A)')'    <td WIDTH="5%" VALIGN="MIDDLE">'
     &       ,cutlearn(3,i),'</td>'
      enddo

      write(lweb,*)'  </tr>'
      write(lweb,*)'  <tr>'

      do i=i1,i2
      write(lweb,*)'    <td WIDTH="3%" VALIGN="MIDDLE">',chw(4)(1:2)
     &     ,'</td>'
      write(lweb,*)'  <td WIDTH="9%" VALIGN="MIDDLE">'
     &       ,elearn(4,i),'</td>'
      write(lweb,'(A,f4.2,A)')'    <td WIDTH="5%" VALIGN="MIDDLE">'
     &       ,cutlearn(4,i),'</td>'
      enddo

      write(lweb,*)'  </tr>'
      write(lweb,*)'  <tr>'

      enddo

      write(lweb,*)'</table>'
      write(lweb,*)'<hr>'
      write(lweb,*)' '
      write(lweb,*)
     &     '<h2><big>Table of Results</big><a NAME="Results"></a></h2>'
      write(lweb,*)' '
      write(lweb,*)'<ul>'
      write(lweb,*)'  <li>N<sub>1</sub>'
      write(lweb,*)
     &'&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'
      write(lweb,*)'   =  number of events selected by NN </li>'
      write(lweb,*)'  <li>N<sub>2</sub>'
      write(lweb,*)
     &'&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'
      write(lweb,*)
     &'   =  number of events selected for mass measurement </li>'
      write(lweb,*)'  <font FACE="Symbol">'
      write(lweb,*)'  <li>D</font>M'
      write(lweb,*)'&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'
      write(lweb,*)'   =  expected statistical error (from samples)'
      write(lweb,*)' in GeV/c<sup>2</sup></li>'
      write(lweb,*)'  <li><font FACE="Symbol">D</font>M<sub>200</sub>'
      write(lweb,*)'&nbsp;&nbsp;&nbsp;'
      write(lweb,*)
     &'   =  expected statistical error (for an integrated luminosity'
      write(lweb,*)' of 200 pb<sup>-1</sup>) in GeV/c<sup>2</sup></li>'
      write(lweb,*)'  <li><sup><big>For the Data, the first error is '
      write(lweb,*)
     &'statistical, then systematic, including the error due to LEP'
      write(lweb,*)' energy</big></sup>. '
      write(lweb,*)'  <sup><big>For the 4q </big></sup><br>'
      write(lweb,*)'  <sup><big>channel, an additional error comes from'
      write(lweb,*)
     &'FSI, which is propagated as the last error on the total of all'
      write(lweb,*)' channels.</big></sup></li>'
      write(lweb,*)'</ul>'
      write(lweb,*)' '
      write(lweb,*)
     &'<p>NB:<br> Present mass numbers for Ecm=182.65 and 188.6 GeV'
      write(lweb,*)' are from the published papers.<br>'
      write(lweb,*)'The other numbers are given'
      write(lweb,*)'by D. Boumediene, G. Leibenguth, A. Moutoussi'
      write(lweb,*)'and J. Nowell.<br>'
      write(lweb,*)'<font FACE="Symbol">D</font>M'
      write(lweb,*)'is normalized to the observed number of events for'
      write(lweb,*)' published numbers and normalized to the'
      write(lweb,*)'theoretically expected<br>' 
      write(lweb,*)'number of events for the'
      write(lweb,*)'unpublished. For instance the expected error in the' 
      write(lweb,*)'4q channel at 189 GeV is 96 MeV/c<sup>2</sup><br>' 
      write(lweb,*)'(90 MeV/c<sup>2</sup> if scaled to'
      write(lweb,*)'200 pb<sup>-1</sup>), but due to the'
      write(lweb,*)'low number of observed events, the expected error'
      write(lweb,*)'is rescaled to 101 MeV/c<sup>2</sup>.<br>'
      write(lweb,*)'Some cells are clickable for more details. <br>'
      write(lweb,*)'Systematic errors '
      write(lweb,*)'are taken from the tables below.<br>'
      write(lweb,*)'The error on M<SMALL>4q</SMALL>-M<SMALL>'
     &          //'<FONTface=Vivante-DTC>l</FONT></SMALL>'
      write(lweb,*)'is the overall stat+syst error'
      write(lweb,*)'(same for the expectation)'        

c-------------------------------------------------------------------------
      write(lweb,*)'<table BORDER="1" CELLSPACING="1" WIDTH="732" '
      write(lweb,*)'height="1476" bordercolor="#000000">'
      write(lweb,*)'  <tr>'
      write(lweb,*)'    <td WIDTH="50" VALIGN="middle" height="25"'
      write(lweb,*)'    style="border-left: 3px solid rgb(0,128,0); '
      write(lweb,*)'border-right: 1px solid rgb(0,0,0); '
      write(lweb,*)'border-top: 3px solid rgb(0,128,0); '
      write(lweb,*)'border-bottom: 3px solid rgb(0,128,0)">'
      write(lweb,*)'<p ALIGN="CENTER">Ecm</td>'

      write(lweb,*)'    <td WIDTH="50" VALIGN="middle" height="25"'
      write(lweb,*)'    style="border-left: 1px solid rgb(0,0,0); '
      write(lweb,*)'border-right: 1px solid rgb(0,0,0); ' 
      write(lweb,*)'border-top: 3px solid rgb(0,128,0); '
      write(lweb,*)'border-bottom: 3px solid rgb(0,128,0)">'
      write(lweb,*)'<p ALIGN="CENTER">L(pb<sup>-1</sup>)</td>'

      write(lweb,*)'    <td WIDTH="65" VALIGN="MIDDLE" height="25" '
      write(lweb,*)'align="center" '
      write(lweb,*)'style="border-left: 1px solid rgb(0,0,0); '
      write(lweb,*)'border-right: 1px solid rgb(0,0,0); '
      write(lweb,*)'border-top: 3px solid rgb(0,128,0); '
      write(lweb,*)'border-bottom: 3px solid rgb(0,128,0)">'
      write(lweb,*)'<p ALIGN="CENTER">ch</td>'

      write(lweb,*)'    <td WIDTH="66" VALIGN="middle" height="25" '
      write(lweb,*)'align="center" '
      write(lweb,*)'style="border-left: 1px solid rgb(0,0,0); '
      write(lweb,*)'border-right: 1px solid rgb(0,0,0); ' 
      write(lweb,*)'border-top: 3px solid rgb(0,128,0); '
      write(lweb,*)'border-bottom: 3px solid rgb(0,128,0)">'
      write(lweb,*)'<p ALIGN="CENTER">N<sub>1</sub></td>'

      write(lweb,*)'    <td WIDTH="66" VALIGN="middle" height="25" '
      write(lweb,*)'align="center" '
      write(lweb,*)'style="border-left: 1px solid rgb(0,0,0); '
      write(lweb,*)'border-right: 1px solid rgb(0,0,0); ' 
      write(lweb,*)'border-top: 3px solid rgb(0,128,0); '
      write(lweb,*)'border-bottom: 3px solid rgb(0,128,0)">'
      write(lweb,*)'<p ALIGN="CENTER">N<sub>2</sub></td>'

      write(lweb,*)'    <td WIDTH="70" VALIGN="middle" height="25" '
      write(lweb,*)'align="center" '
      write(lweb,*)'style="border-left: 1px solid rgb(0,0,0); '
      write(lweb,*)'border-right: 1px solid rgb(0,0,0); ' 
      write(lweb,*)'border-top: 3px solid rgb(0,128,0); '
      write(lweb,*)'border-bottom: 3px solid rgb(0,128,0)">'
      write(lweb,*)
     &'<font FACE="Symbol"><p ALIGN="CENTER">D</font>M</td>'

      write(lweb,*)'    <td WIDTH="70" VALIGN="middle" height="25" '
      write(lweb,*)'align="center" '
      write(lweb,*)'style="border-left: 1px solid rgb(0,0,0); '
      write(lweb,*)'border-right: 1px solid rgb(0,0,0); ' 
      write(lweb,*)'border-top: 3px solid rgb(0,128,0); '
      write(lweb,*)'border-bottom: 3px solid rgb(0,128,0)">'
      write(lweb,*)'<font FACE="Symbol">D</font>M<sub>200</sub></td>'

      write(lweb,*)'    <td WIDTH="244" VALIGN="middle" height="25" '
      write(lweb,*)'align="center" '
      write(lweb,*)'style="border-left: 1px solid rgb(0,0,0); '
      write(lweb,*)'border-right: 1px solid rgb(0,0,0); ' 
      write(lweb,*)'border-top: 3px solid rgb(0,128,0); '
      write(lweb,*)'border-bottom: 3px solid rgb(0,128,0)" '
      write(lweb,*)'bordercolor="#000000">'
      write(lweb,*)'<p ALIGN="CENTER">Data (GeV/c<sup>2</sup>)</td>'

      write(lweb,*)'    <td WIDTH="50" VALIGN="middle" height="25" '
      write(lweb,*)'align="center" '
      write(lweb,*)'style="border-left: 1px solid rgb(0,0,0); '
      write(lweb,*)'border-right: 3px solid rgb(0,128,0); ' 
      write(lweb,*)'border-top: 3px solid rgb(0,128,0); '
      write(lweb,*)'border-bottom: 3px solid rgb(0,128,0)" '
      write(lweb,*)'bordercolor="#000000">'
      write(lweb,*)'<p ALIGN="CENTER">weight</td>'

      write(lweb,*)'  </tr>'
      write(lweb,*)'  <tr>'
c----------------------------------------------------------
      sumlum=0.
      do i=1,nchan
         n1tot(i)=0
         n2tot(i)=0
         n119x(i)=0
         n219x(i)=0
         n120x(i)=0
         n220x(i)=0
      enddo
      necmtot=necm+3
      do iecm=1,necmtot
      write(lweb,*)
     &'    <td WIDTH="50" VALIGN="middle" ROWSPAN="7" height="172" '
      write(lweb,*)'style="border-left: 3px solid rgb(0,128,0); '
      write(lweb,*)'border-right: 1px solid rgb(0,0,0); ' 
      write(lweb,*)'border-top: 1px solid rgb(0,0,0); '
      write(lweb,*)'border-bottom: 3px solid rgb(0,128,0)">'
      if(iecm.le.necm) then
        write(lweb,'(A,f7.3,A)') '<p ALIGN="CENTER">',ecm(iecm),'</td>'
      elseif(iecm.eq.necmtot) then
        write(lweb,*)'<p ALIGN="CENTER">All available energies</td>' 
      elseif(iecm.eq.necmtot-2) then
        write(lweb,*)'<p ALIGN="CENTER">19X energies</td>' 
      elseif(iecm.eq.necmtot-1) then
        write(lweb,*)'<p ALIGN="CENTER">20X energies</td>' 
      endif
      write(lweb,*)
     &'<td WIDTH="50" VALIGN="middle" ROWSPAN="7" height="172" '
      write(lweb,*)'style="border-left: 1px solid;' 
      write(lweb,*)'border-right: 1px solid; border-top: 1px solid; '
      write(lweb,*)'border-bottom: 3px solid rgb(0,128,0)">'
      if(iecm.le.necm) then
        write(lweb,'(A,f7.3,A)') '<p ALIGN="CENTER">',lumi(iecm),'</td>'
        sumlum=sumlum+lumi(iecm)
      elseif(iecm.eq.necmtot) then
        write(lweb,'(A,f7.3,A)') '<p ALIGN="CENTER">',sumlum,'</td>' 
      elseif(iecm.eq.necmtot-2) then
        do j=1,necm
          if(ecm(j).gt.190..and.ecm(j).lt.202.) then
            sumlum19x=sumlum19x+lumi(j)
          endif
        enddo
        write(lweb,'(A,f7.3,A)')'<p ALIGN="CENTER">',sumlum19x,'</td>'
      elseif(iecm.eq.necmtot-1) then
        do j=1,necm
          if(ecm(j).gt.202..and.ecm(j).lt.209.) then
            sumlum20x=sumlum20x+lumi(j)
          endif
        enddo
        write(lweb,'(A,f7.3,A)')'<p ALIGN="CENTER">',sumlum20x,'</td>'
      endif

      do ic=1,7
        weblink=' '
        if(ic.le.3.or.ic.eq.5) then
          i=ic
          if(ic.eq.5)i=4
          titr=chw(i)
          if(iecm.le.9) then
            if(i.eq.1) write(weblink,'(A,i3,A)')
     &           '<a href="http://alephwww.cern.ch/~leibeng/mw_e_'
     &           ,nint(ecm(iecm)),'.html">'
            if(i.eq.2) write(weblink,'(A,i3,A)')
     &           '<a href="http://alephwww.cern.ch/~leibeng/mw_mu_'
     &           ,nint(ecm(iecm)),'.html">'
            if(i.eq.3) write(weblink,'(A,i3,A)')
     &           '<a href="http://alephwww.cern.ch/~boumedi/mw_tau_'
     &           ,nint(ecm(iecm)),'.html">'
            if(i.eq.4) write(weblink,'(A,i3,A)')
     &           '<a href="http://alephwww.cern.ch/~moutouss/mw_4q_'
     &           ,nint(ecm(iecm)),'.html">'
          endif
          if(iecm.le.necm)then
            nn1=n1(i,iecm)
            nn2=n2(i,iecm)
            n1tot(i)=n1tot(i)+nn1
            n2tot(i)=n2tot(i)+nn2
            if(ecm(iecm).gt.190..and.ecm(iecm).lt.202.) then
              n119x(i)=n119x(i)+nn1
              n219x(i)=n219x(i)+nn2
            elseif(ecm(iecm).ge.202.and.ecm(iecm).lt.209.) then
              n120x(i)=n120x(i)+nn1
              n220x(i)=n220x(i)+nn2
            endif
            dmx=dmxpect(i,iecm)
            dmx2=dmx*sqrt(lumi(iecm)/200.)
            x=xmw(i,iecm)
            ex=estat(i,iecm)
            sx1=es1(i,iecm)
            sx2=es2(i,iecm)
            poids=wei_mtot(i,iecm)
          elseif(iecm.eq.necmtot) then
             nn1=n1tot(i)
             nn2=n2tot(i)
             dmx=exm_chan_e(i)
             dmx2=0.
             x=xm_chan(i)
             ex=exm_chan(i)
             sx1=sxm_chan1(i)
             sx2=sxm_chan2(i)
             poids=wei_mall(i)
           elseif(iecm.eq.necmtot-2) then
             nn1=n119x(i)
             nn2=n219x(i)
             dmx=exm_chan_e19x(i)
             dmx2=dmx*sqrt(sumlum19x/200.)
             x=xm_chan19x(i)
             ex=exm_chan19x(i)
             sx1=sxm_chan119x(i)
             sx2=sxm_chan219x(i)
             poids=wei_mall19x(i)
           elseif(iecm.eq.necmtot-1) then
             nn1=n120x(i)
             nn2=n220x(i)
             dmx=exm_chan_e20x(i)
             dmx2=dmx*sqrt(sumlum20x/200.)
             x=xm_chan20x(i)
             ex=exm_chan20x(i)
             sx1=sxm_chan120x(i)
             sx2=sxm_chan220x(i)
             poids=wei_mall20x(i)
          endif      
        elseif(ic.eq.4) then
          titr='<font face="Vivante-DTC">l</font>' 
          if(iecm.le.necm)then
             nn1=n1(1,iecm)+n1(2,iecm)+n1(3,iecm)
             nn2=n2(1,iecm)+n2(2,iecm)+n2(3,iecm)
             dmx=emlvqq_e(iecm)
             dmx2=dmx*sqrt(lumi(iecm)/200.)
             x=xmlvqq(iecm)
             ex=emlvqq(iecm)
             sx1=eslvqq(iecm)
             poids=0.
          elseif(iecm.eq.necmtot) then
             nn1=n1tot(1)+n1tot(2)+n1tot(3)
             nn2=n2tot(1)+n2tot(2)+n2tot(3)
             dmx=em_lvqq_e
             dmx2=0.
             x=xm_lvqq
             ex=em_lvqq
             sx1=es_lvqq1
             sx2=es_lvqq2
             poids=0.
          elseif(iecm.eq.necmtot-2) then
             nn1=n119x(1)+n119x(2)+n119x(3)
             nn2=n219x(1)+n219x(2)+n219x(3)
             dmx=em_lvqq_e19x
             dmx2=dmx*sqrt(sumlum19x/200.)
             x=xm_lvqq19x
             ex=em_lvqq19x
             sx1=es_lvqq119x
             sx2=es_lvqq219x
             poids=0.
          elseif(iecm.eq.necmtot-1) then
             nn1=n120x(1)+n120x(2)+n120x(3)
             nn2=n220x(1)+n220x(2)+n220x(3)
             dmx=em_lvqq_e20x
             dmx2=dmx*sqrt(sumlum20x/200.)
             x=xm_lvqq20x
             ex=em_lvqq20x
             sx1=es_lvqq120x
             sx2=es_lvqq220x
             poids=0.
          endif
        elseif(ic.eq.6) then
          titr='all' 
          if(iecm.le.necm)then
             nn1=n1(1,iecm)+n1(2,iecm)+n1(3,iecm)+n1(4,iecm)
             nn2=n2(1,iecm)+n2(2,iecm)+n2(3,iecm)+n2(4,iecm)
             dmx=emtot_e(iecm)
             dmx2=dmx*sqrt(lumi(iecm)/200.)
             x=xmtot(iecm)
             ex=emtot(iecm)
             sx1=estot1(iecm)
             sx2=estot2(iecm)
             poids=0.
          elseif(iecm.eq.necmtot) then
             nn1=n1tot(1)+n1tot(2)+n1tot(3)+n1tot(4)
             nn2=n2tot(1)+n2tot(2)+n2tot(3)+n2tot(4)
             dmx=0.
             dmx2=0.
             x=xm_all
             ex=em_all
             sx1=es_all1
             sx2=es_all2
             poids=0.
          elseif(iecm.eq.necmtot-2) then
             nn1=n119x(1)+n119x(2)+n119x(3)+n119x(4)
             nn2=n219x(1)+n219x(2)+n219x(3)+n219x(4)
             dmx=0.
             dmx2=0.
             x=xm_all19x
             ex=em_all19x
             sx1=es_all119x
             sx2=es_all219x
             poids=0.
          elseif(iecm.eq.necmtot-1) then
             nn1=n120x(1)+n120x(2)+n120x(3)+n120x(4)
             nn2=n220x(1)+n220x(2)+n220x(3)+n220x(4)
             dmx=0.
             dmx2=0.
             x=xm_all20x
             ex=em_all20x
             sx1=es_all120x
             sx2=es_all220x
             poids=0.
          endif
        elseif(ic.eq.7) then
          titr='M<SMALL>4q</SMALL>-M<SMALL>'
     &          //'<FONTface=Vivante-DTC>l</FONT></SMALL>' 
          if(iecm.le.necm)then
             nn1=n1(1,iecm)+n1(2,iecm)+n1(3,iecm)+n1(4,iecm)
             nn2=n2(1,iecm)+n2(2,iecm)+n2(3,iecm)+n2(4,iecm)
             dmx=emwdiff_e(iecm)
             dmx2=0.
             x=mwdiff(iecm)
             ex=emwdiff(iecm)
             sx1=0.
             sx2=0.
             poids=0.
          elseif(iecm.eq.necmtot) then
             nn1=n1tot(1)+n1tot(2)+n1tot(3)+n1tot(4)
             nn2=n2tot(1)+n2tot(2)+n2tot(3)+n2tot(4)
             dmx=0.
             dmx2=0.
             x=mwdiftot
             ex=emwdiftot
             sx1=0.
             sx2=0.
             poids=0.
          elseif(iecm.eq.necmtot-2) then
             nn1=0
             nn2=0
             dmx=0.
             dmx2=0.
             x=mwdif19x
             ex=emwdif19x
             sx1=0.
             sx2=0.
             poids=0.
          elseif(iecm.eq.necmtot-1) then
             nn1=0
             nn2=0
             dmx=0.
             dmx2=0.
             x=mwdif20x
             ex=emwdif20x
             sx1=0.
             sx2=0.
             poids=0.
          endif
        endif

        if(ic.gt.1) write(lweb,*)'<tr>'
        tit_align=' VALIGN="MIDDLE" height="19" align="center" '
        write(lweb,*)'<td WIDTH="65"'//tit_align
        write(lweb,*)'style="border: 1px solid;'
        if(ic.eq.7)write(lweb,*)
     &  ' border-bottom: rgb(0,128,0) 3px solid; '
        write(lweb,*)' ">',weblink
        write(lweb,*),'<p ALIGN="CENTER">',titr,'</a></p></td>'

        write(lweb,*)'<td WIDTH="66"'//tit_align
        write(lweb,*)'style="border: 1px solid;'
        if(ic.eq.7)write(lweb,*)
     &       'border-bottom: rgb(0,128,0) 3px solid; '
        if(nn1.le.0.) then
          write(lweb,'(A)')'">-</td>'
        else
          write(lweb,'(A,i4,A)')' ">',nn1,'</td>'
        endif

        write(lweb,*)'<td WIDTH="66"'//tit_align
        write(lweb,*)'style="border: 1px solid;'
        if(ic.eq.7)write(lweb,*)
     &       'border-bottom: rgb(0,128,0) 3px solid; '
        if(nn2.le.0.) then
          write(lweb,'(A)')'">-</td>'
        else
          write(lweb,'(A,i4,A)')' ">',nn2,'</td>'
        endif

        write(lweb,*)'<td WIDTH="70"'//tit_align
        write(lweb,*)'style="border: 1px solid;'
        if(ic.eq.7)write(lweb,*)
     &       'border-bottom: rgb(0,128,0) 3px solid; '
        if(dmx.le.0.005) then
          write(lweb,'(A)')'">-</td>'
        else
          write(lweb,'(A,f5.3,A)')'">',dmx,'</td>'
        endif
        write(lweb,*)'<td WIDTH="70"'//tit_align
        write(lweb,*)'style="border: 1px solid;'
        if(ic.eq.7)write(lweb,*)
     &       'border-bottom: rgb(0,128,0) 3px solid; '
        if(dmx2.le.0.005) then
          write(lweb,'(A)')'">-</td>'
        else
          write(lweb,'(A,f5.3,A)')'">',dmx2,'</td>'
        endif

        write(lweb,*)'<td WIDTH="244"'//tit_align
        write(lweb,*)'style="border: 1px solid; '
      if(ic.eq.7)
     &       write(lweb,*)'border-bottom: rgb(0,128,0) 3px solid; '
        if(x.eq.0.) then
          write(lweb,'(A)')'">-</td>'
        else
          write(lweb,'(A,f7.3,A,f5.3)')
     &         ' ">',x,'&plusmn;',ex
          if(ic.le.6) write(lweb,'(A,f5.3)')'&plusmn;',sx1
          if(ic.eq.5.or.ic.eq.6) write(lweb,'(A,f5.3)')'&plusmn;',sx2
        endif
      write(lweb,*)'</td>'
      write(lweb,*)'<td WIDTH="50"'//tit_align
      write(lweb,*)'style="border: 1px solid; '
      write(lweb,*)'border-right: 3px solid rgb(0,128,0); ' 
      if(ic.eq.7)
     &       write(lweb,*)'border-bottom: rgb(0,128,0) 3px solid; '
      if(poids.le.0.) then
        write(lweb,'(A)')'">-</td>'
      else
        write(lweb,'(A,f4.2,A)')'">',poids,'</td>'
      endif
      write(lweb,*)'  </tr>'
      enddo                     !loop on channels at given Ecm
      enddo                     !loop on ecm


      write(lweb,*)'</table>'
c------------------------------------------------------------------------------------
      write(lweb,*)'<hr>'
      write(lweb,*)
     &'<h2><big>Systematics</big><a NAME="Systematics"></a> </h2>'
      write(lweb,*)' '
      write(lweb,*)
     &'<p>LEP energy error is updated using tables below</p>'

      do iecm=1,necm
        if(nrange_sys(iecm).gt.0) then
          write(lweb,*)'<b><i>Valid for '
          do j=1,necm
            if(irange_sys(j,iecm).gt.0.)
     &           write(lweb,'(i4)')nint(ecmc(irange_sys(j,iecm)))
          enddo
          if(nrange_sys(iecm).eq.1) then
            write(lweb,*)' GeV sample</i></b>'
          elseif(nrange_sys(iecm).gt.1) then
            write(lweb,*)' GeV samples</i></b>'
          endif

      write(lweb,*)'<TABLE border=3' 
      write(lweb,*)'style="BORDER-BOTTOM: rgb(0,128,0) 3px solid; '
      write(lweb,*)'BORDER-LEFT: rgb(0,128,0) 3px solid; '
      write(lweb,*)'BORDER-RIGHT: rgb(0,128,0) 3px solid; '
      write(lweb,*)'BORDER-TOP: rgb(0,128,0) 3px solid" width=682>'
      write(lweb,*)'  <TBODY>'
      write(lweb,*)'  <TR>'
      write(lweb,*)'    <TD borderColor=#000000 width="40%">&nbsp;</TD>'
      write(lweb,*)
     &'    <TD borderColor=#000000 colSpan=5><FONT face=Symbol>'
      write(lweb,*)
     &'    <P align=center>D</FONT>M (MeV/c<SUP>2</SUP>)</P></TD></TR>'
      write(lweb,*)'  <TR>'
      write(lweb,*)'    <TD borderColor=#000000 width="40%">Source</TD>'
      write(lweb,*)
     &     '    <TD align=middle borderColor=#000000 width="12%">'
      write(lweb,*)'e<FONT face=Symbol>n</FONT>qq</TD>'
      write(lweb,*)
     &'    <TD align=middle borderColor=#000000 width="12%">'
      write(lweb,*)'<FONT face=Symbol>mn</FONT>qq</TD>'
      write(lweb,*)
     &'    <TD align=middle borderColor=#000000 width="12%">'
      write(lweb,*)'<FONT face=Symbol>tn</FONT>qq</TD>'
      write(lweb,*)
     &'    <TD align=middle borderColor=#000000 width="12%">'
      write(lweb,*)'<FONT face=Vivante-DTC>l</FONT>'
      write(lweb,*)' <FONT face=Symbol>n</FONT>qq</TD>'
      write(lweb,*)
     &'    <TD align=middle borderColor=#000000 width="12%">4q</TD>'
      write(lweb,*)'  </TR>'
      write(lweb,*)'  <TR>'
      write(lweb,*)'    <TD borderColor=#000000' 
      write(lweb,*)'    style="BORDER-BOTTOM: rgb(0,128,0) 2px solid;'
      write(lweb,*)' BORDER-TOP: rgb(0,128,0) 2px solid"' 
      write(lweb,*)'    width="40%">(a) Correlated errors</TD>'
      write(lweb,*)'    <TD borderColor=#000000 colSpan=5' 
      write(lweb,*)'    style="BORDER-BOTTOM: rgb(0,128,0) 2px solid;'
      write(lweb,*)' BORDER-TOP: rgb(0,128,0) 2px solid">&nbsp;'
      write(lweb,*)'</TD></TR>'

      do i=1,5
         est(i)=0.
      enddo
      do j=1,10
         write(lweb,*)
     &'   <TR><TD borderColor=#000000 width="40%">',chsys(j),'</TD>'
         do i=1,5
            if(i.le.3.or.i.eq.5) then
               ii=i
               if(i.eq.5) ii=4
               if(j.eq.10) 
     &              esys(j,ii,iecm)=decmc(iecm)*2./ecmc(iecm)*80.35
                write(lweb,'(A,i3,A)')
     &        '  <TD align=middle borderColor=#000000 width="12%">'
     &        ,nint(esys(j,ii,iecm)*1000.),'</TD>'
                est(i)=est(i)+(esys(j,ii,iecm)*1000.)**2
c                print *,chsys(j),esys(j,ii),int(esys(j,ii))
             else
                esl=es_lq(j,iecm)
                est(i)=est(i)+(esl*1000.)**2
                write(lweb,'(A,i3,A)')
     &'    <TD align=middle borderColor=#000000 width="12%">'
     &               ,nint(esl*1000.),'</TD>'
             endif
          enddo
          write(lweb,*)'  </TR>'
       enddo
      write(lweb,*)'  <TR>'
      write(lweb,*)'<TD borderColor=#000000 '
      write(lweb,*)'    style="BORDER-BOTTOM: rgb(0,128,0) 2px solid;'
      write(lweb,*)' BORDER-TOP: rgb(0,128,0) 2px solid" width="40%">'
      write(lweb,*)'(b) Uncorrelated errors</TD>'
      write(lweb,*)'<TD borderColor=#000000 colSpan=5' 
      write(lweb,*)'    style="BORDER-BOTTOM: rgb(0,128,0) 2px solid;'
      write(lweb,*)' BORDER-TOP: rgb(0,128,0) 2px solid">&nbsp;'
      write(lweb,*)'</TD></TR>'
      write(lweb,*)'  <TR>'
      write(lweb,*)'<TD borderColor=#000000 width="40%">'

      do j=11,14
         write(lweb,*)
     &'   <TR><TD borderColor=#000000 width="40%">',chsys(j),'</TD>'
         do i=1,5
            if(i.le.3.or.i.eq.5) then
               ii=i
               if(i.eq.5) ii=4
                write(lweb,'(A,i3,A)')
     &        '  <TD align=middle borderColor=#000000 width="12%">'
     &        ,int(esys(j,ii,iecm)*1000.),'</TD>'
                est(i)=est(i)+(esys(j,ii,iecm)*1000.)**2
             else
                esl=es_lq(j,iecm)
                est(i)=est(i)+(esl*1000.)**2
                write(lweb,'(A,i3,A)')
     &'    <TD align=middle borderColor=#000000 width="12%">'
     &               ,int(esl*1000.),'</TD>'
             endif
          enddo
          write(lweb,*)'  </TR>'
       enddo

      write(lweb,*)'  <TR>'
      write(lweb,*)'<TD borderColor=#000000 '
      write(lweb,*)'    style="BORDER-BOTTOM: rgb(0,0,0) 1px solid;'
      write(lweb,*)' BORDER-LEFT: rgb(0,0,0) 1px solid;'
      write(lweb,*)' BORDER-RIGHT: rgb(0,0,0) 1px solid;'
      write(lweb,*)' BORDER-TOP: rgb(0,128,0) 2px solid"' 
      write(lweb,*)'    width="40%">Total (a+b)</TD>'
      do i=1,5
      write(lweb,*)'<TD align=middle borderColor=#000000' 
      write(lweb,*)'    style="BORDER-BOTTOM: rgb(0,0,0) 1px solid;'
      write(lweb,*)' BORDER-LEFT: rgb(0,0,0) 1px solid;'
      write(lweb,*)' BORDER-RIGHT: rgb(0,0,0) 1px solid;'
      write(lweb,*)' BORDER-TOP: rgb(0,128,0) 2px solid"' 
      write(lweb,'(A,i3,A)')'    width="12%">',int(sqrt(est(i))),'</TD>'
      enddo
      write(lweb,*)'</TR></TBODY></TABLE>'
      write(lweb,*)'<P align=left>&nbsp;</P>'
      
      endif                     !check validity of table for this energy
      enddo                     !loop on ecm energies
      
c--------------------------------------------------------------------------------------
      write(lweb,*)
     &'<H2 align=left>Lep Energy Error (in MeV)</H2>'
      write(lweb,*)'<TABLE border=2 borderColorLight=#000000' 
      write(lweb,*)'style="BORDER-BOTTOM: rgb(0,128,0) solid;'
      write(lweb,*)'BORDER-LEFT: rgb(0,128,0) solid;'
      write(lweb,*)'BORDER-RIGHT: rgb(0,128,0) solid;'
      write(lweb,*)'BORDER-TOP: rgb(0,128,0) solid; TEXT-ALIGN: left"' 
      write(lweb,*)'width=682>'
      write(lweb,*)'  <TBODY>'
      write(lweb,*)'  <TR>'
      do i=1,necm
         write(lweb,'(A,I3,A)')
     &       ' <TD align=middle width="11%">',nint(ecmc(i)),'</TD>'
      enddo
      write(lweb,*)'</TR>'
      write(lweb,*)'  <TR>'
      do i=1,necm
         write(lweb,'(A,I3,A)')
     &       ' <TD align=middle width="11%">'
     &       ,nint(decmc(i)*1000.),'</TD>'
      enddo
      write(lweb,*)'</TR>'
      write(lweb,*)'</TBODY></TABLE>'
c      write(lweb,*)'<P align=left>&nbsp;</P>'

c--------------------------------------------------------------------------------------
      write(lweb,*)
     &'<H2 align=left>Lep Energy Error Correlation Matrix</H2>'
      write(lweb,*)'<TABLE border=2 borderColorLight=#000000' 
      write(lweb,*)'style="BORDER-BOTTOM: rgb(0,128,0) solid;'
      write(lweb,*)'BORDER-LEFT: rgb(0,128,0) solid;'
      write(lweb,*)'BORDER-RIGHT: rgb(0,128,0) solid;'
      write(lweb,*)'BORDER-TOP: rgb(0,128,0) solid; TEXT-ALIGN: left"' 
      write(lweb,*)'width=682>'
      write(lweb,*)'  <TBODY>'
      write(lweb,*)'  <TR>'
      write(lweb,*)'    <TD width="10%">&nbsp;</TD>'
      do i=1,necm
         write(lweb,'(A,I3,A)')
     &       ' <TD align=middle width="10%">',nint(ecmc(i)),'</TD>'
      enddo
      write(lweb,*)'</TR>'
      write(lweb,*)'  <TR>'
      do i=1,necm
         write(lweb,'(A,I3,A)')
     &       ' <TD align=middle width="10%">',nint(ecmc(i)),'</TD>'
         do j=1,necm
            if(abs(rhoecm(i,j)-1.0).lt.0.001) then
               write(lweb,'(A,i2,A)')
     &        ' <TD align=middle width="10%">',nint(rhoecm(i,j)),'</TD>'
            else
               write(lweb,'(A,f4.2,A)')
     &              ' <TD align=middle width="10%">',rhoecm(i,j),'</TD>'
            endif
         enddo
         write(lweb,*)'</TR>'
      enddo
      write(lweb,*)'</TBODY></TABLE>'
      write(lweb,*)'<P align=left>&nbsp;</P>'
      write(lweb,*)'<P align=left>&nbsp;</P>'

      write(lweb,*)' '
      write(lweb,*)'<p ALIGN="CENTER">'
      write(lweb,*)'<a HREF="'
     &//'http://alephwww.cern.ch/LEP2ANALYSIS/WTF/Analyses/mass.html">'
      write(lweb,*)'<img SRC="'
      write(lweb,*)
     &'http://alephwww.cern.ch/LEP2ANALYSIS/WTF/Picts/color_up.gif"'
      write(lweb,*)'BORDER="0" WIDTH="40" HEIGHT="41"></a></p>'
      write(lweb,*)'</body>'
      write(lweb,*)'</html>'


      end
