LHAPDF/lhapdf5.3.1/wrapEVLCTEQ.f

   1       subroutine EVLCTEQevolve(x,Q,f)
   2       implicit none
   3       include 'parmsetup.inc'
   4       character*16 name(nmxset)
   5       integer nmem(nmxset),ndef(nmxset),mmem
   6       common/NAME/name,nmem,ndef,mmem
   7       real*8 f(-6:6)
   8       real*8 x,Q,xx,qq,up,dn,ubar,dbar,gluon,sbar,cbar,bbar,tbar,Blam
   9       real*8 xmin,qmax,anx,anq,Q2fit,qini,Ahdn,ordi,anf,mass,pi,alfas
  10       real*8 Q0,alfas0
  11       real*8 CtLhpardis,CtLhALPI,parmflavor,CtLhastolam
  12       integer nflavor,Eorder,j,nx,nq
  13       integer nset
  14       data pi / 3.141592653589793d0 /
  15       external parmflavor
  16       save nflavor,xmin,qmax,anx,anq
  17
  18 *
  19         xx = x
  20         qq = Q
  21         up   = xx*CtLhpardis(1,xx,qq)
  22         dn   = xx*CtLhpardis(2,xx,qq)
  23         ubar = xx*CtLhpardis(-1,xx,qq)
  24         dbar = xx*CtLhpardis(-2,xx,qq)
  25         gluon= xx*CtLhpardis( 0,xx,qq)
  26
  27         if(nflavor .ge. 3) then
  28            sbar = xx*CtLhpardis(-3,xx,qq)
  29         else
  30            sbar = 0.d0
  31         endif
  32
  33         if(nflavor .ge. 4) then
  34            cbar = xx*CtLhpardis(-4,xx,qq)
  35         else
  36            cbar = 0.d0
  37         endif
  38
  39         if(nflavor .ge. 5) then
  40            bbar = xx*CtLhpardis(-5,xx,qq)
  41         else
  42            bbar = 0.d0
  43         endif
  44
  45         if(nflavor .eq. 6) then
  46            tbar = xx*CtLhpardis(-6,xx,qq)
  47         else
  48            tbar = 0.d0
  49         endif
  50
  51         f(6) = tbar
  52         f(5) = bbar
  53         f(4) = cbar
  54         f(3) = sbar
  55         f(2) = up
  56         f(1) = dn
  57         f(0) = gluon
  58         f(-1) = dbar
  59         f(-2) = ubar
  60         f(-3) = sbar
  61         f(-4) = cbar
  62         f(-5) = bbar
  63         f(-6) = tbar
  64 *
  65         return
  66 *
  67       entry EVLCTEQread(nset)
  68       call CtLhbldat1
  69       call CtLhbldat2
  70       read(1,*) xmin,qmax,nx,nq
  71       anx=nx
  72       anq=nq-1
  73       return
  74 *
  75       entry EVLCTEQalfa(alfas,Q)
  76       alfas = pi*CtLhALPI(Q)
  77       return
  78 *
  79       entry EVLCTEQinit(nset,Eorder,Q2fit)
  80 *
  81       Ahdn=1d0
  82       call CtLhParPdf(1,'IHDN',Ahdn,j)
  83       call CtLhParPdf(1,'QMAX',qmax,j)
  84       call CtLhParPdf(1,'XMIN',xmin,j)
  85       call GetOrderAsM(nset,j)
  86       ordi=j+1.0
  87       call CtLhParQcd(1,'ORDR',ordi,j)
  88       call GetThresholdM(nset,4,mass)
  89       call CtLhParQcd(1,'M4',mass,j)
  90       call GetThresholdM(nset,5,mass)
  91       call CtLhParQcd(1,'M5',mass,j)
  92       mass=180d0
  93       call CtLhParQcd(1,'M6',mass,j)
  94       call GetNfM(nset,j)
  95       nflavor=j
  96       Anf=nflavor
  97       call CtLhParQcd(1,'NFL',Anf,j)
  98       qini=sqrt(Q2fit)
  99       call CtLhParPdf(1,'QINI',qini,j)
 100       call CtLhEvlpar(1,'NFMX', Anf, j)
 101       ordi=Eorder+1.0
 102       call CtLhParPdf(1,'IKNL',ordi,j)
 103       return
 104 *
 105       entry EVLCTEQpdf(nset)
 106       call GetOrderAsM(nset,j)
 107       call GetAlfas(nset,alfas0,Q0)
 108       Blam=CtLhastolam(alfas0,Q0,j+1,nflavor)
 109       call CtLhSetLam (nflavor,Blam,j+1)
 110       call CtLhParPdf(1,'NX',anx,j)
 111       call CtLhParPdf(1,'NT',anq,j)
 112       call CtLhEvolve (parmflavor,j)
 113       if (j .ne. 0) then
 114          write(*,*) 'EVLCTEQ Evolve Error code :',j
 115          stop
 116       endif
 117       return
 118 *
 119       end
 120 *
 121       double precision function parmflavor(i,x)
 122       implicit none
 123       real*8 x,f(-6:6)
 124       integer i,i0
 125       integer iset
 126 *
 127       call getnset(iset)
 128       call parmPDF(iset,x,f)
 129       i0=i
 130       if (i.eq.-2) i0=-1
 131       if (i.eq.-1) i0=-2
 132       if (i.eq.1) i0=2
 133       if (i.eq.2) i0=1
 134       parmflavor=f(i0)/x
 135       return
 136 *
 137       end
 138
 139       function CtLhastolam(as,q,nloop,nf)
 140       implicit double precision (a-h,o-z)
 141       data pi / 3.141592653589793d0 /
 142       xlp = nloop-1d0
 143       b  = (33.0-2.0*nf)/pi/12.0
 144       bp = (153.0 - 19.0*nf) / pi / 2.0 / (33.0 - 2.0*nf) * xlp
 145       t  = 1.0/b/as
 146 c-----------------------------------------------------------
 147 c Solve the equation
 148 c
 149     1 xlt = log(t)
 150       ot = t
 151 c-----------------------------------------------------------
 152 c Solve the equation
 153 c Value and Derivative of alfa with respect to t
 154 c
 155       as0  = 1/b/t - bp*xlt/(b*t)**2
 156       as1  = - 1/b/t**2 -bp/b**2*(1-2*xlt)/t**3
 157       t  = (as-as0)/as1 + t
 158       if(abs(ot-t)/ot.gt..00001)goto 1
 159       CtLhastolam = q/exp(t/2.0)
 160       return
 161       end
 162