*---------------------------------------------------------------*
* This program calculates  e^- e^+ -> Higgs gamma. The notation *
* and other explanations can be  found at the page:             *
*                                                               *
*       http://porthos.ist.utl.pt/CTQFT                         *  
*      (See also Nuclear Physics B267(1986)509)                 *
*                                                               *
*    Needs: LoopTools, CUBA integration package, gauss.f for    * 
*    gaussian integration.                                      *
*                                                               *
*    Author: Jorge Crispim Romao                                *
*    email:jorge.romao@ist.utl.pt                               *
*    Last Change:20/10/2008                                     *
*---------------------------------------------------------------*

      program eEHG
      implicit none


**************  Standard Model Physical Constants ******
*
      INCLUDE 'smpc.h'
*
********************************************************

*
*** Variables in the program
*
      integer i,j
      real*8 lbd2    ! For LoopTools
      real*8 g,G_F,mW,mZ,thetaW,mH,GW,GZ,pi
      real*8 ee,alpha,alphaz,factor
      real*8 mlep(3),mup(3),mdown(3)
      real*8 mH2,mW2,mZ2,s,kpm,Egamma
      real*8 mf(9),Qf(9),gv(9),ga(9),sw2,cw2,tw2,eps
      real*8 al(4),ar(4)
      real*8 ReGp1,ImGp1,ReGp2,ImGp2,ReGp3,ImGp3,ReGp4,ImGp4
      real*8 ReGm3,ImGm3,ReGm4,ImGm4
      real*8 xseEHG,rs,step
      complex*16 Gp1,Gp2,Gp3,Gp4,Gm3,Gm4
      complex*16 c1GHG,c1GHGpaper,d1ZHG,d1ZHGpaper
      complex*16 cFGHG,cFGHGpaper,cFZHG,cFZHGpaper


      common/data/mW2,mH2,mZ2,qf,gv,mf,tw2,eps
      common/data2/MZ,GZ
      common/data3/ee,g,pi,factor
      common/data4/al,ar
      common/data5/s

*
*** Initialize LoopTools
*
#include "looptools.h"

      call ffini

      lbd2=1.d0
      call setlambda(lbd2)

      open(10,file="xseEHG.dat",status="unknown")

*
***  1GeV^(-2)= 3.893795 10^(8) pb (Conversion factor)
*

      factor=3.893795d+08

      pi=acos(-1.d0)
      mW=SM_MW
      GW=SM_GAW
      mZ=SM_MZ
      GZ=SM_GAZ
      G_F=SM_GFERMI
      mf(1)=SM_EMASS
      mf(2)=SM_MUMASS
      mf(3)=SM_TAUMASS
      mf(4)=SM_UPMASS
      mf(5)=SM_CHARMMASS
      mf(6)=SM_TOPMASS
      mf(7)=SM_DOWNMASS
      mf(8)=SM_STRANGEMASS
      mf(9)=SM_BOTTOMMASS
      alpha=SM_ALPHA
      alphaz=SM_ALPHAZ
      thetaW=asin(sqrt(SM_S2TW))
      sw2=SM_S2TW


*
*** Data at the time of paper
*
      mf(6)=40d0
      mW=80.8d0
      mZ=92.9d0
      sw2=0.21d0
      thetaW=asin(sqrt(sw2))
      tw2=tan(thetaW)**2
      cw2=cos(thetaW)**2
      ee=sqrt(4*pi*alpha)
      g=ee/sin(thetaW)
      GZ=2.953d0

      do i=1,3
         qf(i)=-1
         qf(3+i)=2d0/3d0
         qf(6+i)=-1d0/3d0
      enddo

      do i=1,3
         gv(i)=-1d0/4d0-qf(i)*sw2
         ga(i)=-1d0/4d0
         gv(i+3)=1d0/4d0-qf(i+3)*sw2
         ga(i+3)=1d0/4d0
         gv(i+6)=-1d0/4d0-qf(i+6)*sw2
         ga(i+6)=-1d0/4d0
      enddo

      al(1)=sw2
      al(2)=gv(1)+ga(1)
      ar(2)=gv(1)-ga(1)
      al(3)=1d0
      al(4)=(gv(1)+ga(1))**2
      ar(4)=(gv(1)-ga(1))**2
      
      mW2=mW*mW
      mZ2=mZ*mZ
      call setmudim(mW2)

      write(*,*)mW,mZ,GZ,sw2

      mH=60d0
      mH2=mH*mH
      eps=1d-10
      rs=74d0
      step=0.2d0

 10   if(rs.gt.150d0)goto 20
      s=rs*rs
      write(10,100)rs,xseEHG(s)
      rs=rs+step
      goto 10

 20   if(.FALSE.)then
         s=mZ2
         rs=mZ
         write(10,100)rs,xseEHG(s)
      endif

      close(10)

 98   format(f8.2,4x,2("(",1p,e12.5,",",1p,e12.5,")",2x))

 100  format(f8.2,4x,1p,e14.5)


      end
*-------------------------------------------------------------------*
      real*8 function xseEHG(s)
      implicit none
      real*8 s
      real*8 cst,ee,g,pi,factor,res
      real*8 xmin(1),xmax(1)
      real*8 mW2,mH2,mZ2,qf(9),gv(9),mf(9),tw2,eps

      real*8 ReGp1,ImGp1,ReGp2,ImGp2
      real*8 ReGp3,ImGp3,ReGp4,ImGp4,ReGm3,ImGm3,ReGm4,ImGm4
      real*8 Eg,kpm,sw,alfa


      common/boundary/xmin,xmax
      common/data/mW2,mH2,mZ2,qf,gv,mf,tw2,eps
      common/data3/ee,g,pi,factor


      cst=1d0/(64d0*pi*s)*(s-mH2)/(2d0*s)
      cst=cst*(ee*g**3/(16d0*pi**2))**2
      cst=cst*factor

      xmin(1)=-1d0
      xmax(1)= 1d0

      call Integration(res)

      xseEHG=res*cst

      return

      Call LoopFunctions(s,kpm,ReGp1,ImGp1,ReGp2,ImGp2,
     &     ReGp3,ImGp3,ReGp4,ImGp4,ReGm3,ImGm3,ReGm4,ImGm4)

      sw=sqrt(1d0-1d0/(1d0+tw2))
      alfa=ee**2/(4d0*pi)

      Eg=(s-mH2)/2d0/sqrt(s)
      xseEHG=1d0/(12d0*pi)*(alfa)**4/(sw)**6
      xseEHG=xseEHG*(gv(1)**2+0.25d0**2)*(Eg*s)**2/mW2**3
      xseEHG=xseEHG*Eg/s/sqrt(s)*factor*(ReGp2**2+ImGp2**2)

      return
      end
     

      subroutine Integration(res)
      implicit none
      real*8 res
*
**   Modify the parameter method to choose the method of integration. 
*
      integer i,j,npoints,method
      parameter (npoints = 16)
      parameter (method = 1)
*
***   Variables for all the methods
*
      integer ndim,ncomp
      parameter (ndim = 1)
      parameter (ncomp = 1)
      double precision integral(ncomp)

*
***   Other Variables for vegas,suave,divonne and cuhre
*

      integer verbose,mineval,maxeval
      parameter (verbose = 0)
      parameter (mineval = 0)
      parameter (maxeval = 50000)

      double precision epsrel, epsabs
      parameter (epsrel = 1D-3)
      parameter (epsabs = 1D-12)

      double precision error(ncomp), prob(ncomp)
      integer nregions, neval, fail
      

*
***   Variables only in vegas
*

      integer nstart, nincrease
      parameter (nstart = 5000)
      parameter (nincrease = 500)

      external integrand

      if(method.eq.1)then
         call gauss(ndim, ncomp, npoints, integrand, integral, fail)
         if(fail.eq.1)then
            write(*,*)'Gauss method only works for ndim <= 6'
            stop
         endif

         if(fail.eq.2)then
            write(*,*)'Gauss method only works for
     &           npoints=4,6,8,10,16,20,24,32,64,128'
            stop
         endif
         
      endif

      if(method.eq.2)then
         call vegas(ndim, ncomp, integrand,
     &        epsrel, epsabs, verbose, mineval, maxeval,
     &        nstart, nincrease,
     &        neval, fail, integral, error, prob)
      endif

      res=integral(1)
      
      return
      end

*-------------------------------------------------------------------*
*
*** The integration has to be mapped from the interval [0,1] into
*** [xmin,xmax] 
*
*-------------------------------------------------------------------*
      subroutine integrand(ndim, xx, ncomp, ff)
      implicit none
      integer i, ndim, ncomp
      double precision xx(*), ff(*)
      double precision x(ndim),xmin(1),xmax(1),Jacob
      double precision dfdx

      common/boundary/xmin,xmax

      Jacob=1.d0

      do i=1,ndim
         Jacob=Jacob*(xmax(i)-xmin(i))
         x(i)=xmin(i)+(xmax(i)-xmin(i))*xx(i)
      enddo
      
      ff(1) = dfdx(x(1))*Jacob

      return
      end

*-------------------------------------------------------------------*
*
*** This is the function to integrate. 
*
*-------------------------------------------------------------------*
      double precision function dfdx(x)
      implicit none
      integer idiag
      real*8 mW2,mH2,mZ2,qf(9),gv(9),mf(9),tw2,eps
      real*8 s,x,kpm,Egamma,MSquared
      real*8 al(4),ar(4)
      real*8 ReGp1,ImGp1,ReGp2,ImGp2
      real*8 ReGp3,ImGp3,ReGp4,ImGp4,ReGm3,ImGm3,ReGm4,ImGm4
      real*8 TG,TZ,TGTZ,BW,BZ,BWBZ,TGZBWZ
      real*8 MSquared2

      common/data/mW2,mH2,mZ2,qf,gv,mf,tw2,eps
      common/data4/al,ar
      common/data5/s


      Egamma=(s-MH2)/2d0/sqrt(s)
      kpm=sqrt(s)/2d0*Egamma*(1d0-x)

      Call LoopFunctions(s,kpm,ReGp1,ImGp1,ReGp2,ImGp2,
     &     ReGp3,ImGp3,ReGp4,ImGp4,ReGm3,ImGm3,ReGm4,ImGm4)
        

      idiag=0

        TG=((ImGp1**2 + ReGp1**2)*
     &    (8*kpm**2 + 4*kpm*(mH2 - s) + (mH2-s)**2)*s*al(1)**2)/mW2**3

        TZ=((ImGp2**2 + ReGp2**2)*
     &     (8*kpm**2 + 4*kpm*(mH2 - s) + (mH2 - s)**2)*s*(al(2)**2 
     &     + ar(2)**2))/(2.*mW2**3)

        TGTZ=((ImGp1*ImGp2 + ReGp1*ReGp2)*
     &     (8*kpm**2 + 4*kpm*(mH2 - s) + (mH2 - s)**2)*s*al(1)*(al(2) 
     &       + ar(2)))/mW2**3

        BW=((4*kpm**2*(ImGm3**2 + ImGp3**2 + ReGm3**2 + ReGp3**2) + 
     &       4*kpm*(ImGp3**2 + ReGp3**2)*(mH2 - s) + 
     &       (ImGp3**2 + ReGp3**2)*(mH2 - s)**2)*s*al(3)**2)/(2.*mW2**3)

        BZ=(s*(4*ImGm4**2*kpm**2 + 4*kpm**2*ReGm4**2 + 4*kpm**2*ReGp4**2 
     &       + 4*kpm*mH2*ReGp4**2 + mH2**2*ReGp4**2 + 
     &       ImGp4**2*(2*kpm + mH2 - s)**2 - 4*kpm*ReGp4**2*s - 
     &       2*mH2*ReGp4**2*s + ReGp4**2*s**2)*(al(4)**2 + ar(4)**2))/
     &   (2.*mW2**3)

        BWBZ=(s*(4*ImGm3*ImGm4*kpm**2 + 4*kpm**2*ReGm3*ReGm4 + 
     &       4*kpm**2*ReGp3*ReGp4 + 4*kpm*mH2*ReGp3*ReGp4 + 
     &       mH2**2*ReGp3*ReGp4 + ImGp3*ImGp4*(2*kpm + mH2 - s)**2 - 
     &       4*kpm*ReGp3*ReGp4*s - 2*mH2*ReGp3*ReGp4*s 
     &       + ReGp3*ReGp4*s**2)*al(3)*al(4))/mW2**3

        TGZBWZ=(s*(4*kpm**2*ReGm3*ReGp1*al(1)*al(3) + 
     &       4*kpm**2*ReGp1*ReGp3*al(1)*al(3) + 
     &       4*kpm*mH2*ReGp1*ReGp3*al(1)*al(3) + 
     &       mH2**2*ReGp1*ReGp3*al(1)*al(3) - 
     &       4*kpm*ReGp1*ReGp3*s*al(1)*al(3) - 
     &       2*mH2*ReGp1*ReGp3*s*al(1)*al(3) + 
     &       ReGp1*ReGp3*s**2*al(1)*al(3) + 
     &       4*ImGp2*ImGp3*kpm**2*al(2)*al(3) + 
     &       4*ImGp2*ImGp3*kpm*mH2*al(2)*al(3) + 
     &       ImGp2*ImGp3*mH2**2*al(2)*al(3) + 
     &       4*kpm**2*ReGm3*ReGp2*al(2)*al(3) + 
     &       4*kpm**2*ReGp2*ReGp3*al(2)*al(3) + 
     &       4*kpm*mH2*ReGp2*ReGp3*al(2)*al(3) + 
     &       mH2**2*ReGp2*ReGp3*al(2)*al(3) - 
     &       4*ImGp2*ImGp3*kpm*s*al(2)*al(3) - 
     &       2*ImGp2*ImGp3*mH2*s*al(2)*al(3) - 
     &       4*kpm*ReGp2*ReGp3*s*al(2)*al(3) - 
     &       2*mH2*ReGp2*ReGp3*s*al(2)*al(3) + 
     &       ImGp2*ImGp3*s**2*al(2)*al(3) + 
     &       ReGp2*ReGp3*s**2*al(2)*al(3) + 
     &       4*ImGm3*kpm**2*(ImGp1*al(1) + ImGp2*al(2))*al(3) + 
     &       4*kpm**2*ReGm4*ReGp1*al(1)*al(4) + 
     &       4*kpm**2*ReGp1*ReGp4*al(1)*al(4) + 
     &       4*kpm*mH2*ReGp1*ReGp4*al(1)*al(4) + 
     &       mH2**2*ReGp1*ReGp4*al(1)*al(4) - 
     &       4*kpm*ReGp1*ReGp4*s*al(1)*al(4) - 
     &       2*mH2*ReGp1*ReGp4*s*al(1)*al(4) + 
     &       ReGp1*ReGp4*s**2*al(1)*al(4) + 
     &       4*ImGm4*ImGp2*kpm**2*al(2)*al(4) + 
     &       4*ImGp2*ImGp4*kpm**2*al(2)*al(4) + 
     &       4*ImGp2*ImGp4*kpm*mH2*al(2)*al(4) + 
     &       ImGp2*ImGp4*mH2**2*al(2)*al(4) + 
     &       4*kpm**2*ReGm4*ReGp2*al(2)*al(4) + 
     &       4*kpm**2*ReGp2*ReGp4*al(2)*al(4) + 
     &       4*kpm*mH2*ReGp2*ReGp4*al(2)*al(4) + 
     &       mH2**2*ReGp2*ReGp4*al(2)*al(4) - 
     &       4*ImGp2*ImGp4*kpm*s*al(2)*al(4) - 
     &       2*ImGp2*ImGp4*mH2*s*al(2)*al(4) - 
     &       4*kpm*ReGp2*ReGp4*s*al(2)*al(4) - 
     &       2*mH2*ReGp2*ReGp4*s*al(2)*al(4) 
     &       + ImGp2*ImGp4*s**2*al(2)*al(4) + 
     &       ReGp2*ReGp4*s**2*al(2)*al(4) + 
     &       4*kpm**2*ReGm4*ReGp1*al(1)*ar(4) + 
     &       4*kpm**2*ReGp1*ReGp4*al(1)*ar(4) + 
     &       4*kpm*mH2*ReGp1*ReGp4*al(1)*ar(4) + 
     &       mH2**2*ReGp1*ReGp4*al(1)*ar(4) - 
     &       4*kpm*ReGp1*ReGp4*s*al(1)*ar(4) - 
     &       2*mH2*ReGp1*ReGp4*s*al(1)*ar(4) + 
     &       ReGp1*ReGp4*s**2*al(1)*ar(4) + 
     &       4*ImGm4*ImGp2*kpm**2*ar(2)*ar(4) + 
     &       4*ImGp2*ImGp4*kpm**2*ar(2)*ar(4) + 
     &       4*ImGp2*ImGp4*kpm*mH2*ar(2)*ar(4) + 
     &       ImGp2*ImGp4*mH2**2*ar(2)*ar(4) + 
     &       4*kpm**2*ReGm4*ReGp2*ar(2)*ar(4) + 
     &       4*kpm**2*ReGp2*ReGp4*ar(2)*ar(4) + 
     &       4*kpm*mH2*ReGp2*ReGp4*ar(2)*ar(4) + 
     &       mH2**2*ReGp2*ReGp4*ar(2)*ar(4) - 
     &       4*ImGp2*ImGp4*kpm*s*ar(2)*ar(4) - 
     &       2*ImGp2*ImGp4*mH2*s*ar(2)*ar(4) - 
     &       4*kpm*ReGp2*ReGp4*s*ar(2)*ar(4) - 
     &       2*mH2*ReGp2*ReGp4*s*ar(2)*ar(4) + 
     &       ImGp2*ImGp4*s**2*ar(2)*ar(4) + 
     &       ReGp2*ReGp4*s**2*ar(2)*ar(4) + 
     &       ImGp1*al(1)*(ImGp3*(2*kpm + mH2 - s)**2*al(3) + 
     &       (4*ImGm4*kpm**2 + ImGp4*(2*kpm + mH2 - s)**2)*(al(4) + 
     &       ar(4)))))/mW2**3


      if(idiag.eq.0)then
         MSquared=TG+TZ+TGTZ+BW+BZ+BWBZ+TGZBWZ
      endif
      if(idiag.eq.1)then
         MSquared=TZ
      endif
      if(idiag.eq.2)then
         MSquared=TG
      endif
      if(idiag.eq.3)then
         MSquared=TGTZ
      endif
      if(idiag.eq.4)then
         MSquared=BW+BZ+BWBZ
      endif
      if(idiag.eq.5)then
         MSquared=TGZBWZ
      endif

      dfdx=MSquared


      return
      end
*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      Subroutine LoopFunctions(s,kpm,ReGp1,ImGp1,ReGp2,ImGp2,
     &        ReGp3,ImGp3,ReGp4,ImGp4,ReGm3,ImGm3,ReGm4,ImGm4)
      implicit none
      real*8 s,kpm,ReGp1,ImGp1,ReGp2,ImGp2
      real*8 ReGp3,ImGp3,ReGp4,ImGp4,ReGm3,ImGm3,ReGm4,ImGm4
      complex*16 aux,Gp1,Gp2,Gp3,Gp4,Gm3,Gm4

      aux=Gp1(s)
      ReGp1=dble(aux)
      ImGp1=dimag(aux)
      aux=Gp2(s)
      ReGp2=dble(aux)
      ImGp2=dimag(aux)
      aux=Gp3(s,kpm)
      ReGp3=dble(aux)
      ImGp3=dimag(aux)
      aux=Gp4(s,kpm)
      ReGp4=dble(aux)
      ImGp4=dimag(aux)

      aux=Gm3(s,kpm)
      ReGm3=dble(aux)
      ImGm3=dimag(aux)
      aux=Gm4(s,kpm)
      ReGm4=dble(aux)
      ImGm4=dimag(aux)

      return

      end
*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function Gp1(s)
      implicit none
      real*8 s,mW2,mH2,mZ2,Qf(9),gv(9),mf(9),tw2,eps
      complex*16 CFGHG,c1GHG

      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps

      Gp1=MW2/s*(CFGHG(s)+c1GHG(s))

      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function Gp2(s)
      implicit none
      real*8 s,mW2,mH2,mZ2,Qf(9),gv(9),mf(9),tw2,eps
      real*8 MZ,GZ
      complex*16 DZ2,CFZHG,d1ZHG,aux

      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps
      common/data2/MZ,GZ


      aux=dcmplx(MZ,-0.5d0*GZ)
      
      DZ2=s-aux**2

      Gp2=MW2/DZ2*(CFZHG(s)+d1ZHG(s))

      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function Gp3(s,kpm)
      implicit none
      real*8 mW4
      real*8 s,kpm,mW2,mH2,mZ2,Qf(9),gv(9),mf(9),tw2,eps
      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps

*
*** Initialize LoopTools
*
#include "looptools.h"

      mW4=mW2**2

      
        Gp3=mW4*((2*B0(-2*kpm,eps,mW2))/
     &  ((2*kpm + mH2)*(2*kpm + mH2 - s)) - 
     & (2*B0(mH2,mW2,mW2))/((2*kpm + mH2)*(2*kpm + mH2 - s)) - 
     & (2*kpm*C0(0,0,-2*kpm,eps,mW2,mW2))/(2*kpm + mH2 - s)**2 + 
     & (2*kpm*mW2*C0(0,0,-2*kpm,eps,mW2,mW2))/
     &  ((2*kpm + mH2 - s)**2*s) - 
     & (2*kpm*C0(0,0,-2*kpm,eps,mW2,mW2))/((2*kpm + mH2 - s)*s) + 
     & C0(0,0,2*kpm + mH2 - s,eps,mW2,mW2)/s + 
     & (mW2*C0(0,0,2*kpm + mH2 - s,eps,mW2,mW2))/
     &  (s*(-2*kpm - mH2 + s)) + 
     & (s*C0(0,0,s,mW2,eps,mW2))/(2*kpm + mH2 - s)**2 - 
     & (2*C0(0,-2*kpm,mH2,mW2,eps,mW2))/(2*kpm + mH2) + 
     & (2*mW2*C0(0,-2*kpm,mH2,mW2,eps,mW2))/
     &  ((2*kpm + mH2)*(2*kpm + mH2 - s)) + 
     & (mW2*s*C0(0,-2*kpm,mH2,mW2,eps,mW2))/
     &  ((2*kpm + mH2)*(2*kpm + mH2 - s)**2) - 
     & (3*s*C0(0,-2*kpm,mH2,mW2,eps,mW2))/
     &  ((2*kpm + mH2)*(2*kpm + mH2 - s)) - 
     & (s**2*C0(0,-2*kpm,mH2,mW2,eps,mW2))/
     &  ((2*kpm + mH2)*(2*kpm + mH2 - s)**2) + 
     & (2*kpm*C0(0,-2*kpm,mH2,mW2,eps,mW2))/(2*kpm*s + mH2*s) + 
     & (mH2*C0(0,-2*kpm,mH2,mW2,eps,mW2))/(2*kpm*s + mH2*s) - 
     & (mW2*C0(0,-2*kpm,mH2,mW2,eps,mW2))/(2*kpm*s + mH2*s) - 
     & (4*kpm*mW2*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     &  ((2*kpm - s)*(2*kpm + mH2 - s)**2) + 
     & C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2)/(2*kpm + mH2 - s) + 
     & (kpm*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     &  ((2*kpm - s)*(2*kpm + mH2 - s)) + 
     & (4*kpm**2*mW2*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     &  ((2*kpm - s)*(2*kpm + mH2 - s)**2*s) + 
     & (mW2*s*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     &  ((2*kpm - s)*(2*kpm + mH2 - s)**2) + 
     & (s*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     &  (2.*(2*kpm - s)*(2*kpm + mH2 - s)) - 
     & (s**2*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     &  (4.*kpm*(2*kpm - s)*(2*kpm + mH2 - s)) + 
     & (kpm*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     &  (s*(-2*kpm - mH2 + s)) + 
     & (kpm**2*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     &  ((kpm - s/2.)*s*(-2*kpm - mH2 + s)) - 
     & (s*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     &  (8*kpm**2 + 4*kpm*mH2 - 4*kpm*s) + 
     & (2*kpm*C0(0,s,mH2,mW2,mW2,mW2))/(2*kpm + mH2 - s)**2 - 
     & (2*C0(0,s,mH2,mW2,mW2,mW2))/s - 
     & (4*kpm*mW2*C0(0,s,mH2,mW2,mW2,mW2))/
     &  ((2*kpm + mH2 - s)**2*s) + 
     & (4*kpm*C0(0,s,mH2,mW2,mW2,mW2))/((2*kpm + mH2 - s)*s) + 
     & (2*mW2*C0(0,s,mH2,mW2,mW2,mW2))/((2*kpm + mH2 - s)*s) + 
     & C0(0,s,mH2,mW2,mW2,mW2)/(-2*kpm - mH2 + s) - 
     & (4*kpm*mW2*D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/
     &  (2*kpm + mH2 - s)**2 + 
     & (2*kpm*D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/
     &  (2*kpm + mH2 - s) + 
     & (mW2*D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/s + 
     & (2*kpm*mW4*D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/
     &  ((2*kpm + mH2 - s)**2*s) - 
     & (2*kpm*mW2*D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/
     &  ((2*kpm + mH2 - s)*s) + 
     & (2*kpm*s*D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/
     &  (2*kpm + mH2 - s)**2 + 
     & (mW2*D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/
     &  (-2*kpm - mH2 + s) + 
     & (mW4*D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/
     &  (s*(-2*kpm - mH2 + s)) + 
     & D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,mW2) + 
     & (mW2*D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,mW2))/
     &  s + (2*kpm*mW4*D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,
     &  eps,mW2))/((2*kpm + mH2 - s)**2*s) - 
     & (2*kpm*mW2*D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,
     &  mW2))/((2*kpm + mH2 - s)*s) + 
     & (mW2*D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,mW2))/
     &  (-2*kpm - mH2 + s) + 
     & (mW4*D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,mW2))/
     &  (s*(-2*kpm - mH2 + s)))

      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function Gp4(s,kpm)
      implicit none
      real*8 mZ4
      real*8 s,kpm,mW2,mH2,mZ2,Qf(9),gv(9),mf(9),tw2,eps

      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps

*
*** Initialize LoopTools
*
#include "looptools.h"

      mZ4=mZ2**2

        Gp4=(2*mZ4*((4*kpm*B0(-2*kpm,eps,mZ2))/(2*kpm + mH2) + 
     &   (2*mH2*B0(-2*kpm,eps,mZ2))/(2*kpm + mH2) - 
     &   (2*s*B0(-2*kpm,eps,mZ2))/(2*kpm + mH2) - 
     &   (4*kpm*B0(mH2,mZ2,mZ2))/(2*kpm + mH2) - 
     &   (2*mH2*B0(mH2,mZ2,mZ2))/(2*kpm + mH2) + 
     &   (2*s*B0(mH2,mZ2,mZ2))/(2*kpm + mH2) - 
     &   2*kpm*C0(0,0,-2*kpm,eps,eps,mZ2) + 
     &   (4*kpm**2*C0(0,0,-2*kpm,eps,eps,mZ2))/s + 
     &   (2*kpm*mH2*C0(0,0,-2*kpm,eps,eps,mZ2))/s - 
     &   (2*kpm*mZ2*C0(0,0,-2*kpm,eps,eps,mZ2))/s + 
     &   4*kpm*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2) + 
     &   2*mH2*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2) - 
     &   mZ2*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2) - 
     &   (4*kpm**2*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2))/s - 
     &   (4*kpm*mH2*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2))/s - 
     &   (mH2**2*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2))/s + 
     &   (2*kpm*mZ2*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2))/s + 
     &   (mH2*mZ2*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2))/s - 
     &   s*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2) - 
     &   (12*kpm**2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/(2*kpm + mH2) - 
     &   (12*kpm*mH2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/(2*kpm + mH2) - 
     &   (3*mH2**2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/(2*kpm + mH2) + 
     &   (8*kpm*mZ2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/
     & (2*kpm + mH2) + (4*mH2*mZ2*
     &   C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/(2*kpm + mH2) + 
     &   (8*kpm**3*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/
     & ((2*kpm + mH2)*s) + 
     &   (12*kpm**2*mH2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/
     & ((2*kpm + mH2)*s) + 
     &   (6*kpm*mH2**2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/
     & ((2*kpm + mH2)*s) + 
     &   (mH2**3*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/((2*kpm + mH2)*s) - 
     &   (4*kpm**2*mZ2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/
     & ((2*kpm + mH2)*s) - 
     &   (4*kpm*mH2*mZ2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/
     & ((2*kpm + mH2)*s) - 
     &   (mH2**2*mZ2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/
     & ((2*kpm + mH2)*s) + 
     &   (4*kpm*s*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/(2*kpm + mH2) + 
     &   (2*mH2*s*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/(2*kpm + mH2) - 
     &   (2*mZ2*s*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/(2*kpm + mH2) + 
     &   4*kpm*C0(0,2*kpm + mH2 - s,mH2,mZ2,eps,mZ2) + 
     &   mH2*C0(0,2*kpm + mH2 - s,mH2,mZ2,eps,mZ2) - 
     &   mZ2*C0(0,2*kpm + mH2 - s,mH2,mZ2,eps,mZ2) - 
     &   (4*kpm**2*C0(0,2*kpm + mH2 - s,mH2,mZ2,eps,mZ2))/s - 
     &   (2*kpm*mH2*C0(0,2*kpm + mH2 - s,mH2,mZ2,eps,mZ2))/s + 
     &   (2*kpm*mZ2*C0(0,2*kpm + mH2 - s,mH2,mZ2,eps,mZ2))/s - 
     &   s*C0(0,2*kpm + mH2 - s,mH2,mZ2,eps,mZ2) + 
     &   8*kpm**2*D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,
     &   mZ2) + 4*kpm*mH2*
     &   D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,mZ2)
     &   + 2*mH2*mZ2*D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,
     &   eps,mZ2) - mZ4*
     &   D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,mZ2)
     &   - (8*kpm**3*D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,
     &    eps,mZ2))/s - 
     &   (8*kpm**2*mH2*D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,
     &    eps,mZ2))/s - 
     &   (2*kpm*mH2**2*D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,
     &    eps,mZ2))/s + 
     &   (4*kpm**2*mZ2*D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,
     &    eps,mZ2))/s - 
     &   (mH2**2*mZ2*D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,
     &    eps,mZ2))/s + 
     &   (mH2*mZ4*D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,
     &    eps,mZ2))/s - 
     &   2*kpm*s*D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,
     &  mZ2) - mZ2*s*
     &  D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,mZ2)))
     &  /(2*kpm + mH2 - s)**2
    
      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function Gm3(s,kpm)
      implicit none
      
      real*8 s,kpm,mW4,mZ4
      real*8 mW2,mH2,mZ2,Qf(9),gv(9),mf(9),tw2,eps

      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps

*
*** Initialize LoopTools
*
#include "looptools.h"

      mW4=mW2*mW2

        Gm3=(mW4*((-8*(2*kpm + mH2 - s)*B0(mH2,mW2,mW2))/
     & (kpm*(2*kpm + mH2)*(2*kpm - s)) - 
     &   (8*s*B0(mH2,mW2,mW2))/(kpm*(2*kpm + mH2)*(2*kpm - s)) + 
     &   (8*B0(2*kpm + mH2 - s,eps,mW2))/(2*kpm**2 - kpm*s) + 
     &   (8*C0(0,0,-2*kpm,eps,mW2,mW2))/s + 
     &   (4*mW2*C0(0,0,-2*kpm,eps,mW2,mW2))/(kpm*s) + 
     &   (2*(2*kpm + mH2 - s)*C0(0,0,2*kpm + mH2 - s,eps,mW2,mW2))/
     & kpm**2 - (4*(2*kpm + mH2 - s)*
     &   C0(0,0,2*kpm + mH2 - s,eps,mW2,mW2))/(kpm*s) - 
     &   (2*mW2*(2*kpm + mH2 - s)*
     &   C0(0,0,2*kpm + mH2 - s,eps,mW2,mW2))/(kpm**2*s) + 
     &   (2*s*C0(0,0,s,mW2,eps,mW2))/kpm**2 - 
     &   (4*C0(0,-2*kpm,mH2,mW2,eps,mW2))/kpm - 
     &   (2*(2*kpm + mH2 - s)*C0(0,-2*kpm,mH2,mW2,eps,mW2))/
     & (kpm*(2*kpm + mH2)) - 
     &   (4*mW2*(2*kpm + mH2 - s)*C0(0,-2*kpm,mH2,mW2,eps,mW2))/
     & (kpm**2*(2*kpm + mH2)) - 
     &   (2*(2*kpm + mH2 - s)*C0(0,-2*kpm,mH2,mW2,eps,mW2))/
     & (kpm*s) - (2*(2*kpm + mH2 - s)**2*
     &   C0(0,-2*kpm,mH2,mW2,eps,mW2))/(kpm*(2*kpm + mH2)*s) - 
     &   (2*mW2*(2*kpm + mH2 - s)**2*
     &   C0(0,-2*kpm,mH2,mW2,eps,mW2))/(kpm**2*(2*kpm + mH2)*s)
     & + (2*s*C0(0,-2*kpm,mH2,mW2,eps,mW2))/
     & (2*kpm**2 + kpm*mH2) - 
     &   (2*mW2*s*C0(0,-2*kpm,mH2,mW2,eps,mW2))/
     & (kpm**2*(2*kpm + mH2)) - 
     &   (2*s*C0(0,-2*kpm,mH2,mW2,eps,mW2))/
     & (kpm*(2*kpm + mH2 - s)) + 
     &   (2*s**2*C0(0,-2*kpm,mH2,mW2,eps,mW2))/
     & (kpm*(2*kpm + mH2)*(2*kpm + mH2 - s)) + 
     &   (8*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/(2*kpm - s) - 
     &   (2*mW2*s*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     & (kpm**2*(2*kpm - s)) + 
     &   (s**2*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     & (kpm**2*(kpm - s/2.)) + 
     &   (8*mW2*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     & (2*kpm**2 - kpm*s) - 
     &   (12*s*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     & (2*kpm**2 - kpm*s) + 
     &   (16*kpm*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     & (2*kpm*s - s**2) + 
     &   (8*mW2*C0(0,2*kpm + mH2 - s,mH2,mW2,eps,mW2))/
     & (2*kpm*s - s**2) + (4*C0(0,s,mH2,mW2,mW2,mW2))/kpm - 
     &   (2*(2*kpm + mH2 - s)*C0(0,s,mH2,mW2,mW2,mW2))/kpm**2 - 
     &   (16*C0(0,s,mH2,mW2,mW2,mW2))/s - 
     &   (8*mW2*C0(0,s,mH2,mW2,mW2,mW2))/(kpm*s) + 
     &   (8*(2*kpm + mH2 - s)*C0(0,s,mH2,mW2,mW2,mW2))/(kpm*s) + 
     &   (4*mW2*(2*kpm + mH2 - s)*C0(0,s,mH2,mW2,mW2,mW2))/
     & (kpm**2*s) + 8*D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,
     &  mW2) + (4*mW2*
     &   D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/kpm + 
     &   (8*mW2*D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/s + 
     &   (4*mW4*D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/
     & (kpm*s) - (4*mW2*(2*kpm + mH2 - s)*
     &   D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/(kpm*s) - 
     &   (2*mW4*(2*kpm + mH2 - s)*
     &   D0(0,0,0,mH2,-2*kpm,s,mW2,mW2,eps,mW2))/(kpm**2*s) + 
     &   (4*mW2*D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,
     &    mW2))/kpm + (4*(2*kpm + mH2 - s)*
     &   D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,mW2))/kpm
     & + (4*mW2*(2*kpm + mH2 - s)*
     &   D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,mW2))/
     & kpm**2 + (8*mW2*
     &   D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,mW2))/s + 
     &   (4*mW4*D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,
     &    mW2))/(kpm*s) - 
     &   (4*mW2*(2*kpm + mH2 - s)*
     &   D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,mW2))/
     & (kpm*s) - (2*mW4*(2*kpm + mH2 - s)*
     &   D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,mW2))/
     & (kpm**2*s) - (2*(2*kpm + mH2 - s)*s*
     &   D0(0,0,0,mH2,2*kpm + mH2 - s,s,mW2,mW2,eps,mW2))/
     & kpm**2))/8.

      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function Gm4(s,kpm)
      implicit none
      
      real*8 s,kpm,mZ4
      real*8 mW2,mH2,mZ2,Qf(9),gv(9),mf(9),tw2,eps

      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps

*
*** Initialize LoopTools
*
#include "looptools.h"
      

      mZ4=mZ2*mZ2

        Gm4=(mZ4*((-4*kpm*B0(mH2,mZ2,mZ2))/(2*kpm - s) + 
     &   (4*kpm*B0(2*kpm + mH2 - s,eps,mZ2))/(2*kpm - s) - 
     &   (4*kpm**2*C0(0,0,-2*kpm,eps,eps,mZ2))/s - 
     &   (2*kpm*mZ2*C0(0,0,-2*kpm,eps,eps,mZ2))/s - 
     &   2*kpm*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2) - 
     &   mZ2*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2) + 
     &   (4*kpm**2*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2))/s + 
     &   (2*kpm*mH2*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2))/s + 
     &   (2*kpm*mZ2*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2))/s + 
     &   (mH2*mZ2*C0(0,0,2*kpm + mH2 - s,eps,eps,mZ2))/s - 
     &   (4*kpm**2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/s - 
     &   (2*kpm*mH2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/s - 
     &   (2*kpm*mZ2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/s - 
     &   (mH2*mZ2*C0(0,-2*kpm,mH2,mZ2,eps,mZ2))/s + 
     &   (4*kpm*mZ2*C0(0,2*kpm + mH2 - s,mH2,mZ2,eps,mZ2))/
     & (2*kpm - s) + (8*kpm**3*
     &   C0(0,2*kpm + mH2 - s,mH2,mZ2,eps,mZ2))/((2*kpm - s)*s)
     & + (4*kpm**2*mZ2*C0(0,2*kpm + mH2 - s,mH2,mZ2,eps,mZ2))/
     & ((2*kpm - s)*s) - 
     &   (2*kpm*s*C0(0,2*kpm + mH2 - s,mH2,mZ2,eps,mZ2))/
     & (2*kpm - s) - (mZ2*s*
     &   C0(0,2*kpm + mH2 - s,mH2,mZ2,eps,mZ2))/(2*kpm - s) - 
     &   4*kpm**2*D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,
     &  mZ2) - 4*kpm*mZ2*
     & D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,mZ2)
     & - mZ4*D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,
     &  mZ2) + (8*kpm**3*
     &   D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,mZ2)
     &   )/s + (4*kpm**2*mH2*
     &   D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,mZ2)
     &   )/s + (4*kpm**2*mZ2*
     &   D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,mZ2)
     &   )/s + (4*kpm*mH2*mZ2*
     &   D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,mZ2)
     &   )/s + (mH2*mZ4*
     &   D0(0,0,0,mH2,-2*kpm,2*kpm + mH2 - s,mZ2,eps,eps,mZ2)
     &   )/s))/(2.*kpm**2)
    
      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function c1GHG(s)
      implicit none
      real*8 s
      real*8 mW2,mH2,mZ2,Qf(9),gv(9),mf(9),tw2,eps

      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps

*
*** Initialize LoopTools
*
#include "looptools.h"


        c1GHG=(2*(-((mH2 + 6*mW2)*s*B0(mH2,mW2,mW2)) + 
     &   (mH2 + 6*mW2)*s*B0(s,mW2,mW2) + 
     &   (mH2 - s)*(-mH2 - 6*mW2 - 
     &   2*mW2*(-3*mH2 + 6*mW2 + 4*s)*
     &    C0(0,s,mH2,mW2,mW2,mW2))))/(mH2 - s)**2



      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function d1ZHG(s)
      implicit none
      real*8 s
      real*8 mW2,mH2,mZ2,Qf(9),gv(9),mf(9),tw2,eps

      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps
*
*** Initialize LoopTools
*
#include "looptools.h"


        d1ZHG=(s*(-2*mW2*(-5 + tw2) - mH2*(-1 + tw2))*
     &  B0(mH2,mW2,mW2) + 
     & s*(2*mW2*(-5 + tw2) + mH2*(-1 + tw2))*B0(s,mW2,mW2) + 
     & (mH2 - s)*(mH2 + 10*mW2 - mH2*tw2 - 2*mW2*tw2 + 
     & 2*mW2*(-2*(mW2*(-5 + tw2) + s*(-3 + tw2)) + 
     &    mH2*(-5 + tw2))*C0(0,s,mH2,mW2,mW2,mW2)))/
     &   (mH2 - s)**2

    
      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function CFGHG(s)
      implicit none
      integer i
      real*8 s,nc
      real*8 mW2,mH2,mZ2,Qf(9),gv(9),mf(9),tw2,eps

      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps

*
*** Initialize LoopTools
*
#include "looptools.h"

      CFGHG=0d0
      do i=1,9
         if(i.le.3)then
            nc=1d0
         else
            nc=3d0
         endif

        CFGHG=CFGHG+nc*(-4*mf(i)**2*(-2*s*B0(mH2,mf(i)**2,mf(i)**2) + 
     &   2*s*B0(s,mf(i)**2,mf(i)**2) + 
     &   (mH2 - s)*(-2 + C0(0,s,mH2,mf(i)**2,mf(i)**2,mf(i)**2)*
     &       (mH2 - s - 4*mf(i)**2)))*Qf(i)**2)/(mH2 - s)**2

      enddo

      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function CFZHG(s)
      implicit none
      integer i
      real*8 s,cw2,nc
      real*8 mW2,mH2,mZ2,qf(9),gv(9),mf(9),tw2,eps

      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps
*
*** Initialize LoopTools
*
#include "looptools.h"

      cw2=1d0/(1d0+tw2)

      CFZHG=0d0
      do i=1,9
         if(i.le.3)then
            nc=1d0
         else
            nc=3d0
         endif

       CFZHG=CFZHG+nc*(4*gv(i)*mf(i)**2*
     & (-2*s*B0(mH2,mf(i)**2,mf(i)**2) + 2*s*B0(s,mf(i)**2,mf(i)**2) + 
     &   (mH2 - s)*(-2 + C0(0,s,mH2,mf(i)**2,mf(i)**2,mf(i)**2)*
     &       (mH2 - s - 4*mf(i)**2)))*Qf(i))/(cw2*(mH2 - s)**2)

       enddo

      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function c1GHGpaper(s)
      implicit none
      real*8 s
      real*8 mW2,mH2,mZ2,Qf(9),gv(9),mf(9),tw2,eps
      real*8 beta,betaH
      complex*16 J1,J2
      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps

      beta=s/mW2
      betaH=mH2/mW2

      c1GHGpaper= -16.*J1(beta,betaH) + 4.*(6d0+betaH)*J2(beta,betaH)

      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function d1ZHGpaper(s)
      implicit none
      real*8 s
      real*8 mW2,mH2,mZ2,Qf(9),gv(9),mf(9),tw2,eps
      real*8 beta,betaH
      complex*16 J1,J2
      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps

      beta=s/mW2
      betaH=mH2/mW2

      d1ZHGpaper=4d0*(3d0-tw2)*J1(beta,betaH)
     &     +4d0*(-5d0+tw2-0.5d0*betaH*(1d0-tw2))*J2(beta,betaH)

      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function CFGHGpaper(s)
      implicit none
      integer i
      real*8 s,nc
      real*8 mW2,mH2,mZ2,Qf(9),gv(9),mf(9),tw2,eps
      real*8 betap,betapH
      complex*16 J1,J2

      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps

*
*** Initialize LoopTools
*
#include "looptools.h"

      CFGHGpaper=0d0
      do i=6,6
         if(i.le.3)then
            nc=1d0
         else
            nc=3d0
         endif

         betap=s/mf(i)**2
         betapH=mH2/mf(i)**2

       CFGHGpaper=CFGHGpaper+nc*
     &        qf(i)**2*4d0*(J1(betap,betapH)-4d0*J2(betap,betapH))

      enddo

      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function CFZHGpaper(s)
      implicit none
      integer i
      real*8 s,cw2,nc
      real*8 mW2,mH2,mZ2,qf(9),gv(9),mf(9),tw2,eps
      real*8 betap,betapH
      complex*16 J1,J2

      common/data/mW2,mH2,mZ2,Qf,gv,mf,tw2,eps
*
*** Initialize LoopTools
*
#include "looptools.h"


      CFZHGpaper=0d0
      do i=6,6
         if(i.le.3)then
            nc=1d0
         else
            nc=3d0
         endif
         
         betap=s/mf(i)**2
         betapH=mH2/mf(i)**2

       CFZHGpaper=CFZHGpaper+nc*(1d0+tw2)*(-1d0)*4d0*
     &        qf(i)*gv(i)*(J1(betap,betapH)-4d0*J2(betap,betapH))

       enddo

      return
      end

*-------------------------------------------------------------------*
*-------------------------------------------------------------------*
      complex*16 function J1(x,y)
      implicit none
      real*8 x,y
      complex*16 S

      J1=-(S(x)-S(y))/(x-y)
    
      return
      end

*-------------------------------------------------------------------*
      complex*16 function J2(x,y)
      implicit none
      real*8 x,y
      complex*16 S,L

      J2=0.5d0/(x-y)+1d0/(x-y)**2*(S(x)-S(y)-0.5d0*x*(L(x)-L(y)))
    
      return
      end

*-------------------------------------------------------------------*
      complex*16 function S(x)
      implicit none
      real*8 x,pi,ReS,ImS

      pi=acos(-1d0)

      if(x.lt.4d0)then
         S=-2d0*(asin(sqrt(x)/2d0))**2
      else
         reS=-pi**2/2d0+0.5d0*(log((x-2d0+sqrt(x**2-4d0*x))/2d0))**2
         ImS=-2d0*pi*log(sqrt(x)/2d0+sqrt(x/4d0-1d0))
         S=cmplx(reS,ImS)
      endif
      return
      end



*-------------------------------------------------------------------*
      complex*16 function L(x)
      implicit none
      real*8 x,ReL,ImL,pi
      
      pi=acos(-1d0)

      if(x.lt.4d0)then
         L=-2d0+2d0*sqrt((4d0-x)/x)*asin(sqrt(x)/2d0)
      else
         ReL=-2d0+sqrt((x-4d0)/x)*log((x-2d0+sqrt(x*(x-4d0)))/2d0)
         ImL=-pi*sqrt((x-4d0)/x)
         L=cmplx(ReL,ImL)
      endif
      return
      end



*-------------------------------------------------------------------*


