The VIRCOL website

=The LHAPDF home page= Theory = LHAPDF ready MC's = Downloads = online manual =

LHAPDF version 1 Users Guide

LHAPDF version 1 Users Guide 1

Abstract

The Les Houches Accord PDF (LHAPDF) interface package is designed to work with PDF sets. A PDF set can consist of many individual member PDFs. While the interpretation of the member PDFs depends on the particular set, the LHAPDF interface is specifically designed to accommodate PDFs with uncertainties. For PDFs with uncertainties the PDF set represents one ``fit'' to the data. The individual member of a PDF set are needed to calculate the PDF uncertainty of the observable.

All PDF sets are defined through external files. This means that a new set can be added by simply downloading its file while the LHAPDF interface code does not change. The evolution code is not part of LHAPDF. Currently, QCDNUM [1] is the default evolution code fully interfaced with LHAPDF. Other evolution codes can easily be interfaced with the LHAPDF.

Contents

1  Introduction
2  Analysis using PDF Uncertainties
    2.1  Method of Random Sampling
3  Interfacing LHAPDF with a Code
A  Installing LHAPDF
B  Examples
    B.1  Example 1: A PDF table
    B.2  Example 2: Calculating Uncertainties
    B.3  Example 3: A Convenient Wrap

1  Introduction

The Les Houches Accord Parton Density Function interface originated at the Les Houches 2001 workshop [2] in the PDF working group to enable the usage of Parton Density Functions with uncertainties in a uniform manner. When PDFs with uncertainties are considered, a ``fit'' to the data no longer is described by a single PDF. Instead in its most flexible implementation, a fit is represented by a PDF set consisting of many individual members. Calculating the observable for all the PDF members enables one to reconstruct the uncertainty on the observable. The LHAPDF interface was made with this in mind and manipulates PDF sets.

The LHAPDF interface can be viewed as a successor to PDFLIB [3]. Many improvements were added, to list some of the features:

  • The handling of PDF sets to enable PDF fits that include uncertainties.
  • All PDF sets are defined through external files in parametrized form. This means the files are compact compared to storing the PDFs in a grid format. Also new PDF sets can be defined by constructing the PDF defining files. The actual LHAPDF code does not have to be changed.
  • The LHAPDF code is modular and default choices can be easily altered.

Note that the current ``best fit'' PDFs can be viewed as sets with one member PDF and can be easily defined through the PDF set external file definition. Alternatively one can group these ``fits'' is single sets (e.g. MRST98 set) as they often represent best fits given different model assumptions and as such reflect theoretical modeling uncertainties.

2  Analysis using PDF Uncertainties

As mentioned before using LHAPDF to determine the PDF uncertainty on an observable will involve evaluating N different PDFs in a single set. How to use the N predictions of the observable depends on the method used in the PDF set for propagation of errors. During the fitting of the PDFs many approximations can be made. An overview is given in the Les Houches proceedings [2]. However, in this section we will only discuss the method of error propagation and how to reconstruct the PDF uncertainty of the observable from the available sets.

2.1  Method of Random Sampling

The method of random sampling [4] makes it conceptually easy to perform analysis involving PDF uncertainties. It draws uniformly N PDFs functionals from the PDF probability density function. This leads to an intuitive procedure to construct the PDF uncertainty for the observable. For each of the PDFs F we make a prediction of the observable xt. This will build up a scatter plot representation of the uncertainty density function of the observable
PPDFscatter(xe) =  1

N

N
å
i=1 
 d(xe-xt(Fi)) .
(1)

For a more ``graphical'' prediction one can ``smoothen out'' the scatter representation by a smearing function S(x)

PPDFsmear(xe)
=
ó
õ
d x PPDFscatter(x) S(xe-x)
=
 1

N

N
å
i=1 
S(xe-xt(Fi)) .
(2)

The simplest example is a histogram with a bin width 2×D

PPDFhisto(xe)=  1

N

N
å
i=1 
Q(xt(Fi)-(xe-D))Q((xe-D)-xt(Fi)) .
(3)

If we compare to an experiment the smoothening function itself S can chosen to be the experimental uncertainty function for the observable.

From the above probability density function one can calculate confidence level intervals giving the probability the observable is bracketed between x0 < xe < x1
CL(x0;x1)= ó
õ
x1

x0 
d x PPDF(x) .
(4)

3  Interfacing LHAPDF with a Code

The interface of LHAPDF with an external code is easy. We will describe the basic steps sufficient for most applications. The web site contains more detailed information about additional function calls. The function calls described here will not be altered in any way in future versions. Including the LHAPDF code into a program involves three steps:

  1. First one has to setup the LHAPDF interface code:

    call InitPDFset(name)

     . It is called only once at the beginning of the code. The string variable name is the file name of the external PDF file that defines the PDF set. For the default evolution code QCDNUM it will either calculate or read from file the LO/NLO splitting function weights. The calculation of the weights might take some time depending on the chosen grid size. However, after the first run a grid file is created. Subsequent runs will use this file to read in the weights so that the lengthy calculation of these weights is avoided. The file depends on the grid parameters and flavor thresholds. This means different PDF sets can have different grid files. The name of the grid file is specified in the PDF setup file.

  2. To use a individual PDF member it has to be initialized:

    call InitPDF(mem)

     . The integer mem specifies the member PDF number. This routine only needs to be called when changing to a new PDF member. The time needed to setup depends on the evolution code used. For QCDNUM the grid size is the determining factor. Note that mem=0 selects the ``best fit'' PDF.

  3. Once the PDF member is initialized one can call the evolution codes which will use the selected member.

    The double precision function call

    function alphasPDF(Q)

    returns the value of aS(Q) at double precision scale Q. Note that its value can change between different PDF members.

    The subroutine call

    call evolvePDF(x,Q,f)

    returns the PDF momentum densities f (i.e. x× PDF number density) at double precision momentum fraction x and double precision scale Q. The double precision array f(-6:6) will contain the momentum PDFs using the labeling convention of table 1.

                   

    -6 -5 -4 -3 -2 -1  0 1 2 3 4 5 6
      P `t `b `c `s `u `d g d u s c b t
    `P  t b c s u d g `d `u `s `c `b `t

    Table 1: The flavor enumeration convention used in the LHAPDF interface.

    As long as the member PDF is not changed (by the call InitPDF of step 2) the evolution calls will always use the same PDF member.

A few additional calls can be useful:

  • To get the number of PDF members in the set:

    call numberPDF(Nmem).

    The integer Nmem will contain the number of PDF members (excluding the special ``best fit'' member, i.e. the member numbers run from 0 to Nmem).

  • Optionally the different PDF members can have weights which are obtained by:

    call weightPDF(wgt).

    The double precision variable wgt is for unweighted PDFs set to 1/Nmem such that the sum of all PDF member weights is unity. For weighted sets the use of the weights has to be defined by the method description.

  • To get the evolution order of the PDFs:

    call GetOrderPDF(order).

    The integer variable order is 0 for Leading Order, 1 for Next-to-Leading Order, etc.

  • To get the evolution order of aS:

    call GetOrderAs(order).

    The integer variable order is 0 for Leading Order, 1 for Next-to-Leading Order, etc.

  • It is possible that during the PDF fitting the renormalization scale was chosen different from the factorization scale. The ratio of the renormalization scale over the factorization scale used in the ``fit'' can be obtained by

    call GetRenFac(muf).

    The double precision variable muf contains the ratio. Usually muf is equal to unity.

  • To get a description of the PDF set:

    call GetDesc().

    This call will print the PDF description to the default output stream

  • The quark masses can be obtained by:

    call GetQmass(nf,mass).

    The mass mass is returned for quark flavor nf. The quark masses are used in the aS evolution.

  • The flavor thresholds in the PDF evolution can be obtained by:

    call GetThreshold(nf,Q).

    The flavor threshold Q is returned for flavor nf. If Q=-1d0 the flavor is not in the evolution (e.g. the top quark is usually not included in the evolution). If Q=0d0 the flavor is parametrized at the parametrization scale. For positive non-zero values of Q the value is set to the flavor threshold at which the PDF starts to evolve.

  • The call returns the number of flavors used in the PDF:

    call GetNf(nfmax).

    Usually the returned value for nfmax is equal to five as the top quark is usually not considered in the PDFs.

A  Installing LHAPDF

Installing LHAPDF involves 3 steps:

  1. Download the LHAPDF files. Untarring the file will create the directory LHAPDFv1 which contains all appropriate files and this manual.
  2. Download the evolution code. For example downloading and untarring the QCDNUM file will create the directory QCDNUM with the appropriate files together with the makefile to make the library file libLHAPDF.a. The library file can be linked with other code to make LHAPDF available. Note that it might be necesarry to edit the pathnames in the makefile.
  3. Download the PDF external files. Untarring will create the directory PDFsets which contains the PDF set files (extension .LHpdf).

B  Examples

The example tar file can be downloaded from the website. Untarring the file will create the directory Examples which contains the source code of the 3 examples together with the makefile to generate the executables. Note that the pathnames in the makefile have to be edited for proper execution. When running an example for the first time with the evolution code QCDNUM a grid has to be calculated. This might take some time depending on the grid size. Running with the PDF set after that the initialization time is much shorter as everything can be read in from an external file.

B.1  Example 1: A PDF table

The executable example.x is made by the command make 1. The first example code simply runs through all 100 PDF members of the Alekhin set [6], printing out aS(MZ) and the gluon PDF momentum densities at a few (x,Q) points:

      program example1
      implicit real*8(a-h,o-z)
      character name*50
      real*8 f(-6:6)
*
      name='../PDFsets/Alekhin_100.LHpdf'
      call InitPDFset(name)
*
      QMZ=91.71d0
      write(*,*)
      call numberPDF(N)
      do i=1,N
         write(*,*) '---------------------------------------------'
         call InitPDF(i)
         write(*,*) 'PDF set ',i
         write(*,*)
         a=alphasPDF(QMZ)
         write(*,*) 'alpha_S(M_Z) = ',a
         write(*,*)
         write(*,*) 'x*Gluon'
         write(*,*) '   x     Q=10 GeV     Q=100 GeV    Q=1000 GeV'
         do x=0.01d0,0.095d0,0.01d0
            Q=10d0
            call evolvePDF(x,Q,f)
            g1=f(0)
            Q=100d0
            call evolvePDF(x,Q,f)
            g2=f(0)
            Q=1000d0
            call evolvePDF(x,Q,f)
            g3=f(0)
            write(*,*) x,g1,g2,g3
         enddo
      enddo
*
      end

B.2  Example 2: Calculating Uncertainties

The executable example.x is made by the command make 2. The second example calculates the average value of the gluon PDF and its standard deviation for several x values at Q=100 GeV using the Botje set [7]. Note that the program calculates the quantities in two manners. First by placing the loop over the PDF members inside the loop over the parton momentum fractions. Second, by exchanging this order and storing the first and second moments in an array depending on the x value. The second method is considerably faster in computing time as each PDF member is only initialized once.

Also calculated is the correlation of the gluon momentum PDF between x=0.001 and x=0.01 and the correlation between the strange quark momentum PDF at x=0.001 and aS(Q) at Q=10 GeV.

      program example2
      implicit real*8(a-h,o-z)
      character*32 name
      real*8 f(-6:6),mom1(9),mom2(9)
*
      Q=100d0
      name='../PDFsets/Botje_100.LHpdf'
      call InitPDFset(name)
*
      call numberPDF(Nmem)
      write(*,*)
      write(*,*) 'Calculating the gluon momentum PDF average <g>'
      write(*,*) 'and standard deviaton SD(g) for several x values'
      write(*,*) 'at Q=100 GeV'
      write(*,*)
      write(*,*) '1. The slow way:'
      write(*,*) 
      write(*,*) '   x       <g>         SD(g)'
      do x=0.01d0,0.095d0,0.01d0
         gmom1=0d0
         gmom2=0d0
         do i=1,Nmem
            call InitPDF(i)
            call evolvePDF(x,Q,f)
            gmom1=gmom1+f(0)
            gmom2=gmom2+f(0)**2
         enddo
         av=gmom1/Nmem
         sd=sqrt(gmom2/Nmem-av**2)
         write(*,*) x,av,sd
      enddo
      write(*,*) 
      write(*,*) '2. The fast way:'
      write(*,*) 
      write(*,*) '   x       <g>         SD(g)'
      do i=1,9
         mom1(i)=0d0
         mom2(i)=0d0
      enddo
      do i=1,Nmem
         call InitPDF(i)
         ic=0
         do x=0.01d0,0.095d0,0.01d0
            call evolvePDF(x,Q,f)
            ic=ic+1
            mom1(ic)=mom1(ic)+f(0)
            mom2(ic)=mom2(ic)+f(0)**2
         enddo
      enddo
      ic=0
      do x=0.01d0,0.095d0,0.01d0
         ic=ic+1
         av=mom1(ic)/Nmem
         sd=sqrt(mom2(ic)/Nmem-av**2)
         write(*,*) x,av,sd
      enddo
*
      Q=10d0
      x1=0.001d0
      x2=0.01d0
      write(*,*)
      write(*,*) 'Calculating the normalized correlation coefficient'
      write(*,*) '<g1g2> between g(x=0.001) and g(x=0.01) and'
      write(*,*) '<sAlpha> between the strange quark momentum PDF'
      write(*,*) 'at x=0.001 and alpha_S(Q) for Q=10 GeV'
      write(*,*)
      avg1=0d0
      avg2=0d0
      avs=0d0
      avAs=0d0
      sdg1=0d0
      sdg2=0d0
      sds=0d0
      sdAs=0d0
      Cg1g2=0d0
      CsAs=0d0
      j=3
      do i=1,Nmem
         call InitPDF(i)
         As=alphasPDF(Q)
         call evolvePDF(x1,Q,f)
         g1=f(0)
         s=f(j)
         call evolvePDF(x2,Q,f)
         g2=f(0)
         avAs=avAs+As
         avg1=avg1+g1
         avg2=avg2+g2
         avs=avs+s
         sdAs=sdAs+As**2
         sdg1=sdg1+g1**2
         sdg2=sdg2+g2**2
         sds=sds+s**2
         CsAs=CsAs+s*As
         Cg1g2=Cg1g2+g1*g2
      enddo
      avAs=avAs/Nmem
      avs=avs/Nmem
      avg1=avg1/Nmem
      avg2=avg2/Nmem
      sdAs=sdAs/Nmem-avAs**2
      sds=sds/Nmem-avs**2
      sdg1=sdg1/Nmem-avg1**2
      sdg2=sdg2/Nmem-avg2**2
      CsAs=CsAs/Nmem-avs*avAs
      Cg1g2=Cg1g2/Nmem-avg1*avg2
      write(*,*) '<g1g2>    = ',Cg1g2/sqrt(sdg1*sdg2)
      write(*,*) '<sAlpha> = ',CsAs/sqrt(sds*sdAs)
      end

B.3  Example 3: A Convenient Wrap

The executable example.x is made by the command make 3. The third example prints out the same as the first example. However it uses a wrapper around the LHAPDF code. The single call wrapper takes care of the initialization calls. Note that the number of PDF members is only known after the first call to the subroutine parden. Also the aS value is only known after the first call to the subroutine parden for a specific member PDF.

      program example3
      implicit real*8(a-h,o-z)
      real*8 f(-6:6)
*
      QMZ=91.71d0
      Nmem=100
      write(*,*)
      do i=1,Nmem
         write(*,*) '---------------------------------------------'
         write(*,*) 'PDF set ',i
         write(*,*)
         write(*,*) 'x*Gluon'
         write(*,*) '   x     Q=10 GeV     Q=100 GeV    Q=1000 GeV'
         do x=0.01d0,0.095d0,0.01d0
            Q=10d0
            call parden(x,Q,f,i)
            g1=f(0)
            Q=100d0
            call parden(x,Q,f,i)
            g2=f(0)
            Q=1000d0
            call parden(x,Q,f,i)
            g3=f(0)
            write(*,*) x,g1,g2,g3
         enddo
         a=alphasPDF(QMZ)
         write(*,*)
         write(*,*) 'alpha_S(M_Z) = ',a
         write(*,*)
      enddo
*
      end
*
      subroutine parden(x,Q,f,imem)
      implicit none
      character*32 name/'../PDFsets/Alekhin_100.LHpdf'/
      real*8 x,Q,f(-6:6),alfas
      integer imem,init/0/,lmem/-1/
      save init,lmem
*
      if (init.eq.0) then
         init=1
         call InitPDFset(name)
      endif
      if (imem.ne.lmem) then
         lmem=imem
         call InitPDF(lmem)
      endif
      call evolvePDF(x,Q,f)
      return
*
      end

References

[1]
M. Botje, QCDNUM version 16.12, preprint ZEUS-97-066,
http://www.nikhef.nl/ ~ h24/qcdnum.
[2]
S. Alekhin, M. Botje, W. Giele, J. Pumplin, F. Olness, G. Salam and A. Vogt,
Physics at TeV Colliders II Workshop, Les Houches, France, May 2001.
[3]
H. Plothow-Besch, Comput. Phys. Commun. 75, 396 (1993).
[4]
W. Giele, S. Keller and D. Kosower, hep-ph/0104052;
W. Giele and S. Keller, hep-ph/9803393 [Phys. Rev. D58, 1997, 094023].
[5]
W. Giele and S. Keller, hep-ph/0104053.
[6]
S. I. Alekhin, hep-ph/0011002, [Phys. Rev. D63, 2001, 094022].
[7]
M. Botje, hep-ph/9912439, [Eur. Phys. J. C14, 2000, 285].


Footnotes:

1W. Giele, Fermilab, giele@fnal.gov


File translated from TEX by TTH, version 3.06.
On 5 Apr 2002, 11:34.


Comments and/or Problems
=The LHAPDF home page= Theory = LHAPDF ready MC's = Downloads = online manual =