LHAPDF/lhapdf5.2.2/wrapQCDNUM4.f

   1       subroutine QCDNUM4evolve(x,Q,pdf)
   2       implicit none
   3       character*64 gridname
   4       character*16 s1,dummy
   5       integer Eorder,index
   6       real*8 x,Q,Qlam,Q2,Q2fit,alfas0,scale0,alfas
   7       real*8 mc,mc2,mb,mb2,mt,mt2,tc,tc2,tb,tb2
   8       real*8 singlet,dm,um,dp,up,sp,ub,db,sb,cb,bb,bp,cp,gl,xx
   9       real*8 QALFAS,QPDFXQ,XFROMIX
  10       real*8 f(-6:6),pdf(-6:6)
  11       integer iq0,iqc,iqb,nf,nf0,qnerr,NFLGET,iflag,ix,i,IQFROMQ
  12       integer nx,nq
  13       integer nset,iset,isetlast
  14       data isetlast/-1/
  15       real*8 xmin,xmax,qmin,qmax,S
  16       save iq0,iqc,iqb,nf0,mc2,mb2,tc2,tb2,mt2
  17       save nx,xmin,xmax,nq,qmin,qmax,gridname
  18 *
  19 c      print *,'QCDNUM4evolve'
  20       call getnset(iset)
  21       if (iset.ne.isetlast) then
  22         call get_pdfqcd(iset)
  23         isetlast = iset
  24       endif
  25 c
  26       Q2=Q*Q
  27       nf=4
  28       if (Q2.lt.tc2) nf=3
  29       pdf(0)=QPDFXQ('GLUON'  ,x,Q2,iflag)
  30       singlet=QPDFXQ('SINGLET',x,Q2,iflag)
  31       dm= QPDFXQ('DM',x,Q2,IFLAG)
  32       um= QPDFXQ('UM',x,Q2,IFLAG)
  33       dp= QPDFXQ('DP',x,Q2,IFLAG)
  34       up= QPDFXQ('UP',x,Q2,IFLAG)
  35       sp= QPDFXQ('SP',x,Q2,IFLAG)
  36       ub=0.5d0*(up-um+singlet/dble(nf))
  37       db=0.5d0*(dp-dm+singlet/dble(nf))
  38       sb=0.5d0*(sp+singlet/dble(nf))
  39       cb=0.d0
  40       if (nf.ge.4) then
  41          cp= QPDFXQ('CP',X,Q2,IFLAG)
  42          cb=0.5d0*(cp+singlet/dble(nf))
  43       end if
  44
  45       pdf(1)=dm+db
  46       pdf(2)=um+ub
  47       pdf(3)=sb
  48       pdf(4)=cb
  49       pdf(5)=0d0
  50       pdf(6)=0d0
  51       pdf(-1)=db
  52       pdf(-2)=ub
  53       pdf(-3)=sb
  54       pdf(-4)=cb
  55       pdf(-5)=0d0
  56       pdf(-6)=0d0
  57       return
  58 *
  59       entry QCDNUM4alfa(alfas,Q)
  60       Q2=Q*Q
  61       nf=4
  62       if (Q2.lt.mc2) nf=3
  63       alfas=QALFAS(Q2,Qlam,nf,iflag)
  64       return
  65 *
  66       entry QCDNUM4read(nset)
  67 c      read(1,*) s1
  68 c      print *,s1
  69 c         gridname='large.grid'
  70 c         nx=400
  71 c         xmin=1d-6
  72 c         xmax=1d0
  73 c         nq=112
  74 c         qmin=1d0
  75 c         qmax=1d10
  76 c         read(1,*) dummy
  77 *      else
  78          read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
  79 c*      endif
  80 c      iset = nset
  81       return
  82 *
  83       entry QCDNUM4init(nset,Eorder,Q2fit)
  84       call qninit
  85       call QNISET('ORDER',Eorder+1)
  86       call grxdef(nx,xmin)
  87       call grqdef(nq,qmin,qmax)
  88       call grqinp(Q2fit,1)
  89       call getQMassM(nset,4,mc)
  90       mc2=mc*mc
  91 c      print *,mc
  92       CALL QNRSET('CMASS',mc)
  93       CALL QNRSET('MCALF',mc)
  94       call getThresholdM(nset,4,tc)
  95       tc2=tc*tc
  96       call GRQINP(tc2,1)
  97       CALL GRQINP(tc2-1.0d-4,1)
  98 c      print *,'setting iq0',Q2fit
  99       call qthres(tc2,2d10)
 100       iq0=IQFROMQ(Q2fit)
 101       iqc=IQFROMQ(tc2)
 102       nf0= NFLGET(iq0)
 103 c      print *,iq0,iqc,nf0,q2fit,tc2
 104 c      nf0=3
 105 c      print *,nf0
 106       CALL QNBOOK(2,'dm')
 107       CALL QNBOOK(3,'um')
 108       CALL QNBOOK(4,'dp')
 109       CALL QNBOOK(5,'up')
 110       CALL QNBOOK(6,'sp')
 111       CALL QNBOOK(7,'cp')
 112       CALL QNLSET('W1ANA',.TRUE.)
 113       CALL QNLSET('W2NUM',.TRUE.)
 114       CALL QNLSET('W2STF',.FALSE.)
 115       if (index(gridname,'none').eq.1) then
 116          call qnfilw(0,0)
 117       else
 118          qnerr=-1
 119          open(unit=2,status='old',file=gridname,
 120      .        form='unformatted',err=1)
 121          call QNREAD(2,1,qnerr)
 122  1       close(2)
 123          if (qnerr.ne.0) then
 124             write(*,*) 'Grid file problem: ',gridname
 125             if (qnerr.lt.0) then
 126                write(*,*) 'Grid file does not exist'
 127                write(*,*) 'Calculating and creating grid file'
 128                call qnfilw(0,0)
 129                open(unit=2,status='unknown',file=gridname,
 130      .              form='unformatted')
 131                call QNDUMP(2)
 132                close(2)
 133             else
 134                write(*,*) 'Existing grid file is inconsistent'
 135                if (qnerr.eq.1)
 136      .              write(*,*) 'Defined grid different'
 137                if (qnerr.eq.2)
 138      .              write(*,*) 'Heavy quark weight table different'
 139                if (qnerr.eq.3)
 140      .              write(*,*) 'Charm mass different'
 141                if (qnerr.eq.4)
 142      .              write(*,*) 'Bottom mass different'
 143                stop
 144             endif
 145          endif
 146       endif
 147       return
 148 *
 149       entry QCDNUM4pdf(nset)
 150 c      print *,'entering QCDNUMpdf',nset
 151       call GetAlfas(nset,alfas0,scale0)
 152 c      print *,alfas0,scale0
 153       Q2=scale0*scale0
 154       CALL QNRSET('ALFAS',alfas0)
 155       CALL QNRSET('ALFQ0',Q2)
 156       DO ix = 1,nx
 157          xx = XFROMIX(ix)
 158 c        print *,'calling parmPDF',ix
 159          call parmPDF(nset,xx,f)
 160 c        if(ix.lt.6) print *,nset,xx,f
 161          singlet=0.d0
 162          do i=1,nf0
 163             singlet=singlet+f(i)+f(-i)
 164 c           print  *,i,singlet
 165          end do
 166          gl=f(0)
 167          dm=f(1)-f(-1)
 168          um=f(2)-f(-2)
 169          dp=f(1)+f(-1)-singlet/dble(nf0)
 170          up=f(2)+f(-2)-singlet/dble(nf0)
 171          sp=f(3)+f(-3)-singlet/dble(nf0)
 172 c         print *,ix,iq0
 173          CALL QNPSET('SINGLET',ix,iq0,singlet)
 174          CALL QNPSET('GLUON',ix,iq0,gl)
 175          CALL QNPSET('DM',ix,iq0,DM)
 176          CALL QNPSET('UM',ix,iq0,UM)
 177          CALL QNPSET('DP',ix,iq0,DP)
 178          CALL QNPSET('UP',ix,iq0,UP)
 179          CALL QNPSET('SP',ix,iq0,SP)
 180       ENDDO
 181       CALL EVOLSG(iq0,1,nq)
 182       CALL EVOLNM('DM',iq0,1,nq)
 183       CALL EVOLNM('UM',iq0,1,nq)
 184       CALL EVLSEA4('dp',iq0,iqc,nq)
 185       CALL EVLSEA4('up',iq0,iqc,nq)
 186       CALL EVLSEA4('sp',iq0,iqc,nq)
 187
 188 c      print *,'calling evols - heavy...'
 189
 190 C--   Heavy quark evolution
 191
 192       CALL EVOLCP4('cp',iqc,nq)
 193 c
 194       call getnset(iset)
 195       call save_pdfqcd(iset)
 196       return
 197 *
 198       end
 199       subroutine EVLSEA4(name,IQ0,IQC,NQGRI)
 200       implicit none
 201
 202       CHARACTER*(*) name
 203       integer iq0,iqc,nqgri
 204       real*8 f34,f45,f43,f54,factor
 205       parameter(f34=1.d0/12.d0,f43=-1.d0/12.d0)
 206
 207 c      print *,iq0,iqc,nqgri
 208       If(IQ0.le.IQC)then
 209          CAll EVPLUS(name,IQ0,1,IQC)
 210          CALL QADDSI(name,IQC,f34)
 211          CAll EVPLUS(name,IQC,IQC,NQGRI)
 212       else if(IQ0.gt.IQC)then
 213          CAll EVPLUS(name,IQ0,IQC,NQGRI)
 214          CALL QADDSI(name,IQC,f43)
 215          CAll EVPLUS(name,IQC,1,IQC)
 216       end if
 217
 218       end
 219
 220
 221       subroutine EVOLCP4(name,IQC,NQGRI)
 222       implicit none
 223
 224       CHARACTER*(*) name
 225       integer iq0,iqc,nqgri
 226       real*8 f4
 227       parameter(f4=-1.d0/4.d0)
 228
 229 c
 230 c     First set to zero to avoid adding -1/4Singl at each iteration
 231 c
 232
 233       CAll QNPNUL(name)
 234       CALL QADDSI(name,IQC,f4)
 235       CAll EVPLUS(name,IQC,IQC,NQGRI)
 236
 237       end