C     ----------------------------------------------------------------
      program convolit
C     ----------------------------------------------------------------
C--   Test convolution engine against Gauss integration 
C--   VFNS w/backward NNLO evolution on multiple grids 
C     ----------------------------------------------------------------

C--   Michiel Botje   h24@nikhef.nl

      implicit double precision (a-h,o-z)
      
      data iord/2/, nfin/0/, iosp/3/               !NNLO, VFNS, iospl
      data as0/0.364/, r20/2.D0/                   !alphas 
      data q2c/7.D0/, q2b/25.D0/, q2t/1000.D0/     !thresholds
      data q20/5.D0/                               !q20

      dimension xmin(5),ixwt(5)
      data xmin/1.D-5,0.2D0,0.4D0,0.6D0,0.75D0/    !multiple x grid
      data ixwt/    1,    2,    4,    8,    16/    !xgrid weights
      data nxin/100/                               !x grid pnts

      dimension qq(2),qw(2)
      data qq/2.D0,1.D4/, qw/1.D0,1.D0/, nqin/60/  !mu2 grid

      external func                                !input parton dists
      dimension def(-6:6,12)                       !flavor composition
      data def  /
C--   tb  bb  cb  sb  ub  db   g   d   u   s   c   b   t
C--   -6  -5  -4  -3  -2  -1   0   1   2   3   4   5   6  
     + 0., 0., 0., 0., 0.,-1., 0., 1., 0., 0., 0., 0., 0.,   !dval
     + 0., 0., 0., 0.,-1., 0., 0., 0., 1., 0., 0., 0., 0.,   !uval
     + 0., 0., 0.,-1., 0., 0., 0., 0., 0., 1., 0., 0., 0.,   !sval
     + 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,   !dbar
     + 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,   !ubar
     + 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,   !sbar
     + 0., 0.,-1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,   !cval
     + 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,   !cbar
     + 52*0.    /
     
      real time1, time2
      
      dimension itypes(4)
      data itypes /0,0,0,0/
      dimension wx(9000)

      dimension idw(4)
      data idw   /0,0,0,0/
      
      dimension gxp(3),dbl(3),gxg(3)
            
      external cvol,pvol,dfun,qpdvol,afun,achi
      common /pass/ x, qmu2, iset, ipdf
      
C--   lun = 6 stdout, -6 stdout w/out banner page
      lun    = 6
      lunout = abs(lun) 
      
      call qcinit(lun,' ')
      call gxmake(xmin,ixwt,5,nxin,nx,iosp)
      call gqmake(qq,qw,2,nqin,nq)

C--   Try to read the weight file and create one if that fails      
      call readwt(20,'../weights/unpol.wgt',id1,id2,nw,ierr)
      if(ierr.ne.0) then
        call fillwt(1,id1,id2,nw)
        call dmpwgt(1,20,'../weights/unpol.wgt')
      endif
  
C--   Evolve unpolarised        
      call setord(iord)                      !LO, NLO, NNLO
      call setalf(as0,r20)                   !input alphas
      iqc  = iqfrmq(q2c)                     !charm threshold
      iqb  = iqfrmq(q2b)                     !bottom threshold
      iqt  = iqfrmq(q2t)                     !top threshold
      call setcbt(nfin,iqc,iqb,iqt)          !thesholds in the vfns
      iq0  = iqfrmq(q20)                     !starting scale
      call cpu_time(time1)                   !start timing
      call evolfg(1,func,def,iq0,eps)        !evolve all pdf's
      call cpu_time(time2)
      
C--   Book and fill weight tables
      itypes(1) = 2                          !two type-1 tables (x)
      itypes(2) = 2                          !two type-2 tables (x,nf)
      call BookTab(wx,9000,itypes,nwrds)     !book tables
      call MakeWtX(wx,101)                   !weights for f * f
      call MakeWtA(wx,102,afun,achi)         !convolution kernel K
      idPGQ = idspfun('PGQ',1,1)             !Pgq in QCDNUM memory
      call CopyWgt(wx,idPGQ,201,0)           !copy Pgq to weight store 
      call WcrossW(wx,201,102,202,0)         !weights for Pgq * K 

      iset  = 1                              !unpolarised 
      ipdf  = 0                              !gluon
      ix    = 1
      iq    = iq0
      x     = xfrmix(ix)                     !interpolation point
      qmu2  = qfrmiq(iq)                     !interpolation point
      
C--   Gauss integration 
      gxp(1) = dmb_gauss(pvol,x,1.D0,1.D-5)           !Pgq * gluon
      dbl(1) = dmb_gauss(qpdvol,x,1.D0,1.D-5)         !K * Pgq * gluon
      gxg(1) = dmb_gauss(cvol,x,1.D0,1.D-5)           !gluon * gluon

C--   Slow convolution engine
      gxp(2) = FcrossK(wx,201,iset,ipdf,ix,iq)        !Pgq * gluon
      dbl(2) = FcrossK(wx,202,iset,ipdf,ix,iq)        !K * Pgq * gluon
      gxg(2) = FcrossF(wx,101,iset,ipdf,ipdf,ix,iq)   !gluon * gluon
      
C--   Fast convolution engine 
      call FastIni(x,qmu2,1,1)               !pass x,q to fast engine
      call FastEpm(iset,ipdf,1)              !1 = gluon fast buffer 
          
      idw(1) = 201                           !Pgq indentifier
      call FastFxK(wx,idw,1,2)               !2 = Pgq * gluon
      call FastFxq(2,gxp(3),1)               !Interpolate
      
      idw(1) = 201                           !Pgq weight table id
      call FastFxK(wx,idw,1,-2)              !2 = Pgq * gluon
      idw(1) = 102                           !K weight table id 
      call FastFxK(wx,idw,2,3)               !3 = K * Pgq * gluon
      call FastFxq(3,dbl(3),1)               !Interpolate
      
      call FastFxF(wx,101,1,1,2)             !2 = gluon * gluon
      call FastFxq(2,gxg(3),1)               !Interpolate

C--   Compare
      write(lunout,'('' gauss, qcdnum, fast = '',3F14.5)') gxp
      write(lunout,'('' gauss, qcdnum, fast = '',3F14.5)') dbl
      write(lunout,'('' gauss, qcdnum, fast = '',3F14.5)') gxg
      
      write(lunout,'(/''time spent:'',F7.4)') time2-time1         

      end
      
C     ----------------------------------------------------------------

C     =================================
      double precision function cvol(z)
C     =================================

C--   Gluon convolution

      implicit double precision (a-h,o-z)
      
      common /pass/ x, qmu2, iset, ipdf
      
      r    = x/z
      cvol = fsnsxq(iset,ipdf,z,qmu2,1) * 
     +       fsnsxq(iset,ipdf,r,qmu2,1) / (z)
      
      return
      end
      
C     ----------------------------------------------------------------

C     =================================
      double precision function pvol(z)
C     =================================

C--   Gluon convolution with PGQ

      implicit double precision (a-h,o-z)
      
      common /pass/ x, qmu2, iset, ipdf
      
      r    = x/z
      pvol = fsnsxq(iset,ipdf,z,qmu2,1) * r * dqcP0gfa(r,3) / (z)
      
      return
      end
 
C     ----------------------------------------------------------------

C     ===================================
      double precision function qpdvol(z)
C     ===================================

C--   Gluon convolution with PGQ and dfun

      implicit double precision (a-h,o-z)
      
      common /pass/ x, qmu2, iset, ipdf
      
      r      = x/z
      qpdvol = fsnsxq(iset,ipdf,z,qmu2,1) * pdconv(r) / (z)
      
      return
      end

           
C     ===================================      
      double precision function pdconv(x)
C     ===================================

      implicit double precision (a-h,o-z)
      
      external pdfun
      
      common /pass2/ y 
      
      y      =  x
      pdconv =  dmb_gauss(pdfun,x,1.D0,1.D-5)
      
      return
      end
     
C     ==================================       
      double precision function pdfun(z)
C     ==================================

C--   Dfun convolution with PGQ

      implicit double precision (a-h,o-z)
      
      common /pass2/ y
      
      r     = y/z
      pdfun = z * dfun(z) * r * dqcP0gfa(r,3) /z
      
      return
      end
                  
      
C     =================================      
      double precision function dfun(x)
C     =================================

C--   Mickey Mouse kernel

      implicit double precision (a-h,o-z)
      
      dfun = x*(1.D0-x)
      
      return
      end
      
C     ----------------------------------------------------------------

C     ====================================
      double precision function achi(qmu2)      
C     ====================================

      implicit double precision (a-h,o-z)
      
      ddum = qmu2
      achi = 1.D0
      
      return
      end                        
     
C     =========================================     
      double precision function afun(x,qmu2,nf)
C     =========================================

      implicit double precision (a-h,o-z)
      
      ddum = qmu2
      idum = nf
      afun = dfun(x)
      
      return
      end                        
      
C     ----------------------------------------------------------------

C     ======================================
      double precision function func(ipdf,x)
C     ======================================

      implicit double precision (a-h,o-z)

                     func = 0.D0
      if(ipdf.eq. 0) func = xglu(x)
      if(ipdf.eq. 1) func = xdnv(x)
      if(ipdf.eq. 2) func = xupv(x)
      if(ipdf.eq. 3) func = 0.D0
      if(ipdf.eq. 4) func = xdbar(x)
      if(ipdf.eq. 5) func = xubar(x)
      if(ipdf.eq. 6) func = xsbar(x)
      if(ipdf.eq. 7) func = 0.D0
      if(ipdf.eq. 8) func = 0.D0
      if(ipdf.eq. 9) func = 0.D0
      if(ipdf.eq.10) func = 0.D0
      if(ipdf.eq.11) func = 0.D0
      if(ipdf.eq.12) func = 0.D0

      return
      end
 
C     =================================
      double precision function xupv(x)
C     =================================

      implicit double precision (a-h,o-z)

      data au /5.107200D0/

      xupv = au * x**0.8D0 * (1.D0-x)**3.D0

      return
      end

C     =================================
      double precision function xdnv(x)
C     =================================

      implicit double precision (a-h,o-z)

      data ad /3.064320D0/

      xdnv = ad * x**0.8D0 * (1.D0-x)**4.D0

      return
      end

C     =================================
      double precision function xglu(x)
C     =================================

      implicit double precision (a-h,o-z)

      data ag    /1.7D0 /

      common /msum/ glsum, uvsum, dvsum

      xglu = ag * x**(-0.1D0) * (1.D0-x)**5.D0

      return
      end

C     ==================================
      double precision function xdbar(x)
C     ==================================

      implicit double precision (a-h,o-z)

      data adbar /0.1939875D0/

      xdbar = adbar * x**(-0.1D0) * (1.D0-x)**6.D0

      return
      end

C     ==================================
      double precision function xubar(x)
C     ==================================

      implicit double precision (a-h,o-z)

      xubar = xdbar(x) * (1.D0-x)

      return
      end

C     ==================================
      double precision function xsbar(x)
C     ==================================

      implicit double precision (a-h,o-z)

      xsbar = 0.2D0 * (xdbar(x)+xubar(x))

      return
      end







