      subroutine hggtotal(inputfile,outputfile)
      implicit none 
      character*50 inputfile, outputfile
#include "numcomp.inc"

      double precision mhiggs, tau, Etot, tauT, mtop
      double precision x1, x2, mur, muf, mbot, mtop_pass
      double precision SSC, botscheme, mhiggs_pass
      character*60 pdfset, porder, collider
      double precision asZ, alphas
      integer isetval, ival, iii, imhiggs, im4up, imu
      double precision qcd_lo_v, qcd_lo_p, qcd_lo_m
      double precision qcd_nlo_v, qcd_nlo_p, qcd_nlo_m
      double precision qcd_nnlo_v, qcd_nnlo_p, qcd_nnlo_m
      double precision qcd_pure_v, qcd_pure_p, qcd_pure_m
      double precision sigma_all_v, sigma_all_p, sigma_all_m

      double precision upv,dnv,usea,dsea,str,chm,bot,glu
      double precision Lgg,Lqg,Lqqb,Lqq,Lq1q2
      double precision deltaLO, deltaNLO, deltaNNLO1, deltaNNLO2
      double precision  lumigg, lumiqq, lumiqg, lumiqqb
      double precision lumiq1q2, ris
      double precision z2, z3, NF, L, Lt, Lr
      double complex abornt, abornb
      double precision sborntt, sbornbb, sborntb
      double precision deltaRadja 
      double precision total, etotal, sbornttSM

      double precision pref, rdelta(3), edelta(3)
      double complex bornampl, borntop, bornall, borninter
      double precision conv, tn, ftn,Gf, pi, cnorm
      double precision  mW, mZ
      double precision  GammaW, GammaZ

      double precision rplus, eplus, rhard, ehard
      double precision fitted

      double precision xgg
C------------ Wilson coefficient definitions
      double precision  C1q, C2q, lambda, C1w, C2w, deltaEWK
      double precision  weak1(idm), weak2(idm)
      double precision  weak3(idm)
      double precision  qcd1(idm), qcd2(idm), qcd3(idm)
      double precision  gen4qcd(idm)
      double precision  eweak1(idm), eweak2(idm)
      double precision  eqcd2(idm), eqcd3(idm),eqcd1(idm)
      double precision weak_lo(idm), weak_nlo(idm)
      double precision weak_nnlo_m5(idm), weak_nnlo_p5(idm)
      double precision qcd_lo(idm), qcd_nlo(idm)
      double precision sigma_all(idm),qcd_pure(idm)
      double precision weak_nnlo_known(idm),qcd_bot(idm)
      double precision qcd_nnlo(idm)
      double precision m4up, m4down, L4down, L4up, Lm1, Lm2
      double precision C1q4up, C1q4down, C2q4up, C2q4down
      double precision C2q4uptop, C2q4downtop, C2q4updown
      double precision C2qdiff


      double precision muval(4), Ktb, Kbb

      integer pset, jj, mm, ll, llmu, lliste
      integer Ngrids, igrids
 

      common /pset/ pset
      COMMON /lumiB/ collider, pdfset, muf
      COMMON /kinem/ tau
      common /xsection/ L, Lt, Lr, alphas, NF 
      common /consts/ pi, z2, z3       
      common /Ngrids/ Ngrids
      
      external INTE2, INTE1, INTE0

cccc======================================
      pi=3.1415926535897D0
      z2=pi**2/6.0D0
      z3=1.2020569031595942854D0
      NF=5d0

      call unlink('results.dat')
      open(20,FILE='results.dat',STATUS='NEW')
      write(20, *), "Mhiggs ", " mu ", " m4up  ", " m4down ", " sigma ",
     &              " Dpdf+ ", " Dpdf- "
      close(20)

       collider = 'LHC'
      collider = 'TEVATRON'
      if (collider.eq.'LHC') then
         Etot = 10000d0         
C         Etot = 14000d0
      elseif (collider.eq.'TEVATRON') then 
         Etot = 1960d0         
      endif
      
      do imhiggs =0,26
      mhiggs= 110d0+5d0*imhiggs 
      do imu=0,2
      muf= mhiggs/(2d0**imu)
      mur = muf 

C---- at the moment this should not be touched, in order  to be 
C---- consistent  with the electroweak corrections from  
      
      mtop = 170.9d0 

C------ 4th generation ----------
      do im4up=0,1
      m4down = 300d0+100d0*im4up 
      m4up = m4down+50d0+10d0*dlog(mhiggs/115d0)

C------
      botscheme = 1d0
      if  (botscheme.eq.0d0) then
C     ---- pole  scheme for bottom quark
         mbot = 4.75d0
      else
C     -----  running mass  scheme -----        
C     bottom mass at 10GeV
         mbot = 3.609d0 
         mbot = mbot*dexp(-2d0*0.179d0/pi*dlog(mur/10d0))
      endif
      
      L=dlog(muf**2/mhiggs**2)
      Lr=dlog(mur**2/mhiggs**2)
      Lt=dlog(mhiggs**2/mtop**2)
      L4up=dlog(mhiggs**2/m4up**2)
      L4down=dlog(mhiggs**2/m4down**2)

      tau = mhiggs**2/Etot**2
      

      mW      = 80.398d0
      mZ      = 91.1876d0
      
      GammaW  = 2.141d0
      GammaZ  = 2.4952d0

C----- SM Wilson coefficient
            
      Lm1=Lt
      C1q=11d0/4d0
      C2q = 2777d0/288d0+19d0/16d0*(Lr+Lm1)
     #+5d0*(-67d0/96d0 +1d0/3d0*(Lr+Lm1))


      C1q4up = C1q
      C1q4down = C1q 

      Lm1 = L4up
      C2q4up = 2777d0/288d0+19d0/16d0*(Lr+Lm1)
     #+5d0*(-67d0/96d0 +1d0/3d0*(Lr+Lm1))

      Lm1 = L4down
      C2q4down = 2777d0/288d0+19d0/16d0*(Lr+Lm1)
     #+5d0*(-67d0/96d0 +1d0/3d0*(Lr+Lm1))



      C2q4uptop = -77d0/288d0
      C2q4downtop =  -77d0/288d0
      C2q4updown =  -77d0/288d0
  
      C2qdiff = C2q4up+C2q4down+C2q4uptop
     &+C2q4downtop+C2q4updown
      C2qdiff = C2qdiff - 2d0*C2q 
      C2qdiff = 2d0/3d0*C2qdiff




C---- reads EWK  corrections from Fig. 21 of http://arXiv.org/pdf/0809.3667
C---- by Actis et al
      mtop_pass = 172.7d0 
      mhiggs_pass = mhiggs*mtop/mtop_pass
      if ((mhiggs_pass.gt.12000d0).and.(mhiggs_pass.lt.49800d0)) then 
         deltaEWK=fitted(mhiggs_pass)/100d0
         print*, "P=", deltaEWK 
      else
C----  Radja's light fermion two-loop EWK extended to the complex mass scheme
C----  (analytic formula in  Bonciani et al)
         deltaEWK=deltaRadja(mW,mZ,GammaW,
&                   GammaZ,mtop,mhiggs)
         print*, "R=", deltaEWK
      endif

      lambda=deltaEWK/2d0
C---  Wilson coefficient computed in arXiv:0811.3458 [hep-ph]  
C---  by Anastasiou et al
      C1w=7d0/6d0*1d0
C------  this is not computed but only estimated for the four-loop 
C------  mix QCD-EWK diagrams.  The  result is rather insensitive to 
C------   "reasonable" values (close to similar top-quark Wilson coeffs.) 
         C2w=-10d0
C--------------------------

C---- Top-bottom and bottom-bottom K-factors!      
C--- This is done properly only for TEVATRON and MRST2004 pdfs      
C--- It will be replaced with an exact calculation  soon 

      if (botscheme.eq.0d0)  then 
C---- pole scheme K-factors ----------
       Ktb= 1.15d0+.002503998382d0*mhiggs;
       Kbb= 1.43d0+.001701722336d0*mhiggs;
      else
C---- running mass K-factors ---------
         Kbb= 1.97d0+.004240955309d0*mhiggs;
         Ktb= 1.32d0+.003902054591d0*mhiggs;
      endif
C----  sloppy job here; it will be  replaced with some proper code!
C----------- FOR MSTW pdfs, we have smaller K-factors
      if (pdfset.eq.'MSTW2008') then
         Kbb = Kbb*0.87d0
         Ktb = Ktb*0.87d0
      endif
C--- Conservative very good estimate with top-bottom and bottom-bottom 
C--- only at LO: 
C      Kbb =1d0
C      Ktb =1d0 
C---------------------------------------------------------------------
         do pset = 0,2 


         if (pset.eq.0) then 
            porder = 'LO' 
C     pdfset = 'MRSTLO'
            pdfset = 'MSTW2008LO'
         elseif (pset.eq.1) then 
            porder = 'NLO' 
C     pdfset = 'MRST2001'
C     pdfset = 'MRST2004NLO'
            pdfset = 'MSTW2008NLO'
         elseif (pset.eq.2) then 
            porder = 'NNLO' 
C     pdfset = 'MRSTNNLO'
C     pdfset = 'MRST2004NNLO'
C     pdfset = 'MRST2006NNLO'
            pdfset = 'MSTW2008NNLO'
         endif



         
C---- number of extra grids (not including  the mean)
         Ngrids = 40
         
                 
         if (pdfset .EQ. 'MRSTLO') then       
            asZ = 0.130D0        
         elseif (pdfset .EQ. 'MRST2001') then 
            asZ  =  0.119D0
         elseif (pdfset .EQ. 'MRSTNNLO') then 
            asZ  =  0.1155D0
         elseif (pdfset .EQ. 'MRST2004NLO') then 
            asZ  =  0.1205D0
         elseif (pdfset .EQ. 'MRST2004NNLO') then 
            asZ  =  0.1167D0
         elseif (pdfset .EQ. 'MRST2006NNLO') then 
            asZ  =  0.1191D0
C--------latest pdf-sets -------------------------------
         elseif (pdfset .EQ. 'MSTW2008LO') then 
            asZ  =  0.13939D0
         elseif (pdfset .EQ. 'MSTW2008NLO') then 
            asZ  =  0.12018D0
         elseif (pdfset .EQ. 'MSTW2008NNLO') then 
            asZ  =  0.11707D0
C--------------------------------------------------------
         endif
         alphas = SSC(mur,asZ,porder)         
C-----computes born amplitudes with exact mquark
         abornt=bornampl(mtop,mhiggs)
         sbornttSM= cdabs(abornt)**2d0
C----- add the fourth generation to the top amplitude....
         abornt=bornampl(m4up,mhiggs)+abornt
         abornt=bornampl(m4down,mhiggs)+abornt
C--------------------------------------------------
         abornb=bornampl(mbot,mhiggs)*1d0

C--------------------------------------------------

         sborntt= cdabs(abornt)**2d0
         sbornbb= cdabs(abornb)**2d0
         sborntb= cdabs(abornt+abornb)**2d0
         sborntb = sborntb-sborntt-sbornbb
C------------------------------------------------
         
         pref=35.0309d0*alphas**2*sborntt

         

         
         ndim=2 


         if (pset.eq.2) then 
            ncomp =3*(Ngrids+1)
            call integrate(ndim,ncomp,INTE2)
            ival=0 
            do iii=0,Ngrids
               ival =ival+1
               qcd3(iii+1) = integral(ival)*pref
               qcd_nnlo(iii+1) = qcd3(iii+1)
            enddo
            do iii=0,Ngrids
               ival =ival+1
               qcd2(iii+1) = integral(ival)*pref
            enddo
            do iii=0,Ngrids
               ival =ival+1
               qcd1(iii+1) = integral(ival)*pref
            enddo

            do iii=0,Ngrids 
               qcd2(iii+1) = qcd2(iii+1)-qcd1(iii+1)
               qcd3(iii+1) = qcd3(iii+1)-qcd2(iii+1)-qcd1(iii+1)
            enddo           
         endif
         
         if (pset.eq.1) then 
            ncomp =2*(Ngrids+1)
            call integrate(ndim,ncomp,INTE1)
            ival=0 
            do iii=0,Ngrids
               qcd3(iii+1) = 0d0
            enddo
            do iii=0,Ngrids
               ival =ival+1
               qcd2(iii+1) = integral(ival)*pref
               qcd_nlo(iii+1) = qcd2(iii+1)
            enddo
            do iii=0,Ngrids
               ival =ival+1
               qcd1(iii+1) = integral(ival)*pref
            enddo
            do iii=0,Ngrids 
               qcd2(iii+1) = qcd2(iii+1)-qcd1(iii+1)
            enddo           
         endif

         if (pset.eq.0) then 
            ncomp =(Ngrids+1)
            call integrate(ndim,ncomp,INTE0)
            ival=0 
            do iii=0,Ngrids
               qcd3(iii+1) = 0d0
            enddo
            do iii=0,Ngrids
               qcd2(iii+1) = 0d0
            enddo
            do iii=0,Ngrids
               ival =ival+1
               qcd1(iii+1) = integral(ival)*pref
               qcd_lo(iii+1) = qcd1(iii+1)
            enddo
         endif

C         do iii=0,Ngrids 
C            print*, pset, qcd1(iii+1), qcd2(iii+1), qcd3(iii+1)
C         enddo      

C----    lambda ->lambda/3 forSM4gen
         
         lambda=lambda/3d0
C----

      do iii=1,(Ngrids+1)
         weak1(iii) = 2d0*lambda*qcd1(iii)
         weak2(iii) = 2d0*lambda*(qcd2(iii)
     &             +alphas*qcd1(iii)*(C1w-C1q))
         
         if (pset.eq.0) then 
            gen4qcd(iii) = 0d0
            weak_lo(iii) =weak1(iii)
         elseif (pset.eq.1) then
            gen4qcd(iii) = 0d0
            weak_nlo(iii)=weak1(iii)+weak2(iii)
         elseif (pset.eq.2) then 
            gen4qcd(iii) = alphas**2*qcd1(iii)*C2qdiff
            weak_nnlo_known(iii)=weak1(iii)+weak2(iii)
            weak3(iii) = 2d0*lambda*(qcd3(iii)
     &                    +alphas*qcd2(iii)*(C1w-C1q))
     #+2d0*lambda*(alphas**2)*qcd1(iii)*(C2w-C2q+C1q*(C1q-C1w))
            weak_nnlo_m5(iii)=weak1(iii)+weak2(iii)+weak3(iii)
         endif
      enddo
      enddo
      
      do iii=1,(Ngrids+1)
         qcd_bot(iii) = qcd_lo(iii)/sborntt*(Ktb*sborntb+Kbb*sbornbb)
         qcd_pure(iii) = qcd_nnlo(iii)+qcd_bot(iii)
         sigma_all(iii) = qcd_pure(iii)+1d0*weak_nnlo_m5(iii)
     &     +1d0*gen4qcd(iii)                     
      enddo

      call givevalerr(Ngrids,qcd_lo,qcd_lo_v,qcd_lo_p,qcd_lo_m)
      call givevalerr(Ngrids,qcd_nlo,qcd_nlo_v,qcd_nlo_p,qcd_nlo_m)
      call givevalerr(Ngrids,qcd_nnlo,qcd_nnlo_v,qcd_nnlo_p,
     &                                                    qcd_nnlo_m)
      call givevalerr(Ngrids,qcd_pure,qcd_pure_v,qcd_pure_p,
     &                                                    qcd_pure_m)
      call givevalerr(Ngrids,sigma_all,sigma_all_v,sigma_all_p,
     &                                                    sigma_all_m)


      print*, "-----------------------"
      print*, collider 
      print*,"mhiggs,mur,muf =",mhiggs,mur,muf
      print*,"mtop,mbot =",mtop,mbot
      print*," cross-section with pdf-error ( +xx% -yy%)"
      print*,"     LO =", qcd_lo_v," ( +",qcd_lo_p,"% -",qcd_lo_m,"%)" 
      print*,"    NLO =", qcd_nlo_v,
     &                    " ( +",qcd_nlo_p,"% -",qcd_nlo_m,"%)"       
      print*,"   NNLO =", qcd_nnlo_v,
     &                    " ( +",qcd_nnlo_p,"% -",qcd_nnlo_m,"%)" 
      print*,"NNLO+bot=", qcd_pure_v,
     &                    " ( +",qcd_pure_p,"% -",qcd_pure_m,"%)" 
      print*,"   sigma=", sigma_all_v,
     &                   " ( +",sigma_all_p,"% -",sigma_all_m,"%)" 

      open(20,FILE='results.dat',STATUS='OLD')
      write(20, *), mhiggs, mur , m4up, m4down, sigma_all_v,
     &                  sigma_all_p, sigma_all_m  
      close(20)

      enddo
      enddo
      enddo
      END


