program epem
      implicit none
      integer ndim,nevents
      parameter (ndim=4,nevents=5000000)
      real * 8 xx(ndim),xgrid(0:50,ndim),xint,ymax(50,ndim),sigtot,error
      common/cminthvq/xx,xgrid,xint,ymax,sigtot,error
      real * 8 q0,xlam,ptmin2
      common/physpars/q0,xlam,ptmin2
      real * 8 costh1bar
      common/cbbar/costh1bar
      real * 8 vec1(4),vec2(4),vec3(4)
      common/kine/vec1,vec2,vec3
      integer k
      real * 8 thrust,res
      real * 8 random,btilde,dgauss,a
      external random,btilde,dgauss,a
      q0=100d0
      xlam=0.24d0
      ptmin2=(2*xlam)**2
      call bbinit
      do k=1,2
         costh1bar=1-2*random()
         call testsudakov
      enddo
      call mbook(1,'Thrust',0.005d0,0.5d0,1d0)
      call mbook(2,'Thrust',0.005d0,0.5d0,1d0)
      call mbook(3,'Thrust',0.005d0,0.5d0,1d0)
      do k=1,nevents
         call gen(btilde,ndim,xgrid,ymax,1,xx)
c costh1bar has been generated
         call genptxy
c full POWHEG event has been generated
         thrust=max(vec1(4),vec2(4),vec3(4))
         thrust=thrust*2/q0
         call mfill(1,thrust,(1-thrust)/nevents)
c all kinematics has been generated
      enddo
c compute thrust at fixed order
      do k=101,200
         res=dgauss(a,(k-1)*0.005d0,k*0.005d0,1d-4)
         call mfill(2,(k-1)*0.005d0+0.0025d0,res)
      enddo
      call mfinal(1)
      call mfinal(2)
      call multitop(1,3,1,1,'thrust','thrust','LIN')
      call multitop(2,3,1,1,'thrust','thrust','LIN')
      end

      function a(t)
c formula for thrust distribution at fixed order,
c De Rujula etal, Nucl. Phys. B138 (1978) 387.
      implicit none
      real * 8 a,t
      real * 8 pi
      parameter (pi=3.141592653589793d0)
      real * 8 q0,xlam,ptmin2
      common/physpars/q0,xlam,ptmin2
      real * 8 alfas
      external alfas
c a=1/sigma d sigma /d t * (1-t) at order alpha_s
      a=4d0/3*(2*(3*t**2-3*t+2)/t*log((2*t-1)/(1-t))-3*(3*t-2)*(2-t))
     #  *alfas(q0**2/4,xlam,5)/(2*pi)
      end

      subroutine bbinit
      implicit none
      integer ndim,ncall1,ncall2,itmx1,itmx2,ncall3
      parameter (ndim=4)
      real * 8 xx(ndim),xgrid(0:50,ndim),xint,ymax(50,ndim),sigtot,error
      common/cminthvq/xx,xgrid,xint,ymax,sigtot,error
      real * 8 estimn,errorn,estimp,errorp,y,xi,phi,rrr
      real * 8 pi,cf
      parameter (pi=3.141592653589793d0,cf=4d0/3d0)
      real * 8 q0,xlam,ptmin2
      common/physpars/q0,xlam,ptmin2
      integer ifold(4)
      common/cifold/ifold
      logical negflag
      common/cnegflag/negflag
      real * 8 unorm
      common/cunorm/unorm
      real * 8 costh1bar
      common/cbbar/costh1bar
      real * 8 alphas
      common/calphas/alphas
      integer j
      real * 8 btilde,upper,random,kt2,born
      external btilde,upper,random,kt2,born
      ncall1=10000
      itmx1=5
c set up the grid
      negflag=.false.
      call mint(btilde,ndim,ncall1,itmx1,0,
     #        xgrid,xint,ymax,sigtot,error)
      write(*,*) sigtot*3d0/8,error*3d0/8
c
c setup folding
c Notice that with no folding (all ifold=1)
c we get around 2 per mill negative contribution to
c the integral, for  fixed alphas=0.5; with the
c choice below it disappears.
      ifold(1)=1
      ifold(2)=2
      ifold(3)=2
      ifold(4)=2
      ncall2=10000
      itmx2=5
      negflag=.true.
      call mint(btilde,ndim,ncall2,itmx2,1,
     #        xgrid,xint,ymax,estimn,errorn)
      negflag=.false.
      call mint(btilde,ndim,ncall2,itmx2,1,
     #        xgrid,xint,ymax,estimp,errorp)
c should get 1+alfas/pi:
      write(*,*) (estimp+estimn)*3d0/8,' +- ',
     #  sqrt(errorp**2+errorn**2)*3d0/8
      write(*,*)  'Shoud equal ',1+alphas/pi
      negflag=.false.
c set up the normalization of the upper bounding function
      unorm=-1
      ncall3=10000
      do j=1,ncall3
         costh1bar=1-2*random()
         y=1-2*random()
         xi=random()
         if(kt2(xi,y).ge.ptmin2) then
            phi=2*pi*random()
            call rvals(xi,y,phi,rrr)
            unorm=max((rrr/born())/upper(xi,y),unorm)
         endif
      enddo
c increase unorm, just to be on the safe side ...
      unorm=unorm*2
c initialize gen
      call gen(btilde,ndim,xgrid,ymax,0,xx)
      end

      function born()
      implicit none
      real * 8 born
      real * 8 costh1bar
      common/cbbar/costh1bar
      born=1+costh1bar**2
      end

      function upper(xi,y)
      implicit none
      real * 8 pi
      parameter (pi=3.141592653589793d0)
      real * 8 upper,xi,y
      real * 8 q0,xlam,ptmin2
      common/physpars/q0,xlam,ptmin2
      real * 8 alfas,b0
      b0=(33-2*5)/(12*pi)
      alfas=1/(b0*log(q0**2/4*(1-y**2)*xi**2/xlam**2))
      upper=1/(1-y**2)/xi * alfas
      end


      function kt2(xi,y)
      implicit none
      real * 8 kt2,xi,y
      real * 8 q0,xlam,ptmin2
      common/physpars/q0,xlam,ptmin2
      kt2=(q0/2)**2*(1-y**2)*xi**2
      end

      subroutine kinematics(xi,y,phi,vec1,vec2,vec3)
      implicit none
      real * 8 xi,y,phi,vec1(4),vec2(4),vec3(4)
      real * 8 q0,xlam,ptmin2
      common/physpars/q0,xlam,ptmin2
      real * 8 costh1bar
      common/cbbar/costh1bar
      real * 8 barredk(4,2),xk(4,3),jac
      integer j
      barredk(1,1)=0
      barredk(2,1)=q0/2*sqrt(1-costh1bar**2)
      barredk(3,1)=q0/2*costh1bar
      barredk(4,1)=q0/2
      barredk(4,2)=q0/2
      do j=1,3
         barredk(j,2)=-barredk(j,1)
      enddo
      call barradmap(2,q0,barredk,xi,y,phi,xk,jac)
      do j=1,4
         vec1(j)=xk(j,1)
         vec2(j)=xk(j,2)
         vec3(j)=xk(j,3)
      enddo
      end

      function btilde(xxx,weight,ifl)
      implicit none
      real * 8 btilde,xxx(4),weight
      integer ifl
      real * 8 pi,cf
      parameter (pi=3.141592653589793d0,cf=4d0/3d0)
      real * 8 d31,d32
      real * 8 q0,xlam,ptmin2
      common/physpars/q0,xlam,ptmin2
      real * 8 a,b,q
c q is the Ellis-Sexton scale; is set to arbitrary values
      parameter (a=2d0,b=1d0,q=1)
      real * 8 vec1(4),vec2(4),vec3(4),y,xi,phi,x1,x2,x3,cth1,cth2
      real * 8 bbb,svirt,vfin,qcal,tcal12,accum,xjac0,xjac1
      real * 8 Axiy,Axi1,A01,A0y,r31,r32,r12,s32
      logical negflag
      common/cnegflag/negflag
      real * 8 alphas
      common/calphas/alphas
c we will need to leave in a common block the generated value of costh1bar
c after each call to gen
      real * 8 costh1bar
      common/cbbar/costh1bar
      save bbb,svirt,accum
      real * 8 alfas,born
      external alfas,born
c fix flavours to 5 in this example
      alphas=alfas(q0**2,xlam,5)
      costh1bar=1-2*xxx(1)
      xjac0=2
      xjac1=xjac0
c Do importance sampling of y around y=-1
      y=-1+2*xxx(2)**2
      xjac1=xjac1*4*xxx(2)
c Do importance sampling of xi around xi=1
      xi=1-xxx(3)**2
      xjac1=xjac1*2*xxx(3)
c
      phi=2*pi*xxx(4)
      xjac1=xjac1*2*pi
c
      if(ifl.eq.2) then
         btilde=accum
         accum=0
         if(negflag) then
            if(btilde.gt.0) btilde=0
         else
            if(btilde.lt.0) btilde=0
         endif
         return
      endif
      if(ifl.eq.0) then
         bbb=born()
         qcal=2*cf*((13d0/2-2*pi**2/3)-3d0/2*log(q0**2/q**2))
         tcal12=1d0/2*log(q0**2/q**2)**2-pi**2/6
         vfin = cf * (pi**2-8+3*log(q0**2/Q**2) - log(q0**2/Q**2)**2)
         svirt=alphas/(2*pi)*(qcal+2*cf*tcal12+vfin)*bbb
         accum=0
      endif
c R32 contribution
      call kinematics(xi,y,phi,vec1,vec2,vec3)
      x1=2*vec1(4)/q0
      x2=2*vec2(4)/q0
      x3=2*vec3(4)/q0
      d32=(x3*x2)**(a-b)*(1-x1)**b
      d31=(x3*x1)**(a-b)*(1-x2)**b
      cth1=vec1(3)/vec1(4)
      cth2=vec2(3)/vec2(4)
      s32=(1/d32)/(1/d32+1/d31)
      Axiy=(x1**2*(1+cth1**2)+x2**2*(1+cth2**2))
     #*(2-xi*(1-y))**2/((1+y)*(1-xi))
     # * s32
c collinear kinematics, y=1
      cth1=costh1bar
      cth2=-costh1bar
      s32=1
      Axi1=((1+cth1**2)+(1-xi)**2*(1+cth2**2))*
     #     4/(2*(1-xi))*s32
c soft collinear kinematics,y=1 + xi=0      
      A01=((1+cth1**2)+(1+cth2**2))*2 * s32
c soft kinematics, xi->0, y#1
      s32=(1+y)**b/((1-y)**b+(1+y)**b)
      A0y=((1+cth1**2)+(1+cth2**2))
     #     *4/(1+y)*s32
c full contribution
      r32=(8*pi*alphas*cf)/q0**2* q0**2/(4*pi)**3/(xi*(1-y))*(
     # (x2**2/(1-xi)*Axiy-(1-xi)*Axi1)
     #-(A0y-A01))
c Should repeat the same for r31; just double r32
      btilde = ((bbb + svirt)*xjac0 + 2*r32 * xjac1)*weight
      accum=accum+btilde
      end

      subroutine rvals(xi,y,phi,r32)
      implicit none
      real * 8 btilde,xi,y,phi
      real * 8 pi,cf
      parameter (pi=3.141592653589793d0,cf=4d0/3d0)
      real * 8 d31,d32
      real * 8 a,b,alphas,q
      parameter (a=2d0,b=1d0,q=2)
      real * 8 vec1(4),vec2(4),vec3(4),x1,x2,x3,cth1,cth2
      real * 8 born,svirt,vfin,qcal,tcal12
      real * 8 Axiy,r31,r32,s32
      real * 8 q0,xlam,ptmin2
      common/physpars/q0,xlam,ptmin2
      real * 8 costh1bar
      common/cbbar/costh1bar
      real * 8 alfas,kt2
      external alfas,kt2
      alphas=alfas(kt2(xi,y),xlam,5)
c R32 contribution
      call kinematics(xi,y,phi,vec1,vec2,vec3)
      x1=2*vec1(4)/q0
      x2=2*vec2(4)/q0
      x3=2*vec3(4)/q0
      d32=(x3*x2)**(a-b)*(1-x1)**b
      d31=(x3*x1)**(a-b)*(1-x2)**b
      cth1=vec1(3)/vec1(4)
      cth2=vec2(3)/vec2(4)
      s32=(1/d32)/(1/d32+1/d31)
      Axiy=(x1**2*(1+cth1**2)+x2**2*(1+cth2**2))
     #*(2-xi*(1-y))**2/((1+y)*(1-xi))
     # * s32
c full contribution
      r32=(8*pi*alphas*cf)/q0**2* q0**2/(4*pi)**3/(xi*(1-y))*(
     # (x2**2/(1-xi)*Axiy))
c double it to account for the 31 region
      r32=2*r32
      end



      function pt2solve(pt2,i)
c Returns  xlr - log(\Delta^{(\tilde{V})}) , see eq. D14, D15 in ZZ paper
c We use it to find its zero in pt2.
      implicit none
      real * 8 pt2solve,pt2
c i set by dzero: 1 for first call, 2 for subsequent calls, 3 for last call
c before a normal exit; not used here
      integer i
      real * 8 pi
      parameter (pi=3.141592653589793d0)
      real * 8 xlr,xlam2c,ktmax2,cunorm
      common/cpt2solve/xlr,ktmax2,xlam2c,cunorm
      real * 8 b0
      b0=(33-2*5)/(12*pi)
      pt2solve=cunorm*pi/b0
     #    *(log(4*ktmax2/xlam2c)*log(log(ktmax2/xlam2c)/log(pt2/xlam2c))
     #     - log(ktmax2/pt2)) + xlr
      end


      subroutine genptxy
      implicit none
      real * 8 pi
      parameter (pi=3.141592653589793d0)
c generate radiation kinematics
      real * 8 q0,xlam,ptmin2
      common/physpars/q0,xlam,ptmin2
      real * 8 unorm
      common/cunorm/unorm
      real * 8 xlr,xlam2c,kmax2,cunorm
      common/cpt2solve/xlr,kmax2,xlam2c,cunorm
      real * 8 pt2,rprime,rsecond,rv,lyfr,xi,y,phi,rrr,err
      common/genpt2/pt2
      real * 8 costh1bar
      common/cbbar/costh1bar
      real * 8 random,upper,pt2solve,born
      external random,upper,pt2solve,born
c the following is to prevent error messages
c for out of range in dzero calls
      call KERSET('C205. ',0,0,101)
c
      kmax2=q0**2/4
      xlam2c=xlam**2
      cunorm=unorm
      xlr=0
 1    continue
c this is log(r), r being the random number
      xlr=xlr+log(random())
c dzero(xmin,xmax,x(return value),err(return value),requested accuracy,
c max number of calls, f(x))
      call dzero(ptmin2,kmax2,pt2,err,1d-8,1000000,pt2solve)
c error conditions
      if(pt2.eq.0.and.err.lt.0d0
     # .and.err.gt.ptmin2-kmax2) then
         write(*,*) 'DZERO fails'
         write(*,*) ' number of calls exceeded'
         stop
      endif
      if(pt2.lt.ptmin2) then
c generate a 2-body event
         call twobody
         return
      endif
      rprime=random()
c check the upper bounds for safety!
      if(log(4*kmax2/pt2)
     #     .lt.log((1+sqrt(1-pt2/kmax2))/(1-sqrt(1-pt2/kmax2)))) then
         write(*,*) ' upper bound!'
         stop
      endif
      if(rprime*log(4*kmax2/pt2).gt.
     # log((1+sqrt(1-pt2/kmax2))/(1-sqrt(1-pt2/kmax2))) ) goto 1
      rv=random()
      lyfr=(1-2*rv)
     #   *log((1+sqrt(1-pt2/kmax2))/(1-sqrt(1-pt2/kmax2)))
      y=(exp(lyfr)-1)/(exp(lyfr)+1)
      xi=2*sqrt(pt2)/(q0*sqrt(1-y**2))
c last veto
      rsecond=random()
      phi=2*pi*random()
      call rvals(xi,y,phi,rrr)
      rrr=rrr/born()
c check the upper bound
      if(unorm*upper(xi,y).lt.rrr) then
         write(*,*) ' upper bound!'
         stop
      endif
      if(rsecond*unorm*upper(xi,y).gt.rrr) goto 1
c veto passed; generate the event
      call threebody(xi,y,phi)
      end

      subroutine increasecnt
      end

      subroutine twobody
      implicit none
      real * 8 q0,xlam,ptmin2
      common/physpars/q0,xlam,ptmin2
      real * 8 vec1(4),vec2(4),vec3(4)
      common/kine/vec1,vec2,vec3
      real * 8 costh1bar
      common/cbbar/costh1bar
      integer j
      vec1(1)=0
      vec1(2)=q0/2*sqrt(1-costh1bar**2)
      vec1(3)=q0/2*costh1bar
      vec1(4)=q0/2
      vec2(4)=q0/2
      do j=1,3
         vec2(j)=-vec1(j)
         vec3(j)=0
      enddo
      vec3(4)=0
      end

      subroutine threebody(xi,y,phi)
      implicit none
      real * 8 xi,y,phi
      real * 8 vec1(4),vec2(4),vec3(4)
      common/kine/vec1,vec2,vec3
      call kinematics(xi,y,phi,vec1,vec2,vec3)
      end

      subroutine testsudakov
c at fixed costh1bar, computes the sudakov form factor
c by computing directly int R/B \theta(kt-pt), and by
c calling genptxy
      implicit none
      integer ncalls,nbins
      parameter (ncalls=50000,nbins=40)
      real * 8 pi
      parameter (pi=3.141592653589793d0)
      real * 8 q0,xlam,ptmin2
      common/physpars/q0,xlam,ptmin2
      real * 8 costh1bar
      common/cbbar/costh1bar
      real * 8 pt2
      common/genpt2/pt2
      real * 8 hist1(nbins),hist2(nbins)
      real * 8 xi,y,phi,kmax2,rrr
      integer j,k
      real * 8 random
      external random
      do j=1,nbins
         hist1(j)=0
         hist2(j)=0
      enddo
      kmax2=(q0/2)**2
c First compute using direct integration
      do k=1,ncalls
         xi=random()
         y=1-2*random()
         phi=2*pi*random()
         pt2=(q0/2)**2*(1-y**2)*xi**2
         if(pt2.gt.ptmin2) then
            call rvals(xi,y,phi,rrr)
c supply phase space volume
            rrr=4*pi*rrr/(1+costh1bar**2)
            do j=1,nbins
               if(log(pt2/ptmin2).gt.j*log(kmax2/ptmin2)/nbins)
     #              hist1(j)=hist1(j)+rrr
            enddo
         endif
      enddo
c Now compute it using genptxy
      do k=1,ncalls
         call genptxy
         do j=1,nbins
            if(pt2.gt.0.and.
     #           log(pt2/ptmin2).gt.j*log(kmax2/ptmin2)/nbins)
     #           hist2(j)=hist2(j)+1d0
         enddo
      enddo
      write(42,*) ' newplot'
      write(42,*) ' set scale x log'
      do k=1,nbins
         write(42,*) exp(k*log(kmax2/ptmin2)/nbins)*ptmin2,
     #                   exp(-hist1(k)/ncalls)
      enddo
      write(42,*) ' join'
      do k=1,nbins
         write(42,*) exp(k*log(kmax2/ptmin2)/nbins)*ptmin2,
     #                   1-hist2(k)/ncalls
      enddo
      write(42,*) ' join dashes'
      end