LHAPDF/lhapdf5.3.1/wrapQCDNUM.f

   1       subroutine QCDNUMevolve(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,isetlast
  18 *
  19       call getnset(iset)
  20 c fix 14.7.2006 to speed up sue of QCDNUM
  21       if(iset.ne.isetlast) then
  22         call get_pdfqcd(iset)
  23         isetlast = iset
  24       endif
  25 c
  26       Q2=Q*Q
  27       nf=5
  28       if (Q2.lt.tb2) nf=4
  29       if (Q2.lt.tc2) nf=3
  30       pdf(0)=QPDFXQ('GLUON'  ,x,Q2,iflag)
  31       singlet=QPDFXQ('SINGLET',x,Q2,iflag)
  32       dm= QPDFXQ('DM',x,Q2,IFLAG)
  33       um= QPDFXQ('UM',x,Q2,IFLAG)
  34       dp= QPDFXQ('DP',x,Q2,IFLAG)
  35       up= QPDFXQ('UP',x,Q2,IFLAG)
  36       sp= QPDFXQ('SP',x,Q2,IFLAG)
  37       ub=0.5d0*(up-um+singlet/dble(nf))
  38       db=0.5d0*(dp-dm+singlet/dble(nf))
  39       sb=0.5d0*(sp+singlet/dble(nf))
  40       cb=0.d0
  41       if (nf.ge.4) then
  42          cp= QPDFXQ('CP',X,Q2,IFLAG)
  43          cb=0.5d0*(cp+singlet/dble(nf))
  44       end if
  45       bb=0.d0
  46       if (nf.ge.5) then
  47          bp= QPDFXQ('BP',X,Q2,IFLAG)
  48          bb=0.5d0*(bp+singlet/dble(nf))
  49       end if
  50
  51       pdf(1)=dm+db
  52       pdf(2)=um+ub
  53       pdf(3)=sb
  54       pdf(4)=cb
  55       pdf(5)=bb
  56       pdf(6)=0d0
  57       pdf(-1)=db
  58       pdf(-2)=ub
  59       pdf(-3)=sb
  60       pdf(-4)=cb
  61       pdf(-5)=bb
  62       pdf(-6)=0d0
  63       return
  64 *
  65       entry QCDNUMalfa(alfas,Q)
  66       Q2=Q*Q
  67       nf=6
  68       if (Q2.lt.mt2) nf=5
  69       if (Q2.lt.mb2) nf=4
  70       if (Q2.lt.mc2) nf=3
  71       alfas=QALFAS(Q2,Qlam,nf,iflag)
  72       return
  73 *
  74       entry QCDNUMread(nset)
  75 c      read(1,*) s1
  76 c      print *,s1
  77 c         gridname='large.grid'
  78 c         nx=400
  79 c         xmin=1d-6
  80 c         xmax=1d0
  81 c         nq=112
  82 c         qmin=1d0
  83 c         qmax=1d10
  84 c         read(1,*) dummy
  85 *      else
  86          read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
  87 c*      endif
  88 c      iset = nset
  89       return
  90 *
  91       entry QCDNUMinit(nset,Eorder,Q2fit)
  92       call qninit
  93       call QNISET('ORDER',Eorder+1)
  94       call grxdef(nx,xmin)
  95       call grqdef(nq,qmin,qmax)
  96       call grqinp(Q2fit,1)
  97       call getQMassM(nset,4,mc)
  98       mc2=mc*mc
  99       CALL QNRSET('CMASS',mc)
 100       CALL QNRSET('MCALF',mc)
 101       call getThresholdM(nset,4,tc)
 102       tc2=tc*tc
 103       call GRQINP(tc2,1)
 104       CALL GRQINP(tc2-1.0d-4,1)
 105       call getQMassM(nset,5,mb)
 106       mb2=mb*mb
 107       CALL QNRSET('BMASS',mb)
 108       CALL QNRSET('MBALF',mb)
 109       call getThresholdM(nset,5,tb)
 110       tb2=tb*tb
 111       call GRQINP(tb2,1)
 112       CALL GRQINP(tb2-1.0d-4,1)
 113       call getQMassM(nset,6,mt)
 114       mt2=mt*mt
 115       CALL QNRSET('TMASS',mt)
 116       CALL QNRSET('MTALF',mt)
 117       call qthres(tc2,tb2)
 118       iq0=IQFROMQ(Q2fit)
 119       iqc=IQFROMQ(tc2)
 120       iqb=IQFROMQ(tb2)
 121       nf0= NFLGET(iq0)
 122 c      print *,nf0
 123       CALL QNBOOK(2,'dm')
 124       CALL QNBOOK(3,'um')
 125       CALL QNBOOK(4,'dp')
 126       CALL QNBOOK(5,'up')
 127       CALL QNBOOK(6,'sp')
 128       CALL QNBOOK(7,'cp')
 129       CALL QNBOOK(8,'bp')
 130       CALL QNLSET('W1ANA',.TRUE.)
 131       CALL QNLSET('W2NUM',.TRUE.)
 132       CALL QNLSET('W2STF',.FALSE.)
 133       if (index(gridname,'none').eq.1) then
 134          call qnfilw(0,0)
 135       else
 136          qnerr=-1
 137          open(unit=2,status='old',file=gridname,
 138      .        form='unformatted',err=1)
 139          call QNREAD(2,1,qnerr)
 140  1       close(2)
 141          if (qnerr.ne.0) then
 142             write(*,*) 'Grid file problem: ',gridname
 143             if (qnerr.lt.0) then
 144                write(*,*) 'Grid file does not exist'
 145                write(*,*) 'Calculating and creating grid file'
 146                call qnfilw(0,0)
 147                open(unit=2,status='unknown',file=gridname,
 148      .              form='unformatted')
 149                call QNDUMP(2)
 150                close(2)
 151             else
 152                write(*,*) 'Existing grid file is inconsistent'
 153                if (qnerr.eq.1)
 154      .              write(*,*) 'Defined grid different'
 155                if (qnerr.eq.2)
 156      .              write(*,*) 'Heavy quark weight table different'
 157                if (qnerr.eq.3)
 158      .              write(*,*) 'Charm mass different'
 159                if (qnerr.eq.4)
 160      .              write(*,*) 'Bottom mass different'
 161                stop
 162             endif
 163          endif
 164       endif
 165       return
 166 *
 167       entry QCDNUMpdf(nset)
 168 c      print *,'entering QCDNUMpdf',nset
 169       call GetAlfas(nset,alfas0,scale0)
 170 c      print *,alfas0,scale0
 171       Q2=scale0*scale0
 172       CALL QNRSET('ALFAS',alfas0)
 173       CALL QNRSET('ALFQ0',Q2)
 174       DO ix = 1,nx
 175          xx = XFROMIX(ix)
 176 c        print *,'calling parmPDF',ix
 177          call parmPDF(nset,xx,f)
 178 c        if(ix.lt.6) print *,nset,xx,f
 179          singlet=0.d0
 180          do i=1,nf0
 181             singlet=singlet+f(i)+f(-i)
 182 c           print  *,i,singlet
 183          end do
 184          gl=f(0)
 185          dm=f(1)-f(-1)
 186          um=f(2)-f(-2)
 187          dp=f(1)+f(-1)-singlet/dble(nf0)
 188          up=f(2)+f(-2)-singlet/dble(nf0)
 189          sp=f(3)+f(-3)-singlet/dble(nf0)
 190 c      print *,'calling qpnset sg...',ix,iq0,singlet
 191          CALL QNPSET('SINGLET',ix,iq0,singlet)
 192 c      print *,'calling qpnset gl...',ix,iq0,gl
 193          CALL QNPSET('GLUON',ix,iq0,gl)
 194 c      print *,'calling qpnset dm...',ix,iq0,dm
 195          CALL QNPSET('DM',ix,iq0,DM)
 196 c      print *,'calling qpnset um...',ix,iq0,um
 197          CALL QNPSET('UM',ix,iq0,UM)
 198 c      print *,'calling qpnset dp...',ix,iq0,dp
 199          CALL QNPSET('DP',ix,iq0,DP)
 200 c      print *,'calling qpnset up...',ix,iq0,up
 201          CALL QNPSET('UP',ix,iq0,UP)
 202 c      print *,'calling qpnset sp...',ix,iq0,sp
 203          CALL QNPSET('SP',ix,iq0,SP)
 204       ENDDO
 205 c      print *,'calling evols...'
 206       CALL EVOLSG(iq0,1,nq)
 207 c      print *,'calling evols dm...'
 208       CALL EVOLNM('DM',iq0,1,nq)
 209 c      print *,'calling evols um...'
 210       CALL EVOLNM('UM',iq0,1,nq)
 211 c      print *,'calling evols dp...',iq0,iqc,iqb,nq
 212       CALL EVLSEA('dp',iq0,iqc,iqb,nq)
 213 c      print *,'calling evols up...'
 214       CALL EVLSEA('up',iq0,iqc,iqb,nq)
 215 c      print *,'calling evols sp...'
 216       CALL EVLSEA('sp',iq0,iqc,iqb,nq)
 217
 218 c      print *,'calling evols - heavy...'
 219
 220 C--   Heavy quark evolution
 221
 222       CALL EVOLCP('cp',iqc,iqb,nq)
 223       CALL EVOLBP('bp',iqb,nq)
 224 c
 225       call getnset(iset)
 226       call save_pdfqcd(iset)
 227       return
 228 *
 229       end
 230       subroutine EVLSEA(name,IQ0,IQC,IQB,NQGRI)
 231       implicit none
 232
 233       CHARACTER*(*) name
 234       integer iq0,iqc,iqb,nqgri
 235       real*8 f34,f45,f43,f54,factor
 236       parameter(f34=1.d0/12.d0,f45=1.d0/20.d0,f43=-1.d0/12.d0,
 237      >          f54=-1.d0/20.d0)
 238
 239
 240       If(IQ0.le.IQC)then
 241          CAll EVPLUS(name,IQ0,1,IQC)
 242          CALL QADDSI(name,IQC,f34)
 243          CAll EVPLUS(name,IQC,IQC,IQB)
 244          CALL QADDSI(name,IQB,f45)
 245          CAll EVPLUS(name,IQB,IQB,NQGRI)
 246       else if(IQ0.gt.IQC.and.IQ0.le.IQB)then
 247          CAll EVPLUS(name,IQ0,IQC,IQB)
 248          CALL QADDSI(name,IQC,f43)
 249          CAll EVPLUS(name,IQC,1,IQC)
 250          CALL QADDSI(name,IQC,f34)
 251          CALL QADDSI(name,IQB,f45)
 252          CAll EVPLUS(name,IQB,IQB,NQGRI)
 253       else if(IQ0.gt.IQB)then
 254          CAll EVPLUS(name,IQ0,IQB,NQGRI)
 255          CALL QADDSI(name,IQB,f54)
 256          CAll EVPLUS(name,IQB,IQC,IQB)
 257          CALL QADDSI(name,IQB,f45)
 258          CALL QADDSI(name,IQC,f43)
 259          CAll EVPLUS(name,IQC,1,IQC)
 260          CALL QADDSI(name,IQC,f34)
 261       end if
 262
 263       end
 264
 265
 266       subroutine EVOLCP(name,IQC,IQB,NQGRI)
 267       implicit none
 268
 269       CHARACTER*(*) name
 270       integer iq0,iqc,iqb,nqgri
 271       real*8 f4,f45
 272       parameter(f4=-1.d0/4.d0,f45=1.d0/20.d0)
 273
 274 c
 275 c     First set to zero to avoid adding -1/4Singl at each iteration
 276 c
 277
 278       CAll QNPNUL(name)
 279       CALL QADDSI(name,IQC,f4)
 280       CAll EVPLUS(name,IQC,IQC,IQB)
 281       CALL QADDSI(name,IQB,f45)
 282       CAll EVPLUS(name,IQB,IQB,NQGRI)
 283
 284       end
 285
 286       subroutine EVOLBP(name,IQB,NQGRI)
 287       implicit none
 288
 289       CHARACTER*(*) name
 290       integer iqb,nqgri
 291       real*8 f5
 292       parameter(f5=-1.d0/5.d0)
 293
 294 c
 295 c     First set to zero to avoid adding -1/5Singl at each iteration
 296 c
 297
 298       CAll QNPNUL(name)
 299       CALL QADDSI(name,IQB,f5)
 300       CAll EVPLUS(name,IQB,IQB,NQGRI)
 301
 302       end
 303 c
 304 c routine to swap in/out the array PDFQCD(MXX,MQ2,0:10) in QCDNUM as required
 305 c by nset
 306 c
 307       subroutine get_pdfqcd(nset)
 308 c
 309       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 310 c
 311       include 'parmsetup.inc'
 312 c
 313       PARAMETER ( MXX = 410 )
 314       PARAMETER ( MQ2 =  120 )
 315 C--   Do not set the following parameter to zero!
 316       PARAMETER ( NDFMAX = 20)
 317 c
 318       COMMON/QCPASS/
 319      +ALPHA0, Q0ALFA, ASLAST, QALAST,
 320      +ALFASQ(MQ2), ALFAPQ(MQ2), ALFA2Q(MQ2),
 321      +DELUP(MQ2), DELDN(MQ2), PDFQCD(MXX,MQ2,0:10),
 322      +FNSQCD(MXX,MQ2),DNSQCD(MXX,MQ2),
 323      +FSIQCD(MXX,MQ2),DSIQCD(MXX,MQ2),
 324      +FGLQCD(MXX,MQ2),DGGQCD(MXX,MQ2),
 325      +FSTORE(MXX,MQ2,31:30+NDFMAX),IDFAST(7,30),NDFAST,
 326      +MARKFF(MXX,MQ2),MARKFH(MXX,MQ2),MARKQQ(MQ2),
 327      +ISTFID(31:30+NDFMAX),IPDFID(31:30+NDFMAX),IEALFA(MQ2),
 328      +IQL_LAST(10),IQ0_LAST(10),IQH_LAST(10)
 329 c
 330       real*8  PDFSAV(MXX,MQ2,0:10,nmxset)
 331       save PDFSAV
 332 c
 333       do k=0,10
 334         do j=1,mq2
 335           do i=1,mxx
 336              PDFQCD(i,j,k)=PDFSAV(i,j,k,nset)
 337           enddo
 338         enddo
 339       enddo
 340 c
 341       return
 342 c
 343       entry save_pdfqcd(nset)
 344 c
 345       do k=0,10
 346         do j=1,mq2
 347           do i=1,mxx
 348              PDFSAV(i,j,k,nset)=PDFQCD(i,j,k)
 349           enddo
 350         enddo
 351       enddo
 352 c
 353       return
 354       end
 355
 356