NNPDF is hosted by Hepforge, IPPP Durham

PDF Sets: Tutorial, using NNPDF sets

In this tutorial we show how to compute correctly the central value and the 1-sigma error band for the gluon PDF using a generic NNPDF set.

We propose the same code in C++ and FORTRAN languages, by using the LHAPDF interface.

Downloading and running the tutorial

1. Download
wget http://nnpdf.hepforge.org/nnpdf_tutorial.tar.gz
2. Unpack
tar -xvzf nnpdf_tutorial.tar.gz
3. Compile
cd nnpdf_tutorial
make
4. Run
# run C++ example
./tutorial_cxx

# run Fortran example
./tutorial_fortran

C++ code example:

/**
 * NNPDF tutorial: How to compute the gluon PDF central value 
 * and 1-sigma error band in C++, using LHAPDF interface.
 *
 * Compile with:
 * g++ tutorial.cxx -o tutorial_cxx -I $(shell lhapdf-config --incdir)
 * -L $(shell lhapdf-config --libdir) -lLHAPDF
 * 
 * Author: The NNPDF Collaboration - 2012
 */

#include "iostream"
#include "cmath"
#include "LHAPDF/LHAPDF.h"
using namespace std;

#define NxPoints 50

// Computes the standard deviation
double ComputeStdDev(int n, double *x, double cv)
{
  double sum = 0;
  for (int i = 0; i < n; i++)
    sum += pow(x[i] - cv, 2.0);

  sum /= n-1;
  return sqrt(sum);
}

int main()
{ 
  // Declaring variables
  const double Q = pow(2.0, 0.5);
  const double xmin = log(1e-5);
  const double xmax = log(1.0);
  const double delta = (xmax - xmin)/NxPoints;
    
  // Loading the NNPDF set
  LHAPDF::initPDFSet("NNPDF21_100.LHgrid");
  const int NumberReplicas = LHAPDF::numberPDF();

  // Getting the central values
  for (int i = 0; i < NxPoints; i++)
    {
      LHAPDF::initPDF(0); // <- Replica 0 == Central Value (AVG)
      const double x = exp(xmin + i*delta); 
      const double xgCV = LHAPDF::xfx(x,Q,0);      

      // Computing 1-sigma error band
      double ErrorBand[NumberReplicas];
      for (int n = 1; n <= NumberReplicas; n++)
			   {
			     LHAPDF::initPDF(n);
			     ErrorBand[n-1] = LHAPDF::xfx(x,Q,0);
			   }
      const double xgError = ComputeStdDev(NumberReplicas,ErrorBand,xgCV);

      // Printing results
      cout.precision(5);
      cout << scientific << "x = " << x << "\txg(x,Q) = (" << xgCV
           << "\t+/-\t" << xgError << ")" << endl;
    }
  
  return 0;
}

FORTRAN code example:

!
!     NNPDF tutorial: How to compute the gluon PDF central value
!     and 1-sigma error band in Fortran, using LHAPDF interface.
!
!     Compile with:
!     gfortran tutorial.f -o tutorial_fortran 
!     -L $(shell lhapdf-config --libdir) -lLHAPDF
!
!     Author: The NNPDF Collaboration - 2012
!
      program tutorial
      implicit none

!     Declaring variables
      integer i, n
      integer NxPoints
      integer NumberReplicas
      double precision Q, x
      double precision xpdfCV(-6:6)
      double precision xpdfError(-6:6) 
      double precision xpdfErr
      double precision xmin, xmax
      double precision delta
      double precision ErrorBand(1000)
      double precision ComputeStdDev
      
!     Initializing variables
      NxPoints = 50
      Q = 2.d0**0.5
      xmin = dlog(1d-5)
      xmax = dlog(1.d0)
      delta = (xmax - xmin)/NxPoints

!     Loading the NNPDF set
      call InitPDFsetByName("NNPDF21_100.LHgrid")
      call numberPDF(NumberReplicas)


!     Getting the central values
      do i = 0, NxPoints-1

         ! Replica 0 == Central Value (AVG)
         call InitPDF(0)
         x = exp(xmin + i*delta)
         call evolvePDF(x,Q,xpdfCV)

         ! Computing 1-sigma error band
         do n = 1, NumberReplicas
            call InitPDF(n)
            call evolvePDF(x,Q,xpdfError)
            ErrorBand(n) = xpdfError(0)
         enddo
         
         xpdfErr = ComputeStdDev(NumberReplicas,
     1        ErrorBand, xpdfCV(0))

         ! Printing results
         write(6,441) "x = ",x, " xg(x,Q) = (",
     1        xpdfCV(0)," +/- ",xpdfErr,")"
      enddo

 441  format(a,f15.10,a,f15.10,a,f15.10,a)

      stop
      end

!     Computes the standard deviation
      function ComputeStdDev(n,x,cv)
      integer i, n
      double precision sum, cv
      double precision x(n)
      double precision ComputeStdDev

      sum = 0.d0
      do i= 1, n
         sum = sum + (x(i)-cv)**2
      enddo
      
      sum = sum/(n-1)
      ComputeStdDev = sqrt(sum)

      return 
      end