LHAPDF/lhapdf5.5.1/src/wrapQCDNUM4.f

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