      program SN
* extraction of S/B from *.sn files
      implicit none
      integer i,iunit
      integer nfiles
      parameter (nfiles=11)
      real xsec(nfiles)
      real effsel(nfiles),deltaeffsel(nfiles)
      real effpresel(nfiles),deltaeffpresel(nfiles)
      real xsecpre(nfiles),deltaxsecpre(nfiles)
      real xsecpost(nfiles),deltaxsecpost(nfiles)
      integer nini(nfiles)
      integer npre(nfiles),deltanpre(nfiles)
      integer npost(nfiles),deltanpost(nfiles)
      character*16 channel(nfiles)
      real intlumi
      parameter(intlumi=30000) ! 3 years integrated luminosity (pb^-1)
      real s,b,sb
      real delta_s,delta_b,delta_sb
      real deltavtb
c---  
      open(unit=21,file='ewtop_s.sn',status='OLD')
      open(unit=22,file='ewtop_t.sn',status='OLD')
      open(unit=23,file='wt.sn',status='OLD')
      open(unit=24,file='wt_2l.sn',status='OLD')
      open(unit=25,file='ttbar.sn',status='OLD')
      open(unit=26,file='tt_2l.sn',status='OLD')
      open(unit=27,file='wbb.sn',status='OLD')
      open(unit=28,file='wcc.sn',status='OLD')
      open(unit=29,file='ww_incl.sn',status='OLD')
      open(unit=30,file='zw_incl.sn',status='OLD')
      open(unit=31,file='zz_incl.sn',status='OLD')
c      open(unit=32,file='zw_bbln.sn',status='OLD')

* Cross sections times BR(e,mu), picobarns:
c$$$      xsec(1)=2.22
c$$$      xsec(2)=54.2
c$$$      xsec(3)=17.8
c$$$      xsec(4)=246
c$$$      xsec(5)=66.7
      xsec(1)=2.38
      xsec(2)=54.4
      xsec(3)=16.6
      xsec(4)=2.77
      xsec(5)=244.5
      xsec(6)=40.8
      xsec(7)=77.8
      xsec(8)=311.2
* Wbb: 300 pb, Wcc: 1200 pb according to Sergey Slabospitsky
      xsec(9)=74.2
      xsec(10)=28.7
* from: http://cmsdoc.cern.ch/~spitskyc/sigtab/node2.html
      xsec(11)=28.7*28.7/74.2 ! in realta' non so quanto valga!
c      xsec(12)=28.7*0.222*0.21

* Efficiencies:
      do i=1,nfiles
         iunit=20+i
         read(iunit,*) effpresel(i)
         read(iunit,*) deltaeffpresel(i)
         read(iunit,*) effsel(i)
         read(iunit,*) deltaeffsel(i)
         xsecpre(i)=xsec(i)*effpresel(i)
         xsecpost(i)=xsec(i)*effsel(i)
         deltaxsecpre(i)=xsec(i)*deltaeffpresel(i)
         deltaxsecpost(i)=xsec(i)*deltaeffsel(i)
         nini(i)=xsec(i)*intlumi
         npre(i)=xsecpre(i)*intlumi
         npost(i)=xsecpost(i)*intlumi
         deltanpre(i)=deltaxsecpre(i)*intlumi
         deltanpost(i)=deltaxsecpost(i)*intlumi
      enddo

      channel(1) = 'W* (s-channel):' 
      channel(2) = 'Wg (t-channel):' 
      channel(3) = 'Wt->1l:'
      channel(4) = 'Wt->2l:'
      channel(5) = 'tt->1l:'
      channel(6) = 'tt->2l:'
      channel(7) = 'Wbb:'
      channel(8) = 'Wcc:'
      channel(9) = 'WW (incl):'
      channel(10) = 'ZW (incl):'
      channel(11) = 'ZZ (incl):'
c      channel(12) = 'ZW->bbln:'

      write(6,*) '--------------------------'
      write(6,*) 'xsec*BR(e,mu):' 
      write(6,*) '--------------------------'
      do i=1,nfiles
         write(6,*) channel(i),xsec(i)," pb (",nini(i)," evts)"
      enddo
      s=xsec(1)+xsec(2)+xsec(3)+xsec(4)
      b=xsec(5)+xsec(6)+xsec(7)+xsec(8)+xsec(9)+xsec(10)+xsec(11)
c     &     +xsec(12)
      write(6,*) 'S(single top)/B = ',s/b
      s=xsec(1)
      b=xsec(2)+xsec(3)+xsec(4)+xsec(5)+xsec(6)+xsec(7)+xsec(8)+xsec(9)
     &     +xsec(10)+xsec(11)
c     &     +xsec(12)
      write(6,*) 'S(s-channel)/B = ',s/b
      s=xsec(2)
      b=xsec(1)+xsec(3)+xsec(4)+xsec(5)+xsec(6)+xsec(7)+xsec(8)+xsec(9)
     &     +xsec(10)+xsec(11)
c     &     +xsec(12)
      write(6,*) 'S(t-channel)/B = ',s/b
      write(6,*) '--------------------------'
      write(6,*) 'After preselection:'
      write(6,*) '--------------------------'
      do i=1,nfiles
         write(6,*) channel(i),xsecpre(i)," +-",deltaxsecpre(i)," pb ("
     &        ,npre(i)," +-",deltanpre(i)," evts)"
      enddo
      s=xsecpre(1)+xsecpre(2)+xsecpre(3)+xsecpre(4)
      b=xsecpre(5)+xsecpre(6)+xsecpre(7)+xsecpre(8)+xsecpre(9)
     &     +xsecpre(10)+xsecpre(11)
c     &     +xsecpre(12)
      write(6,*) 'S(single top)/B = ',s/b
      s=xsecpre(1)
      b=xsecpre(2)+xsecpre(3)+xsecpre(4)+xsecpre(5)+xsecpre(6)
     &     +xsecpre(7)+xsecpre(8)+xsecpre(9)+xsecpre(10)+xsecpre(11)
c     &     +xsecpre(12)
      write(6,*) 'S(s-channel)/B = ',s/b
      s=xsecpre(2)
      b=xsecpre(1)+xsecpre(3)+xsecpre(4)+xsecpre(5)+xsecpre(6)
     &     +xsecpre(7)+xsecpre(8)+xsecpre(9)+xsecpre(10)+xsecpre(11)
c     &     +xsecpre(12)
      write(6,*) 'S(t-channel)/B = ',s/b
      call latex(effpresel,deltaeffpresel,xsecpre,deltaxsecpre,npre,
     &     deltanpre,1)
      write(6,*) '--------------------------'
      write(6,*) 'After selection:'
      write(6,*) '--------------------------'
      do i=1,nfiles
         write(6,*) channel(i),xsecpost(i)," +-",deltaxsecpost(i),
     &        " pb (",npost(i)," +-",deltanpost(i)," evts)"
      enddo
      s=xsecpost(1)+xsecpost(2)+xsecpost(3)+xsecpost(4)
      delta_s=sqrt(deltaxsecpost(1)**2+deltaxsecpost(2)**2+
     &     deltaxsecpost(3)**2+deltaxsecpost(4)**2)
      b=xsecpost(5)+xsecpost(6)+xsecpost(7)+xsecpost(8)+xsecpost(9)
     &     +xsecpost(10)+xsecpost(11)
c     &     +xsecpost(12)
      delta_b=sqrt(deltaxsecpost(5)**2+deltaxsecpost(6)**2
     &     +deltaxsecpost(7)**2+deltaxsecpost(8)**2+deltaxsecpost(9)**2
     &     +deltaxsecpost(10)**2+deltaxsecpost(11)**2)
c     &     +deltaxsecpost(12)**2)
      sb=s/b
      delta_sb=sb*sqrt((delta_s/s)**2+(delta_b/b)**2)
      write(6,*) 'S(single top)/B = ',sb," +-",delta_sb
      deltavtb=0.5*sqrt((s+b)*intlumi)/(s*intlumi)
      write(6,*) 'Delta Vtb / Vtb (3 yr)= ',deltavtb*100,' %'
      s=xsecpost(1)
      delta_s=deltaxsecpost(1)
      b=xsecpost(2)+xsecpost(3)+xsecpost(4)+xsecpost(5)+xsecpost(6)
     &     +xsecpost(7)+xsecpost(8)+xsecpost(9)+xsecpost(10)
     &     +xsecpost(11)
c     &     +xsecpost(12)
      delta_b=sqrt(deltaxsecpost(2)**2+deltaxsecpost(3)**2+
     &     deltaxsecpost(4)**2+deltaxsecpost(5)**2+deltaxsecpost(6)**2
     &     +deltaxsecpost(7)**2+deltaxsecpost(8)**2+deltaxsecpost(9)**2
     &     +deltaxsecpost(10)**2+deltaxsecpost(11)**2)
c     &     +deltaxsecpost(12)**2)
      sb=s/b
      delta_sb=sb*sqrt((delta_s/s)**2+(delta_b/b)**2)
      write(6,*) 'S(s-channel)/B = ',sb," +-",delta_sb
      deltavtb=0.5*sqrt((s+b)*intlumi)/(s*intlumi)
      write(6,*) 's/t = ',xsecpost(1)/xsecpost(2)
      write(6,*) 'Delta Vtb / Vtb (3 yr)= ',deltavtb*100,' %'
      write(6,*) '(Eff. s-channel:',effsel(1)*100.,'%)'
      s=xsecpost(2)
      delta_s=deltaxsecpost(2)
      b=xsecpost(1)+xsecpost(3)+xsecpost(4)+xsecpost(5)+xsecpost(6)
     &     +xsecpost(7)+xsecpost(8)+xsecpost(9)+xsecpost(10)
     &     +xsecpost(11)
c     &     +xsecpost(12)
      delta_b=sqrt(deltaxsecpost(1)**2+deltaxsecpost(3)**2+
     &     deltaxsecpost(4)**2+deltaxsecpost(5)**2+deltaxsecpost(6)**2
     &     +deltaxsecpost(7)**2+deltaxsecpost(8)**2+deltaxsecpost(9)**2
     &     +deltaxsecpost(10)**2+deltaxsecpost(11)**2)
c     &     +deltaxsecpost(12)**2)
      sb=s/b
      delta_sb=sb*sqrt((delta_s/s)**2+(delta_b/b)**2)
      write(6,*) 'S(t-channel)/B = ',sb," +-",delta_sb
      write(6,*) 't/s = ',xsecpost(2)/xsecpost(1)
      deltavtb=0.5*sqrt((s+b)*intlumi)/(s*intlumi)
      write(6,*) 'Delta Vtb / Vtb (3 yr)= ',deltavtb*100,' %'
      write(6,*) '(Eff. t-channel:',effsel(2)*100.,'%)'
      call latex(effsel,deltaeffsel,xsecpost,deltaxsecpost,npost,
     &     deltanpost,2)
      
c      open(unit=20, file='output.opt', status='UNKNOWN')
c      write(20,*) effsel(1)*100.,sb,deltavtb*100.

      end
c----------------------------------------------------------

      subroutine latex(eff,deltaeff,x,deltax,n,deltan,iflag)
      implicit none
      integer iflag
      integer nfiles
      parameter (nfiles=11)
      integer i
      real x(nfiles),deltax(nfiles)
      real eff(nfiles),deltaeff(nfiles)
      integer n(nfiles),deltan(nfiles)
      character*32 text(nfiles)
      character*3 label

c--
      text(1) = '$W^* \\rightarrow 1l$' 
      text(2) = '$Wg \\rightarrow 1l$' 
      text(3) = '$Wt \\rightarrow 1l$'
      text(4) = '$Wt \\rightarrow 2l$'
      text(5) = '$t\\bar t \\rightarrow 1l$'
      text(6) = '$t\\bar t \\rightarrow 2l$'
      text(7) = '$Wb\\bar b \\rightarrow 1l$'
      text(8) = '$Wc\\bar c \\rightarrow 1l$'
      text(9) = '$WW (incl.)$'
      text(10) = '$WZ (incl.)$'
      text(11) = '$ZZ (incl.)$'
c      text(12) = '$WZ \\rightarrow l\\nu b\\bar b$'

      if (iflag.eq.1) label='pre'
      if (iflag.eq.2) label='sel'
      open(unit=20, file='table_'//label//'.tex', status='UNKNOWN')
      write(20,*) '\\documentclass[12pt,a4paper]{report}'
      write(20,*) '\\usepackage{rotating}'
      write(20,*) ''
      write(20,*) ''
      write(20,*) ''
      write(20,*) '\\begin{document}'
      write(20,*) '\\centering'
      write(20,*) '\\begin{table}[htb]'
      write(20,*) '\\centering'
      write(20,*) '{\\Large'
c      write(20,*) '\\rotatebox{90}{'
      write(20,*) '\\begin{tabular}{|l|c|c|c|}'
      write(20,*) '\\hline'
      write(20,*) ' & Selected/Generated & After selection'//
     *     '& Evts @ 30 fb$^{-1}$\\\\\\hline'
      do i=1,nfiles
         write(20,5) text(i),' & ',eff(i)*100,' $\\pm$',deltaeff(i)*100,
     &        ' ($\\%$) & ',x(i),' $\\pm$',deltax(i),' (pb) & ',n(i),
     &        ' $\\pm$',deltan(i),'\\\\'
         write(20,*) '\\hline'
      enddo
      write(20,*) '\\end{tabular}'
c      write(20,*) '}'
      write(20,*) '}'
      write(20,*) '\\end{table}'
      write(20,*) ''
      write(20,*) ''
      write(20,*) '\\end{document}'

 5    format(A,A,F12.3,A,F12.3,A,F12.4,A,F12.4,A,I6,A,I4,A)
      return
      end
