NNPDFDriver.f

Go to the documentation of this file.

***
*
* NNPDF Fortran Driver
*
* Stefano Carrazza for the NNPDF Collaboration
* email: [email protected]
*
* February 2013
*
* Usage:
*
*  NNPDFDriver("gridname.LHgrid");
*
*  NNinitPDF(0); // select replica [0,Mem]
*
*  NNevolvePDF(x,Q,pdf); // -> returns double array (-6,7)
*     
      subroutine nnpdfdriver(gridfilename,lenfilename)
      implicit none
      integer nfl,nx,nq2,mem,rep,lenfilename
      double precision alphas
      double precision xgrid(100),logxgrid(100)
      double precision q2grid(50),logq2grid(50)
      double precision pdfgrid(0:100,14,100,50)
      logical hasphoton
      common /nnpdf/nfl,nx,nq2,mem,rep,hasphoton,alphas,xgrid,logxgrid,
     1     q2grid,logq2grid,pdfgrid
 
      character(len=lenfilename) gridfilename
*
      nfl = 13
      nx = 100
      nq2 = 50
      mem = 1
      rep = 0
      alphas = 0
*
*     Logo
      write(6,*) " ****************************************"
      write(6,*) ""
      write(6,*) "      NNPDFDriver version 1.0.3"
      write(6,*) "  Grid: ", gridfilename
      write(6,*) " ****************************************"
 
      call readpdfset(gridfilename,lenfilename)
      
      end subroutine
 
      subroutine nninitpdf(irep)
      implicit none
      integer irep
 
      integer nfl,nx,nq2,mem,rep
      double precision alphas
      double precision xgrid(100),logxgrid(100)
      double precision q2grid(50),logq2grid(50)
      double precision pdfgrid(0:100,14,100,50)
      logical hasphoton
      common /nnpdf/nfl,nx,nq2,mem,rep,hasphoton,alphas,xgrid,logxgrid,
     1     q2grid,logq2grid,pdfgrid
 
      if (irep.gt.mem.or.irep.lt.0d0) then
         write(6,*) "Error: replica out of range [0,",mem,"]"
      else
         rep = irep
      endif
      
      end subroutine
 
      subroutine readpdfset(gridfilename,lenfilename)
      implicit none
      integer i,ix,iq,fl,imem,lenfilename
      character(len=lenfilename) gridfilename
      character(len=100) line
*
      integer nfl,nx,nq2,mem,rep
      double precision alphas
      double precision xgrid(100),logxgrid(100)
      double precision q2grid(50),logq2grid(50)
      double precision pdfgrid(0:100,14,100,50)
      logical hasphoton
      common /nnpdf/nfl,nx,nq2,mem,rep,hasphoton,alphas,xgrid,logxgrid,
     1     q2grid,logq2grid,pdfgrid
*
      open (20, file=gridfilename, status='OLD')
*     Read header
      do i=1,1000
         read(20,*) line
         if (line(1:14).eq.'Parameterlist:') then
            read(20,*) line, mem, line, alphas                        
            exit
         endif
      enddo
 
*     Select driver
      do i=1,1000
         read(20,*) line
         if (line(1:13).eq.'NNPDF20intqed') then
            hasphoton = .true.
            nfl = nfl + 1   
            read(20,*) line,line            
            exit
         endif
         if (line(1:13).eq.'NNPDF20int') then
            hasphoton = .false.            
            read(20,*) line,line            
            exit
         endif
      enddo
*
      read(20,*) nx
      do ix=1,nx
         read(20,*) xgrid(ix)
         logxgrid(ix) = dlog(xgrid(ix))
      enddo
*
      read(20,*) nq2
      read(20,*) line
      do iq=1,nq2
         read(20,*) q2grid(iq)
         logq2grid(iq) = dlog(q2grid(iq))
      enddo
*            
      read(20,*) line      
      do imem=0,mem
         do ix=1,nx
            do iq=1,nq2               
               read(20,*) ( pdfgrid(imem,fl,ix,iq), fl=1,nfl,1)            
            enddo
         enddo
      enddo
 
      close(20)
 
      end subroutine
 
      subroutine nnevolvepdf(x,Q,xpdf)
      implicit none
      
      integer i,j,ix,iq2,M,N,ipdf,fmax
      integer minx,maxx,midx
      integer minq,maxq,midq
      double precision x,Q,xpdf(-6:7),Q2
      double precision xmingrid,xch,x2,x1,dy,y
      parameter(m=4, n=2)
      parameter(xmingrid=1d-7, xch=1d-1) 
 
      integer nmax,mmax 
      parameter(nmax=1e3,mmax=1e3)            
      integer ix1a(mmax), ix2a(nmax)
      double precision x1a(mmax), x2a(nmax)
      double precision ya(mmax,nmax)
 
      integer nfl,nx,nq2,mem,rep
      double precision alphas
      double precision xgrid(100),logxgrid(100)
      double precision q2grid(50),logq2grid(50)
      double precision pdfgrid(0:100,14,100,50)
      logical hasphoton
      common /nnpdf/nfl,nx,nq2,mem,rep,hasphoton,alphas,xgrid,logxgrid,
     1     q2grid,logq2grid,pdfgrid
      
 
      q2 = q*q
*     check bounds
      if (x.lt.xmingrid.or.x.lt.xgrid(1).or.x.gt.xgrid(nx)) then
         write(6,*) "Parton interpolation: x out of range -- freezed"
         if (x.lt.xgrid(1)) x = xgrid(1)
         if (x.lt.xmingrid) x = xmingrid
         if (x.gt.xgrid(nx))x = xgrid(nx)
      endif
      if (q2.lt.q2grid(1).or.q2.gt.q2grid(nq2)) then
         write(6,*) "Parton interpolation: Q2 out of range -- freezed"
         if (q2.lt.q2grid(1)) q2 = q2grid(1)
         if (q2.gt.q2grid(nq2)) q2 = q2grid(nq2)
      endif
 
      minx = 1
      maxx = nx+1
 10   continue
      midx = (minx+maxx)/2
      if (x.lt.xgrid(midx)) then
         maxx=midx
      else
         minx=midx
      endif
      if ((maxx-minx).gt.1) go to 10
      ix = minx
 
      minq = 1
      maxq = nq2+1
 20   continue
      midq = (minq+maxq)/2
      if (q2.lt.q2grid(midq)) then
         maxq=midq
      else
         minq=midq
      endif 
      if ((maxq-minq).gt.1) go to 20
      iq2 = minq
 
*     Assign grid for interpolation. M, N -> order of polyN interpolation      
      do i=1,m
         if(ix.ge.m/2.and.ix.le.(nx-m/2)) ix1a(i) = ix - m/2 + i
         if(ix.lt.m/2) ix1a(i) = i
         if(ix.gt.(nx-m/2)) ix1a(i) = (nx - m) + i
         
*     Check grids
         if(ix1a(i).le.0.or.ix1a(i).gt.nx) then
            write(6,*) "Error in grids! "
            write(6,*) "I, IXIA(I) = ",i, ix1a(i)
            call exit(-10)
         endif
      enddo
 
      do j=1,n
         if(iq2.ge.n/2.and.iq2.le.(nq2-n/2)) ix2a(j) = iq2 - n/2 + j
         if(iq2.lt.n/2) ix2a(j) = j
         if(iq2.gt.(nq2-n/2)) ix2a(j) = (nq2 - n) + j
*     Check grids
         if(ix2a(j).le.0.or.ix2a(j).gt.nq2) then
            write(6,*) "Error in grids! "
            write(6,*) "J, IXIA(J) = ",j,ix2a(j)
            call exit(-10)
         endif
      enddo
            
*     Define points where to evaluate interpolation
*     Choose between linear or logarithmic (x,Q2) interpolation
 
      IF(x.LT.xch)THEN
         x1=dlog(x)          
      ELSE
         x1=x
      ENDIF
      x2=dlog(q2)
 
*     initialize output vector
      do i=-6,7 
         xpdf(i) = 0
      enddo
            
      fmax = 6
      if (nfl.eq.14) fmax=7
 
      DO ipdf = -6,fmax,1                 
*     Choose between linear or logarithmic (x,Q2) interpolation        
         DO i=1,m
            IF(x.LT.xch)THEN
               x1a(i)= logxgrid(ix1a(i))
            ELSE
               x1a(i)= xgrid(ix1a(i))
            ENDIF
            DO j=1,n
               x2a(j) = logq2grid(ix2a(j))
               ya(i,j) = pdfgrid(rep,ipdf+7,ix1a(i),ix2a(j))               
            enddo
         enddo
         
!     2D polynomial interpolation
         call lh_polin2(x1a,x2a,ya,m,n,x1,x2,y,dy)
         xpdf(ipdf) = y
      enddo                 
 
      end subroutine
 
      subroutine lh_polin2(x1a,x2a,ya,m,n,x1,x2,y,dy) 
      implicit none 
!                                                                       
      integer m,n,nmax,mmax 
      integer j,k 
      parameter(nmax=1e3,mmax=1e3) 
                                                                        
      real(8) dy,x1,x2,y,x1a(mmax),x2a(nmax),ya(mmax,nmax) 
      real(8) ymtmp(nmax),yntmp(nmax) 
                                                                        
      do j=1,m 
         do k=1,n 
            yntmp(k)=ya(j,k) 
         enddo 
         call lh_polint(x2a,yntmp,n,x2,ymtmp(j),dy) 
      enddo 
      call lh_polint(x1a,ymtmp,m,x1,y,dy) 
!                                                                       
      return 
      END                                           
 
      subroutine lh_polint(xa,ya,n,x,y,dy) 
      implicit none 
!                                                                       
      integer n,NMAX 
!     Largest anticipated value of n                                    
      parameter(nmax=1e3) 
      real(8) dy,x,y,xa(nmax),ya(nmax) 
      integer i,m,ns 
      real(8) den,dif,dift,ho,hp,w,c(nmax),d(nmax) 
      ns=1 
      dif=abs(x-xa(1)) 
      do 11 i=1,n 
         dift=abs(x-xa(i)) 
         if(dift.lt.dif) then 
            ns=i 
            dif=dift 
         endif 
         c(i)=ya(i) 
         d(i)=ya(i) 
   11 enddo 
      y=ya(ns) 
      ns=ns-1 
      do m=1,n-1 
         do i=1,n-m 
            ho=xa(i)-x 
            hp=xa(i+m)-x 
            w=c(i+1)-d(i) 
            den=ho-hp 
            if(den.eq.0) then 
               write(*,*)'failure in polint' 
               stop 
            endif 
            den=w/den 
            d(i)=hp*den 
            c(i)=ho*den 
         enddo 
         if(2*ns.lt.(n-m)) then 
            dy=c(ns+1) 
         else 
            dy=d(ns) 
            ns=ns-1 
         endif 
         y=y+dy 
      enddo 
                                                                        
      return 
      END