diff -Naur POWHEG-hvq/hvqpdfpho.f POWHEG-hvq-new/hvqpdfpho.f
--- POWHEG-hvq/hvqpdfpho.f	2007-05-25 16:02:39.000000000 +0200
+++ POWHEG-hvq-new/hvqpdfpho.f	2007-12-18 09:31:34.000000000 +0100
@@ -6626,8 +6626,8 @@
       Implicit Double Precision (A-H,O-Z)
       Character Flnm(10)*11
       Common
-     > / K719CtqPar2 / Nx, Nt, NfMx
-     > / K719QCDtable /  Alambda, Nfl, Iorder
+     > / K718CtqPar2 / Nx, Nt, NfMx
+     > / K718QCDtable /  Alambda, Nfl, Iorder
       Data (Flnm(I), I=1,10)
      > / 'cteq4m', 'cteq4d', 'cteq4l'
      > , 'cteq4a1', 'cteq4a2', 'cteq4m', 'cteq4a4'
@@ -6694,11 +6694,11 @@
       PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
       PARAMETER (MXPQX = (MXF *2 +2) * MXQ * MXX)
       Common 
-     > / K719CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
-     > / K719CtqPar2 / Nx, Nt, NfMx
-     > / K719XQrange / Qini, Qmax, Xmin
-     > / K719QCDtable /  Alambda, Nfl, Iorder
-     > / K719Masstbl / Amass(6)
+     > / K718CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
+     > / K718CtqPar2 / Nx, Nt, NfMx
+     > / K718XQrange / Qini, Qmax, Xmin
+     > / K718QCDtable /  Alambda, Nfl, Iorder
+     > / K718Masstbl / Amass(6)
       
       Read  (Nu, '(A)') Line     
       Read  (Nu, '(A)') Line
@@ -6746,9 +6746,9 @@
       PARAMETER (M= 2, M1 = M + 1)
 C
       Common 
-     > / K719CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
-     > / K719CtqPar2 / Nx, Nt, NfMx
-     > / K719XQrange / Qini, Qmax, Xmin
+     > / K718CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
+     > / K718CtqPar2 / Nx, Nt, NfMx
+     > / K718XQrange / Qini, Qmax, Xmin
       Dimension Fq(M1), Df(M1)
       data ixrange/0/
       data iqmnrng/0/
diff -Naur POWHEG-hvq/lhapdfif.f POWHEG-hvq-new/lhapdfif.f
--- POWHEG-hvq/lhapdfif.f	1970-01-01 01:00:00.000000000 +0100
+++ POWHEG-hvq-new/lhapdfif.f	2007-12-18 09:31:34.000000000 +0100
@@ -0,0 +1,97 @@
+
+      subroutine genericpdfset(ndns)
+c wrap for pdfset; avoids subsequent
+c calls to pdfset (you never know)
+      integer ndns
+      character * 20 parm(20)
+      real * 8 val(20)
+      integer ndns0
+      data ndns0/-1/
+      save ndns0
+      if(ndns.ne.ndns0) then
+         ndns0=ndns
+         parm(1)='DEFAULT'
+         val(1)=ndns
+         call pdfset(parm,val)
+      endif
+      end
+
+      subroutine genericpdf(ndns,ih,xmu2,x,fx)
+c Interface to lhapdf package.
+      implicit none
+      integer ndns,ih
+      real * 8 xmu2,x,fx(-6:6)
+      integer j
+      real * 8 tmp
+      call genericpdfset(ndns)
+      call pftopdg(x,sqrt(xmu2),fx)
+c pftopdg returns density times x
+      do j=-6,6
+         fx(j)=fx(j)/x
+      enddo
+c 1 is proton, -1 is antiproton, 3 is pi+, -3 is pi-
+      if(ih.eq.1) then
+         return
+      elseif(ih.eq.-1) then
+         do j=1,6
+            tmp=fx(j)
+            fx(j)=fx(-j)
+            fx(-j)=tmp
+         enddo
+      elseif(ih.eq.3) then
+         tmp=fx(1)
+         fx(1)=fx(-1)
+         fx(-1)=tmp
+      elseif(ih.eq.-3) then
+         tmp=fx(2)
+         fx(2)=fx(-2)
+         fx(-2)=tmp
+      elseif(ih.eq.0) then
+c 0 is deuteron
+        fx(1)  = 0.5 * ( fx(1)+fx(2) )
+        fx(-1) = 0.5 * ( fx(-1)+fx(-2) )
+        fx(2)  = fx(1)
+        fx(-2) = fx(-1)         
+      elseif(ih.eq.4) then
+c photon pdf
+         continue
+      else
+         write(*,*) ' genericpdf: unimplemented hadron type ',ih
+         stop
+      endif
+c Bug fixes for version 5.3 of lhapdf
+      if(ndns.eq.363) then
+         do j=1,6
+            fx(j)=fx(j)/2
+            fx(-j)=fx(-j)/2
+         enddo
+      endif
+      
+      end
+
+      subroutine genericpdfpar(ndns,ih,xlam,scheme,iret)
+      implicit none
+      integer ndns,ih
+      real * 8 xlam
+      character * 2 scheme
+      integer iret
+      real * 8 QCDL4,QCDL5
+      COMMON/W50512/QCDL4,QCDL5
+      call genericpdfset(ndns)
+      xlam=qcdl5
+      scheme='MS'
+      iret=0
+c bug fix: lambda=0 CTEQ65 and CTEQ65c sets.
+      if((ndns.ge.10350.and.ndns.le.10390).or.
+     #   (ndns.ge.10450.and.ndns.le.10456)) then
+         if(qcdl5.lt.0.001) then
+            QCDL5 = 0.226d0
+            QCDL4 = 0.326d0
+         endif
+      endif
+      end
+
+      function whichpdfpk()
+      character * 3 whichpdfpk
+      whichpdfpk='lha'
+      end
diff -Naur POWHEG-hvq/Makefile POWHEG-hvq-new/Makefile
--- POWHEG-hvq/Makefile	2007-06-15 12:34:40.000000000 +0200
+++ POWHEG-hvq-new/Makefile	2007-12-18 09:31:34.000000000 +0100
@@ -41,31 +41,39 @@
 %.o: %.c
 	$(CC) $(DEBUG) -c $^
 
+# for mlm pdf
+PDFPACK=mlmpdfif.o hvqpdfpho.o
+
+# for lha pdf uncomment the two following lines, and set
+# the appropriate library location
+#PDFPACK=lhapdfif.o
+#LIBS=-L/home/nason/Pheno/PDFpacks/lib -lLHAPDF -static
+
 # target to generate hard event files
-main-lhef:  main-lhef.o lhefwrite.o POWHEG-hvq.o cernroutines.o powheginput.o hvqpdfpho.o $(SYSOBJ)
-	$(FF) $^ -o $@ $(LIBS) $(CERNLIBS)
+main-lhef:  main-lhef.o lhefwrite.o POWHEG-hvq.o cernroutines.o powheginput.o $(PDFPACK) $(SYSOBJ)
+	$(FF) $^ -o $@ $(LIBS)
 
 # target to generate hard event files in MC@NLO format
-main-mcatnlofl:  main-mcatnlofl.o mcatnlofile.o POWHEG-hvq.o cernroutines.o powheginput.o hvqpdfpho.o $(SYSOBJ)
-	$(FF) $^ -o $@ $(LIBS) $(CERNLIBS)
+main-mcatnlofl:  main-mcatnlofl.o mcatnlofile.o POWHEG-hvq.o cernroutines.o powheginput.o $(PDFPACK) $(SYSOBJ)
+	$(FF) $^ -o $@ $(LIBS)
 
 # target to link to herwig, + generic analysis
-main-HERWIG: main-HERWIG.o herwig6510.o analize-hvq.o mbook.o POWHEG-hvq.o cernroutines.o powheginput.o hvqpdfpho.o $(SYSOBJ)
-	$(FF) $^ -o $@ $(LIBS) $(CERNLIBS)
+main-HERWIG: main-HERWIG.o herwig6510.o analize-hvq.o mbook.o POWHEG-hvq.o cernroutines.o powheginput.o $(PDFPACK) $(SYSOBJ)
+	$(FF) $^ -o $@ $(LIBS)
 
 # target to link to pythia, + generic analysis
-main-PYTHIA: main-PYTHIA.o pythia6326.o analize-hvq.o mbook.o POWHEG-hvq.o cernroutines.o powheginput.o hvqpdfpho.o $(SYSOBJ)
-	$(FF) $^ -o $@ $(LIBS) $(CERNLIBS)
+main-PYTHIA: main-PYTHIA.o pythia6326.o analize-hvq.o mbook.o POWHEG-hvq.o cernroutines.o powheginput.o $(PDFPACK) $(SYSOBJ)
+	$(FF) $^ -o $@ $(LIBS)
 
 # target to read event file, shower them with HERWIG, + generic analysis
 main-HERWIG-lhef: main-HERWIG-lhef.o lhefread.o herwig6510.o analize-hvq.o mbook.o utils.o \
         cernroutines.o $(SYSOBJ)
-	$(FF) $^ -o $@ $(LIBS) $(CERNLIBS)
+	$(FF) $^ -o $@
 
 # target to read event file, shower them with PYTHIA, + generic analysis
 main-PYTHIA-lhef: main-PYTHIA-lhef.o lhefread.o pythia6326.o analize-hvq.o mbook.o utils.o \
         cernroutines.o $(SYSOBJ)
-	$(FF) $^ -o $@ $(LIBS) $(CERNLIBS)
+	$(FF) $^ -o $@
 
 
 all: main-lhef main-mcatnlofl main-HERWIG main-PYTHIA main-HERWIG-lhef \
diff -Naur POWHEG-hvq/mcatnlofile.f POWHEG-hvq-new/mcatnlofile.f
--- POWHEG-hvq/mcatnlofile.f	2007-06-07 11:44:13.000000000 +0200
+++ POWHEG-hvq-new/mcatnlofile.f	2007-12-18 09:31:34.000000000 +0100
@@ -55,7 +55,7 @@
          write(iun,814)sqrt(xm2),'--> M_b'
       elseif(ihwproc.eq.-1706+iue) then
          itopdecmode=powheginput('#topdecaymode') 
-         if(itopdecmode.ne.0.and.itopdecaymode.ne.-1000000) then
+         if(itopdecmode.ne.0.and.itopdecmode.ne.-1000000) then
             il1hw=1
             il2hw=1
             twidth=powheginput('tdec/twidth')
diff -Naur POWHEG-hvq/mlmpdfif.f POWHEG-hvq-new/mlmpdfif.f
--- POWHEG-hvq/mlmpdfif.f	1970-01-01 01:00:00.000000000 +0100
+++ POWHEG-hvq-new/mlmpdfif.f	2007-12-18 09:31:34.000000000 +0100
@@ -0,0 +1,34 @@
+      subroutine genericpdf(ndns,ih,xmu2,x,fx)
+c Interface to mlmpdf package.
+      implicit none
+      integer ndns,ih
+      real * 8 xmu2,x,fx(-6:6)
+      real * 4 sxmu2,sx,sfx(-5:5)
+      integer j
+      sx=x
+      sxmu2=xmu2
+      call mlmpdf(ndns,ih,sxmu2,sx,sfx,5)
+      do j=-5,5
+         fx(j)=sfx(j)
+      enddo
+      fx(1)=sfx(2)
+      fx(-1)=sfx(-2)
+      fx(2)=sfx(1)
+      fx(-2)=sfx(-1)
+      fx(6)=0
+      fx(-6)=0
+      end
+
+      subroutine genericpdfpar(ndns,ih,xlam,scheme,iret)
+      implicit none
+      integer ndns,ih
+      real * 8 xlam
+      character * 2 scheme
+      integer iret
+      call pdfpar(ndns,ih,xlam,scheme,iret)
+      end
+
+      function whichpdfpk()
+      character * 3 whichpdfpk
+      whichpdfpk='mlm'
+      end
diff -Naur POWHEG-hvq/POWHEG-hvq.f POWHEG-hvq-new/POWHEG-hvq.f
--- POWHEG-hvq/POWHEG-hvq.f	2007-12-18 09:31:36.000000000 +0100
+++ POWHEG-hvq-new/POWHEG-hvq.f	2007-12-18 09:33:22.000000000 +0100
@@ -425,7 +425,7 @@
       x1b=sqrt(m2qq/sh)*exp(yqq)
       x2b=m2qq/sh/x1b
 
-      call getxmin(y,xmin)
+      call pwggetxmin(y,xmin)
       x=xmin+xscaled*(1-xmin)
 
       end
@@ -518,7 +518,7 @@
       enddo
       end
 
-      subroutine getxmin(y,xmin)
+      subroutine pwggetxmin(y,xmin)
       implicit none
       real * 8 y,xmin
       real * 8 x1b,x2b
@@ -642,9 +642,9 @@
       real * 8 r0(-5:5),rp(-5:5),rm(-5:5)
       integer j,kb
       include 'fixvar-hvq.h'
-      call getxmin(y,xmin)
-      call getxmin(1d0,xminp)
-      call getxmin(-1d0,xminm)
+      call pwggetxmin(y,xmin)
+      call pwggetxmin(1d0,xminp)
+      call pwggetxmin(-1d0,xminm)
       x=xmin + xscaled*(1-xmin)
       xp=xminp + xscaled*(1-xminp)
       xm=xminm + xscaled*(1-xminm)
@@ -711,9 +711,9 @@
 c returns the flavour of the underlying born: 0 for gg, j for q qbar
       ijfl(j)=(kb-1)*j
 c
-      call getxmin(y,xmin)
-      call getxmin(1d0,xminp)
-      call getxmin(-1d0,xminm)
+      call pwggetxmin(y,xmin)
+      call pwggetxmin(1d0,xminp)
+      call pwggetxmin(-1d0,xminm)
       x=xmin + xscaled*(1-xmin)
       xp=xminp + xscaled*(1-xminp)
       xm=xminm + xscaled*(1-xminm)
@@ -948,6 +948,8 @@
       real * 8 xmq
       real * 8 powheginput
       external powheginput
+      character * 3 whichpdfpk
+      external whichpdfpk
 c store Les Houches process information
       iupper=1
 c Set default value for iupper, ignore alternative small
@@ -962,8 +964,16 @@
 c     
       ih1=powheginput('ih1')
       ih2=powheginput('ih2')
-      ndns1=powheginput('ndns1')
-      ndns2=powheginput('ndns2')
+      if(whichpdfpk().eq.'lha') then
+         ndns1=powheginput('lhans1')
+         ndns2=powheginput('lhans2')
+      elseif(whichpdfpk().eq.'mlm') then
+         ndns1=powheginput('ndns1')
+         ndns2=powheginput('ndns2')
+      else
+         write(*,*) ' unimplemented pdf package',whichpdfpk()
+         stop
+      endif
       ebmup(1)=powheginput('ebeam1')
       ebmup(2)=powheginput('ebeam2')
       xmq=powheginput('qmass')
@@ -2598,8 +2608,15 @@
          ans=vtot
          err=etot
       else
-         ans=(ans/err+vtot/etot)/(1/err+1/etot)
-         err=1/sqrt(1/err**2+1/etot**2)
+c protect against zero error!
+c         ans=(ans/err+vtot/etot)/(1/err+1/etot)
+         ans=(ans*etot+vtot*err)/(err+etot)
+         if(err.ne.0.and.etot.ne.0) then
+            err=1/sqrt(1/err**2+1/etot**2)
+         else
+c if one is zero, pick the non-zero one.
+            err=err+etot
+         endif
       endif
       goto 10
       end
@@ -6693,7 +6710,7 @@
       implicit none
       character * 2 prc
       common/process/prc
-      real*4 fh1x1(-5:5),fh2x2(-5:5)
+      real*8 fh1x1(-6:6),fh2x2(-6:6)
       real * 8 x1,x2,sf
       dimension sf(-5:5)
       include 'fixvar-hvq.h'
@@ -6703,8 +6720,8 @@
       xmuf2h1=xmufct2
       xmuf2h2=xmufct2
 c
-      call mlmpdf0(ndns1,ih1,sngl(xmuf2h1),sngl(x1),fh1x1,5)
-      call mlmpdf0(ndns2,ih2,sngl(xmuf2h2),sngl(x2),fh2x2,5)
+      call genericpdf0(ndns1,ih1,xmuf2h1,x1,fh1x1)
+      call genericpdf0(ndns2,ih2,xmuf2h2,x2,fh2x2)
 c
       do j=-5,5
          sf(j)=0
@@ -6751,31 +6768,31 @@
 
 
 
-c Front end to mlmpdf; it stores the arguments and return values of
-c the nrec most recent calls to mlmpdf. When invoked it looks in the
+c Front end to genericpdf; it stores the arguments and return values of
+c the nrec most recent calls to genericpdf. When invoked it looks in the
 c stored calls; if a match its found, its return value is used.
 c In this framework it is found that nrec=8 would be enough.
 c This provides a remarkable increase in spead (better than a factor of 3)
 c when cteq6 pdf are used.
-      subroutine mlmpdf0(ns,ih,xmu,x,fx,nf)
+      subroutine genericpdf0(ns,ih,xmu2,x,fx)
       implicit none
-      integer ns,ih,nf
-      real * 4 xmu,x,fx(-nf:nf)
+      integer ns,ih
+      real * 8 xmu2,x,fx(-6:6)
       integer nrec
       parameter (nrec=10)
-      real * 4 oxmu(nrec),ox(nrec),ofx(-5:5,nrec)
+      real * 8 oxmu2(nrec),ox(nrec),ofx(-6:6,nrec)
       integer ons(nrec),oih(nrec)
       integer irec
-      save oxmu,ox,ofx,ons,oih,irec
+      save oxmu2,ox,ofx,ons,oih,irec
 c set to impossible values to begin with
       data ox/nrec*-1d0/
       data irec/0/
       integer j,k
       do j=1,nrec
          if(x.eq.ox(j)) then
-            if(xmu.eq.oxmu(j)) then
+            if(xmu2.eq.oxmu2(j)) then
                if(ns.eq.ons(j).and.ih.eq.oih(j)) then
-                  do k=-nf,nf
+                  do k=-6,6
                      fx(k)=ofx(k,j)
                   enddo
                   return
@@ -6787,19 +6804,19 @@
       if(irec.gt.nrec) irec=1
       ons(irec)=ns
       oih(irec)=ih
-      oxmu(irec)=xmu
+      oxmu2(irec)=xmu2
       ox(irec)=x
-      call mlmpdf(ns,ih,xmu,x,ofx(-5,irec),5)
+      call genericpdf(ns,ih,xmu2,x,ofx(-6,irec))
 c Flavour thresholds:
-      if(xmu.lt.10) then
+      if(xmu2.lt.25) then
          ofx(5,irec)=0
          ofx(-5,irec)=0
       endif
-      if(xmu.lt.3) then
+      if(xmu2.lt.2.25) then
          ofx(4,irec)=0
          ofx(-4,irec)=0
       endif
-      do k=-nf,nf
+      do k=-6,6
          fx(k)=ofx(k,irec)
       enddo
       end
@@ -6825,12 +6842,12 @@
       common/nl/nlhvq
       real * 8 xlam3,xlam4,xlam5,xlam6
       integer iret
-      call pdfpar(ndns1,ih1,xlam,scheme1,iret)
+      call genericpdfpar(ndns1,ih1,xlam,scheme1,iret)
       if(iret.ne.0) then
          write(*,*) ' error for pdf 1, iret=',iret
          stop
       endif
-      call pdfpar(ndns2,ih2,xlam,scheme2,iret)
+      call genericpdfpar(ndns2,ih2,xlam,scheme2,iret)
       if(iret.ne.0) then
          write(*,*) ' error for pdf 1, iret=',iret
          stop
@@ -7080,7 +7097,7 @@
          call genjflborn
          y=2*random()-1
          th2=2*pi*random()
-         call getxmin(y,xmin)
+         call pwggetxmin(y,xmin)
          xscaled=random()
          x=(1-xmin)*xscaled+xmin
          call maxrat(xnorm)
@@ -7707,19 +7724,17 @@
 c the top and anti-top invariant masses are off the pole mass
       implicit none
       include 'LesHouches.h'
+      include 'strfun.h'
       real * 8 pdfcorr,corfac
       real * 8 mufc2,x(2),xc(2),r
-      real * 4 fx(-5:5)
+      real * 8 fx(-6:6)
       real * 8 powheginput
       external powheginput
-      integer ini,ndns(2),ipart(2),ih(2),k
-      data ini/0/
-      save ini,ndns
-      if(ini.eq.0) then
-         ndns(1)=powheginput('ndns1')
-         ndns(2)=powheginput('ndns2')
-         ini=1
-      endif
+      integer ndns(2),ipart(2),ih(2),k
+      ndns(1)=ndns1
+      ndns(2)=ndns2
+      ih(1)=ih1
+      ih(2)=ih2
 c Factorization scale (relay on 5 being the light parton)
       if(nup.eq.4) then
          mufc2=5d0
@@ -7745,24 +7760,20 @@
             stop
          endif
       enddo
-c parton types, from pdg to mlmpdf
+c parton types, from pdg to generic pdf
       do k=1,2
          if(idup(k).eq.21) then
             ipart(k)=0
-         elseif(abs(idup(k)).eq.1) then
-            ipart(k)=sign(2,idup(k))
-         elseif(abs(idup(1)).eq.2) then
-            ipart(k)=sign(1,idup(k))
          else
             ipart(k)=idup(k)
          endif
       enddo
       r=1
       do k=1,2
-         call mlmpdf(ndns(k),ih(k),sngl(mufc2),sngl(x(k)),fx,5)
+         call genericpdf(ndns(k),ih(k),mufc2,x(k),fx)
          if(fx(ipart(k)).ne.0) then
             r=r/fx(ipart(k))
-            call mlmpdf(ndns(k),ih(k),sngl(mufc2),sngl(xc(k)),fx,5)
+            call genericpdf(ndns(k),ih(k),mufc2,xc(k),fx)
             r=r*fx(ipart(k))
          endif
       enddo