LHAPDF/lhapdf5.5.1/src/wrapQCDNUM.f

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