      subroutine qcdfunc(xdummy)

      write(6,*)'**** initialize qcd reweighting functions ***'

      return
      end

*.....**************************************
      real function reweignc(x,q2,iop,inpdf)
*.....**************************************
* Function to reweight NC cross section from PDFLIB structure fucntions
* to that of H1 fit. 
* May also be used to reweight MCs produced with different beam
* energies or lepton charges
* INPUT x,q2 - x and q2 kinematic variables
*       iop  =0  reweignc= SIG(H1 fit) / SIG(pdflib) (same lep Q same Ep)
*            =1  reweignc= SIG(H1 fit Q=-1) / SIG(pdflib Q=+1) (same  Ep)
*            =-1 reweignc= SIG(H1 fit Q=+1) / SIG(pdflib Q=-1) (same  Ep)
*            =2  reweignc= SIG(H1 fit Ep=820) / SIG(pdflib Ep=920) (same  Q)
*            =3  reweignc= SIG(H1 fit Ep=920) / SIG(pdflib Ep=820) (same  Q)
*            =4  reweignc= SIG(H1 fit Ep=920 Q=-1) / SIG(pdflib Ep=820 Q=+1)
* where Q is lepton beam charge and Ep is proton beam energy in GeV
* Unless stated otherwise these quanties are taken from the 
* vectors qlep and epin
      implicit none
      real x,y,q2,y1,y2,scms,scms1,scms2,redxsec,redxsecm
      real signc
      integer iop,inpdf
      vector elepin,epin,qlep
      integer icount,iprint

*.....The variable inpdf used to describe input pdf to django
*     2 = mrsh
*     3 = grvlo
*     4 = cteq4m

      reweignc=0.0

      if (inpdf.le.0) then
         inpdf=2
         icount=icount+1
         iprint=1
         if (icount.le.10) then
            iprint=1
         elseif(icount.le.100) then
            if (mod(icount,10).eq.0) iprint=1
         elseif(icount.le.1000) then
            if (mod(icount,100).eq.0) iprint=1
         elseif(icount.le.10000) then
            if (mod(icount,1000).eq.0) iprint=1
         endif
         if (iprint.eq.1) then
            write(6,*) 'REWEIGNC:bad inpdf; set to 2; # errors=',icount
         endif
      endif

      if(x.le.0.0) return
      if(iop.le.1) then
         scms   = 4.0*elepin(1)*epin(1)
         y      = q2/scms/x

         if (y.gt.1.0) return
         if(iop.eq.0) then
            redxsec = signc(x,q2,  1  ,1,0.0,0.0)
            redxsecm= signc(x,q2,inpdf,1,0.0,0.0)
         elseif(iop.eq.1) then 
            redxsec = signc(x,q2,  1  ,1,0.0,-1.0)
            redxsecm= signc(x,q2,inpdf,1,0.0, 1.0)
         elseif(iop.eq.-1) then 
            redxsec = signc(x,q2,  1  ,1,0.0, 1.0)
            redxsecm= signc(x,q2,inpdf,1,0.0,-1.0)
         endif
      elseif(iop.eq.2) then
         scms1   = 4.0*elepin(1)*820.0
         scms2   = 4.0*elepin(1)*920.0
         y1      = q2/scms1/x
         y2      = q2/scms2/x
         if (y1.gt.1.0) return
         if (y2.gt.1.0) return
         redxsec  = signc(x,q2,  1  ,1,scms1,0.0)
         redxsecm = signc(x,q2,inpdf,1,scms2,0.0)
      elseif(iop.eq.3) then
         scms1   = 4.0*elepin(1)*920.0
         scms2   = 4.0*elepin(1)*820.0
         y1      = q2/scms1/x
         y2      = q2/scms2/x
         if (y1.gt.1.0) return
         if (y2.gt.1.0) return
         redxsec  = signc(x,q2,  1  ,1,scms1,0.0)
         redxsecm = signc(x,q2,inpdf,1,scms2,0.0)
      elseif(iop.eq.4) then
         scms1   = 4.0*elepin(1)*920.0
         scms2   = 4.0*elepin(1)*820.0
         y1      = q2/scms1/x
         y2      = q2/scms2/x
         if (y1.gt.1.0) return
         if (y2.gt.1.0) return
         redxsec  = signc(x,q2,  1  ,1,scms1,1.0)
         redxsecm = signc(x,q2,inpdf,1,scms2,-1.0)
      endif
      if (redxsec.gt.1e-8.AND.redxsecm.gt.1e-8) then
         reweignc=redxsec/redxsecm
      else
         reweignc=1.0
      endif

      return
      end

*.....***********************************
      real function signc(x,q2,iset,ired,scms,qin)
*.....***********************************
      vector elepin,epin,qlep,igamma(1)
      parameter (pi=3.14159267,alem=1.0/137.0359895)
      parameter (conv=0.3893834e+09,xmw=80.41,gmuer=1.16807E-05)
      real scms,qin
*     iset:  1=hiq2 fit
*            2=MRSH
*            3=GRV-LO
*            4=CTEQ 4M
*
*     ired:  0= signc is the full diff.l x-section
*            1= signc is the reduced x-section
*
*      scms  0.0/>0.1  ep centre of mass energy squared
*                      calculated from vectors epin and elepin
*                      / used from set vaule 
*      qin (real)   0.0/+-1.0   charge of lepton . Taken from 
*                           vector qlep if set to 0.0
*
*      igamma(1)   0/1  full NC or gamma only
*
      signc=0.0
      if (scms.lt.0.1) then
         scms   = 4.0*elepin(1)*epin(1)
      endif
      if (abs(qin).lt.0.1) then
         qin=qlep(1)
      endif

      y      = q2/scms/x
      yplus  = 1.0+(1.0-y)**2
      yminus = 1.0-(1.0-y)**2
      fac = x*q2*q2/(2.0*pi*alem*alem*yplus*conv)

      if (y.gt.1.0) return

*.....Warn that gamma only exits for iset=1
      if (igamma(1).eq.1.AND.iset.ne.1) then
         write(6,*) 'SIGNC: gamma only not implemented for iset=',iset
         return
      endif

*.....Remind people that they are calling for gamma only!
      if (igamma(1).eq.1) then
         icount=icount+1
         if (mod(icount,50).eq.0) then
            write(6,*) 'SIGNC: call for GAMMA only'
         endif
      endif

      if (iset.eq.1) then
*........For this fit we can make the gamma-exchange only x-section
         if (igamma(1).eq.1) then
            f2  = f2emhiq200(x,q2)
            fl  = flemhiq200(x,q2)
            xf3 = 0.0
         else
            f2  = f2nchiq200(x,q2)
            fl  = flnchiq200(x,q2)
            xf3 = f3nchiq200(x,q2)
         endif
      elseif (iset.eq.2) then
         f2  = f2ncmrsh00(x,q2)
         fl  = 0.0
         xf3 = f3ncmrsh00(x,q2)
      elseif (iset.eq.3) then
         f2  = f2ncgrvl00(x,q2)
         fl  = 0.0
         xf3 = f3ncgrvl00(x,q2)
      elseif (iset.eq.4) then
         f2  = f2ncct4m00(x,q2)
         fl  = 0.0
         xf3 = f3ncct4m00(x,q2)
      endif

*.....Note: xF3 doesnt swap sign for e+/- but only xG3 does!
      if (qin.gt.0.0) then
         redxsec =f2  - (y**2/yplus)*fl  - (yminus/yplus)*xf3
      else
         redxsec =f2  - (y**2/yplus)*fl  + (yminus/yplus)*xf3
      endif
      dsdxdq2  =  redxsec/fac
      
      if (ired.eq.1) then
         signc=redxsec
      else
         signc = dsdxdq2
      endif

      return
      end

*.....********************************************
      real function sigcc(x,q2,iset,ired,scms,qin)
*.....********************************************
      vector elepin,epin,qlep
      parameter (pi=3.14159267,alem=1.0/137.0359895)
      parameter (conv=0.3893834e+09,xmw=80.41,gmuer=1.16807E-05)
      real scms
*     iset:  1=hiq2 fit
*            2=MRSH
*            5=CTEQ-4L

*     ired:  0= sigcc is the full diff.l x-section
*            1= sigcc is the reduced x-section

*      scms  0.0/>0.1  ep centre of mass energy squared
*                      calculated from vectors epin and elepin
*                      / used from set vaule 
*      qin (real)   0.0/+-1.0   charge of lepton . Taken from 
*                           vector qlep if set to 0.0

      sigcc=0.0
      if (scms.lt.0.1) then
         scms   = 4.0*elepin(1)*epin(1)
      endif
      if (abs(qin).lt.0.1) then
         qin=qlep(1)
      endif

      y      = q2/scms/x
      yplus  = 1.0+(1.0-y)**2
      yminus = 1.0-(1.0-y)**2

      if (y.gt.1.0) return

* ER  Note the propagator defined here has an extra factor of 2 so as to
*     match the definitions as used in HECTOR (routine HECSET) for the
*     structure functions.
      wprop=gmuer**2/4d0/pi/x*(xmw**2/(q2+xmw**2))**2

      if (qin.gt.0.0) then
         if (iset.eq.1) then
            w2 = w2cphiq200(x,q2)
            wl = wlcphiq200(x,q2)
            xw3= w3cphiq200(x,q2)
         elseif (iset.eq.2) then
            w2  = w2cpmrsh00(x,q2)
            wl  = 0.0
            xw3 = w3cpmrsh00(x,q2)
         elseif (iset.eq.5) then
            w2  = w2cpct4l00(x,q2)
            wl  = 0.0
            xw3 = w3cpct4l00(x,q2)
         endif
      else
         if (iset.eq.1) then
            w2 = w2cehiq200(x,q2)
            wl = wlcehiq200(x,q2)
            xw3= w3cehiq200(x,q2)
         elseif (iset.eq.2) then
            w2 = w2cemrsh00(x,q2)
            wl = 0.0
            xw3= w3cemrsh00(x,q2)
         elseif (iset.eq.5) then
            w2 = w2cect4l00(x,q2)
            wl = 0.0
            xw3= w3cect4l00(x,q2)
         endif
         
      endif

* ER  Note there is an extra factor of 2 in the definition of the
*     reduced cross-section, this is because the structure functions are
*     defined using the hector prescription where e.g. W2=2d0*(db+sb+u+c)
*     (routine STRUFC). That definition is used to create the
*     parametrisation functions of the S.F. in QCDNUM, and so MUST be
*     adopted here.
      pol = 0.0
      sfterm  = (yplus*W2  - qin*yminus*XW3  - y**2*WL )

      dsdxdq2  = sfterm  * wprop * conv * (1.0+qin*pol)/2.0 
      redxs = dsdxdq2 / (2.0 * wprop * conv ) 

c     redpap is the red.x.sec defined in DESY 99-107
c     The differencce is just the values taken for M_w and G_f
c         gf=1.16639E-5
c         xmw_pap= 80.4
c         wprop_pap = gf**2/2d0/pi/xx*(xmw_pap**2/(qq2+xmw_pap**2))**2
c         redpap = fparton / wprop_pap / conv

      if (ired.eq.1) then
         sigcc=redxs
      else
         sigcc = dsdxdq2
      endif

      return
      end

      real function sigccy(x,y)
      vector epin, elepin

      q2=x*y*4*epin(1)*elepin(1)
      sigccy=sigcc(x,q2,1,1,0.0,0.0)
      return
      end

*.....***********************************
      real function sigy(y,q2,iset,ired,if2)
*.....***********************************

      vector elepin,epin,qlep

*.....plot NC cross-section/red-x-sec/or F2 as func. of (y,q2)
c     if2 = 1 ---> plot F2 only
c           2 ---> plot red-x-sec or diff'l x-sect.
c           3 ---> plot Fl only

      x=q2/(4.0*elepin(1)*epin(1))/y
      if (if2.eq.1) then
         sigy=sigf2(x,q2,iset,ired)
      elseif (if2.eq.0) then
         sigy=signc(x,q2,iset,ired,0.0,0.0)
      elseif (if2.eq.2) then
         sigy=sigfl(x,q2,iset)
      endif
      return
      end

*.....***********************************
      real function sigfl(x,q2,iset)
*.....***********************************
      vector elepin,epin,qlep
      parameter (pi=3.14159267,alem=1.0/137.0359895)
      parameter (conv=0.3893834e+09,xmw=80.41,gmuer=1.16807E-05)

*.....plot Fl from hiq2 fit

*     iset:  1=hiq2 fit
*            2=MRSH
*
      sigfl=0.0
      scms   = 4.0*elepin(1)*epin(1)
      y      = q2/scms/x

      if (y.gt.1.0) return

      if (iset.eq.1) then
         fl  = flnchiq200(x,q2)
      elseif (iset.eq.2) then
         fl  = 0.0
      endif

      sigfl=fl

      return
      end
*.....***********************************
      real function sigf2(x,q2,iset,ired)
*.....***********************************
      vector elepin,epin,qlep
      parameter (pi=3.14159267,alem=1.0/137.0359895)
      parameter (conv=0.3893834e+09,xmw=80.41,gmuer=1.16807E-05)

*.....plot the red-x-sec or full diff'l x-section with ONLY
*     F2 contribution, i.e. fl = xf3 = 0

*     iset:  1=hiq2 fit
*            2=MRSH
*
*     ired:  0= sigf2 is the full diff.l x-section
*            1= sigf2 is the reduced x-section

      sigf2=0.0
      scms   = 4.0*elepin(1)*epin(1)
      y      = q2/scms/x
      yplus  = 1.0+(1.0-y)**2
      yminus = 1.0-(1.0-y)**2
      fac = x*q2*q2/(2.0*pi*alem*alem*yplus*conv)

      if (y.gt.1.0) return

      if (iset.eq.1) then
         f2  = f2nchiq200(x,q2)
      elseif (iset.eq.2) then
         f2  = f2ncmrsh00(x,q2)
      endif
      fl=0.0
      xf3=0.0

*.....Note: xF3 doesnt swap sign for e+/- but only xG3 does!
      if (qlep(1).gt.0.0) then
         redxsec =f2  - (y**2/yplus)*fl  - (yminus/yplus)*xf3
      else
         redxsec =f2  - (y**2/yplus)*fl  + (yminus/yplus)*xf3
      endif
      dsdxdq2  =  redxsec/fac
      
      if (ired.eq.1) then
         sigf2=redxsec
      else
         sigf2 = dsdxdq2
      endif

      return
      end

*.....***********************************
      real function sigflcor(x,q2,iset)
*.....***********************************
      vector elepin,epin,qlep
      parameter (pi=3.14159267,alem=1.0/137.0359895)
      parameter (conv=0.3893834e+09,xmw=80.41,gmuer=1.16807E-05)

*     iset:  1=hiq2 fit
*            2=MRSH
*
*     Fl = sigflcor - redxsec

      sigflcor=0.0
      scms   = 4.0*elepin(1)*epin(1)
      y      = q2/scms/x
      yplus  = 1.0+(1.0-y)**2
      yminus = 1.0-(1.0-y)**2

      if (y.gt.1.0) return

      if (iset.eq.1) then
         f2  = f2nchiq200(x,q2)
         fl  = flnchiq200(x,q2)
         xf3 = f3nchiq200(x,q2)
      elseif (iset.eq.2) then
         f2  = f2ncmrsh00(x,q2)
         fl  = 0.0
         xf3 = f3ncmrsh00(x,q2)
      endif

      sigflcor = (yplus/y**2)*f2 + (yminus/y**2)*xf3
      
      return
      end

*.....***************************
      real function djccwg(x,q2,inpdf)
*.....***************************
      vector elepin,epin,qlep
      parameter (pi=3.14159267,alem=1.0/137.0359895)
      parameter (conv=0.3893834e+09,xmw=80.41,gmuer=1.16807E-05)

*.... Reweight function to reweight MRSH CC x-sec to
*     the CC x-sec from the hiq2 fit

*.....The variable inpdf used to describe input pdf to django
*     2 = mrsh
*     5 = cteq4l

      djccwg=0.0
      scms   = 4.0*elepin(1)*epin(1)
      y      = q2/scms/x
      yplus  = 1.0+(1.0-y)**2
      yminus = 1.0-(1.0-y)**2

      if (y.gt.1.0) return

      redxs  = sigcc(x,q2,1,1,0.0,0.0)
      redxsm = sigcc(x,q2,inpdf,1,0.0,0.0)

      if (redxs.gt.1e-8.AND.redxsm.gt.1e-8) then
         djccwg=redxs/redxsm
      else
         djccwg=1.0
      endif

      return
      end
*.....***************************
      real function wnctocc(x,q2)
*.....***************************

*.....AM warning minumum and maximum x,y, and q2
*.....AM only valid for y<0.975

*     weighting function to reweight NC x-sect to CC x-sect.

      vector elepin,epin,qlep
      parameter (pi=3.14159267,alem=1.0/137.0359895)
      parameter (conv=0.3893834e+09,xmw=80.41,gmuer=1.16807E-05)

      scms   = 4.0*elepin(1)*epin(1)
      if(x.gt.0.99) x=0.99
      if(x.lt.0.0000001) x=0.0000001
      if(q2.lt.0.1) q2=0.1
      if(q2.gt.45000) q2=45000

      y      = q2/scms/x
      yplus  = 1.0+(1.0-y)**2
      yminus = 1.0-(1.0-y)**2

      wnctocc=0.0
      fac = x*q2*q2/(2.0*pi*alem*alem*yplus*conv)

      f2  = f2nchiq200(x,q2)
      fl  = flnchiq200(x,q2)
      xf3 = f3nchiq200(x,q2)

*.....Note: xF3 doesnt swap sign for e+/- but only xG3 does!
      if (qlep(1).gt.0.0) then
         redxsec =f2  - (y**2/yplus)*fl  - (yminus/yplus)*xf3
      else
         redxsec =f2  - (y**2/yplus)*fl  + (yminus/yplus)*xf3
      endif
      signc  =  redxsec/fac


      wprop=gmuer**2*xmw**4*scms/4.0/pi/(q2+xmw**2)**2      
      if (qlep(1).gt.0.0) then
         w2 = w2cphiq200(x,q2)
         wl = wlcphiq200(x,q2)
         xw3= w3cphiq200(x,q2)
      else
         w2 = w2cehiq200(x,q2)
         wl = wlcehiq200(x,q2)
         xw3= w3cehiq200(x,q2)
      endif
      redxs  = (yplus*W2  - yminus*XW3  - y**2*WL ) / 2.0
      sigcc  = redxs  * wprop * conv * y/q2
      if(signc.gt.0.0) then
          wnctocc=sigcc/signc
      else
         wnctocc=1.0
         write(6,*) 'wnctocc out of range  x,q2,y',x,q2,y
      endif
      return
      end

*.....*****************************************************
      real function getsigtot(icc,q2min,ymin,ymax,ptmin,eemin)
*.....*****************************************************
*
*     This is just a wrapper for sigtot - makes it into a function
*
      getsigtot=0.0
      call sigtot(icc,q2min,ymin,ymax,ptmin,eemin,val)
      getsigtot=val
      return
      end

*.....******************************************************
      subroutine sigtot(icc,q2min,ymin,ymax,ptmin,eemin,val)
*.....******************************************************


*.....Calculate total x-section for NC or for CC with different cuts
*     icc = 0/1 for NC/CC
*     sigma is given in range  ymin < y < ymax
*     ptmin cut applied for CC only
*     eemin cut applied for NC only
*     theta max cut implemented but not used for NC

      implicit none

      vector epin(1),elepin(1),qlep(1)
      
      real scms,q2min,q2max,xmin,xmax,ymax
      real dq2,dx,q2width,xwidth
      real xc,q2c,ec,thc,yc,ptc
      real xsec,signc,sigcc,sigtit
      real dsdq2,bxmin,bxmax,ybmax
      integer nn,i,j,icc
      real ptmin,ymin,eemin,val,thmax

      scms = 4.0*elepin(1)*epin(1)

*.....500 bins is a good number, accurate to 1.10^-4
      nn=500
      sigtit = 0.0

      thmax=180.0
      q2max=50000.0
      xmin=0.00001
      xmax=0.99

      dq2=(log(q2max)-log(q2min))/float(nn)
      dx =(log(xmax )-log(xmin ))/float(nn)
      do 10 i=1,nn
         dsdq2=0.0
         bxmin=1.0
         bxmax=-1.0
         do 20 j=1,nn
            q2c = exp(log(q2min)+(float(i)-0.5)*dq2)
            xc  = exp(log(xmin )+(float(j)-0.5)*dx)
            yc  = q2c/(scms*xc)
            ptc = sqrt(q2c*(1.0-yc))
            ec  = elepin(1)*(1.0-yc)+epin(1)*xc*yc
            thc = 2.0*acos(sqrt(q2c/(4.0*elepin(1)*ec)))
            if (yc.lt.ymin)       goto 20
            if (yc.gt.ymax)       goto 20
            if (q2c.lt.q2min)     goto 20
            if (icc.eq.1) then
               if (ptc.lt.ptmin)  goto 20
               xsec=sigcc(xc,q2c,1,0,scms,qlep(1))
            else
               if (thc*57.3.gt.thmax) goto 20
               if (ec.lt.eemin)      goto 20
               xsec=signc(xc,q2c,1,0,scms,qlep(1))
            endif
            sigtit=sigtit+xsec*xc*dx*q2c*dq2  ! jacobian applied for log integral!
            dsdq2=dsdq2+xsec*xc*dx
            bxmin=min(bxmin,xc)
            bxmax=max(bxmax,xc)
 20      continue
         ybmax = q2c/scms/bxmin
 10   continue
      val = sigtit
      return
      end

*.....******************************************************
      real function dsdq2(icc,q2c,ymax,emin,ptmin,scms,qin)
*.....******************************************************

      implicit none
      
*     integration for dsigma/dq2
*     ptmin cut applied for CC only
*     eemin cut applied for NC only
*     ymax cut applied for NC and CC
*     theta max cut implemented but not used for NC

*     these inputs overide vectors epin(1),elepin(1),qlep(1)
*     scms = centre of mass energy
*     qin  = lepton charge; or qin=5 means gamma exchange only

      vector epin(1),elepin(1),qlep(1)

      real scms,sigtot,q2min,q2max,xmin,xmax,ymax,emin,thmin,ptmin
      real dq2,dx,q2width,xwidth,qin,ptc
      
      real xc,q2c,ec,thc,yc
      real xsec,signc,sigcc,sigtit,sigf2
      real bxmin,bxmax,ybmax
      integer nn,i,j,icc

      if (scms.lt.0.1) then
         scms   = 4.0*elepin(1)*epin(1)
      endif
      if (abs(qin).lt.0.1) then
         qin=qlep(1)
      endif

*.....500 bins is a good number, accurate to 1.10^-4
      nn=500
      xmin=q2c/scms
      xmax=0.95

      dx =(log(xmax )-log(xmin ))/float(nn)
      dsdq2=0.0
      bxmin=1.0
      bxmax=-1.0
      do 20 j=1,nn
         xc  = log(xmin )+(float(j)-0.5)*dx   
         xc  = exp(xc)
         yc=q2c/(scms*xc)
         if (yc.gt.ymax)       goto 20
         if (icc.eq.1) then
            ptc  = sqrt( q2c*(1.0-yc) )
            if (ptc.lt.ptmin)  goto 20
            xsec = sigcc(xc,q2c,1,0,scms,qin)
         else
            ec=elepin(1)*(1.0-yc)+epin(1)*xc*yc
c            thc=2.0*acos(sqrt(q2c/(4.0*elepin(1)*ec)))
c            if (thc*57.3.gt.thmax) goto 20
            if (ec.lt.emin)       goto 20
* gamma exchange only if qin=5
            if (qin.ne.5) then
               xsec=signc(xc,q2c,1,0,scms,qin)
            else
               xsec=sigf2(xc,q2c,1,0)
            endif
         endif
         dsdq2=dsdq2+xsec*xc*dx
         bxmin=min(bxmin,xc)
         bxmax=max(bxmax,xc)
 20   continue
      return
      end



*.....*********************************************************
      real function dsdx(icc,xc,q2min,ymax,emin,ptmin,scms,qin)
*.....*********************************************************

*     integration for dsigma/dx
*     ptmin cut applied for CC only
*     eemin cut applied for NC only
*     ymax cut applied for NC and CC

*     these inputs overide vectors epin(1),elepin(1),qlep(1)
*     scms = centre of mass energy
*     qin  = lepton charge

      implicit none
      
      vector epin(1),elepin(1),qlep(1)

      real scms,sigtot,q2min,q2max,xmin,xmax,ymax,emin,thmin,ptmin
      real dq2,dx,q2width,xwidth,qin,ptc
      
      real xc,q2c,ec,thc,yc
      real xsec,signc,sigcc,sigtit
      real bxmin,bxmax,ybmax
      integer nn,i,j,icc

      if (scms.lt.0.1) then
         scms   = 4.0*elepin(1)*epin(1)
      endif
      if (abs(qin).lt.0.1) then
         qin=qlep(1)
      endif

*.....500 bins is a good number, accurate to 1.10^-4
      nn=5000
      q2max=scms*xc*0.9
      q2max=46334.0

      dq2 =(log(q2max)-log(q2min))/float(nn)
      dsdx=0.0
      do 20 j=1,nn
         q2c  = log(q2min )+(float(j)-0.5)*dq2   
         q2c  = exp(q2c)
         yc=q2c/(scms*xc)
         if (yc.gt.ymax)       goto 20

         if (icc.eq.1) then
            ptc  = sqrt( q2c*(1.0-yc) )
            if (ptc.lt.ptmin)  goto 20
            xsec = sigcc(xc,q2c,1,0,scms,qin)
         else
            ec=elepin(1)*(1.0-yc)+epin(1)*xc*yc
c            thc=2.0*acos(sqrt(q2c/(4.0*elepin(1)*ec)))
c            if (thc*57.3.gt.thmax) goto 20
            if (ec.lt.emin)       goto 20
            xsec=signc(xc,q2c,1,0,scms,qin)
         endif

         dsdx=dsdx+xsec*q2c*dq2
 20   continue
*      write(6,*) q2c,xc,yc,dsdx
      return
      end

*.....******************************************************
      real function q2integ(q2c,xmin,xmax)
*.....******************************************************

      implicit none
      
      vector epin(1),elepin(1),qlep(1)

      real scms,sigtot,q2min,q2max,xmin,xmax,ymax,emin,thmin,ptmin
      real dq2,dx,q2width,xwidth,qin,ptc
      real dsdq2,dvalhiq2ms,uvalhiq2ms,glu0hiq2ms
      real xc,q2c,ec,thc,yc
      real xsec,signc,sigcc,sigtit,sigf2
      real bxmin,bxmax,ybmax
      integer nn,i,j,icc

      scms   = 4.0*elepin(1)*epin(1)
      qin=qlep(1)

*.....500 bins is a good number, accurate to 1.10^-4
      nn=500
c      xmin=0.0001
*      xmin=q2c/scms
c      xmax=0.001
*      xmax=0.95
*      ymax=0.99

      dx =(log(xmax )-log(xmin ))/float(nn)
      dsdq2=0.0
      bxmin=1.0
      bxmax=-1.0
      do 20 j=1,nn
         xc  = log(xmin )+(float(j)-0.5)*dx   
         xc  = exp(xc)
         yc=q2c/(scms*xc)
*         if (yc.gt.ymax)       goto 20
         xsec = glu0hiq2ms(xc,q2c)/xc
         dsdq2=dsdq2+xsec*xc*dx
         bxmin=min(bxmin,xc)
         bxmax=max(bxmax,xc)
 20   continue
      q2integ=dsdq2
*      write(*,*) 'q2integ',q2c,dsdq2
      return
      end

*.....******************************************************
      DOUBLE PRECISION FUNCTION DCHSUM(MODE,C,N,X)
*.....******************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION C(0:N)

      IF(MODE .EQ. 1) THEN
       H=X
       F=H
       V=1
      ELSE IF(MODE .EQ. 2) THEN
       H=2*X**2-1
       F=H
       V=1
      ELSE IF(MODE .EQ. 3) THEN
       H=2*X**2-1
       F=1
       V=X
      ELSE IF(MODE .EQ. 4) THEN
       H=2*X-1
       F=H
       V=1
      END IF

      ALFA=H+H
      B1=0
      B2=0
      DO 1 I = N,0,-1
      B0=C(I)+ALFA*B1-B2
      B2=B1
    1 B1=B0
      DCHSUM=V*(B0-F*B2)
      RETURN
      END

*.....****************
      subroutine f2vec
*.....****************

*.....Make reweight jpam vectors for NC and CC
*     Note that for CC the reduced x-section is now made
*     with electroweak params defined as below

*.....Also Note that integrated x-sections are not filled 
*     in the jpam vector

      parameter (pi=3.14159267,alem=1.0/137.0359895)
      parameter (conv=0.3893834e+09,xmw=80.41,gmuer=1.16807E-05)

      vector q2f(800),xf(800),yf(800)
      vector sigera(800)

*.... For NC first
      do i=1,400
         if(yf(i).lt.1.0.and.yf(i).gt.0.0) then
            sigera(i)=signc(xf(i),q2f(i),1,1,0.0,0.0)
         else
            sigera(i)=0.0
         endif
      enddo

*.... For NC ds/dq2  
      do i=401,424
         if(q2f(i).gt.1.0) then
            sigera(i)=dsdq2(0,q2f(i),0.9,0.1,0.1,0.0,0.0)
         else
            sigera(i)=0.0
         endif
      enddo

*.... For NC ds/dx  q2>1000
      do i=430,439
         if(xf(i).gt.1e-7) then
            sigera(i)=dsdx(0,xf(i),1000.0,0.9,0.1,0.1,0.0,0.0)
            write(6,*) 'nc dsdx: q2>1000  ',xf(i),sigera(i)
         else
            sigera(i)=0.0
         endif
      enddo

*.... For NC ds/dx  q2>10000
      do i=450,454
         if(xf(i).gt.1e-7) then
            sigera(i)=dsdx(0,xf(i),10000.0,0.9,0.1,0.1,0.0,0.0)
            write(6,*) 'nc dsdx: q2>10000 ',xf(i),sigera(i)
         else
            sigera(i)=0.0
         endif
      enddo

*.... For NC high y
      do i=456,490
         if(yf(i).lt.1.0.and.yf(i).gt.0.0) then
            sigera(i)=signc(xf(i),q2f(i),1,1,0.0,0.0)
         else
            sigera(i)=0.0
         endif
      enddo

*.... Now do CC cross sections
      do i=501,572
         if(yf(i).lt.1.0.and.yf(i).gt.0.0) then
            sigera(i)=sigcc(xf(i),q2f(i),1,0,0.0,0.0)
         else
            sigera(i)=0.0
         endif
      enddo

*.... Now do CC sigtot
*      print *, 'sigtot'
      call sigtot(1,  200.0,0.0,0.9,0.01,0.01,val1)
      call sigtot(1, 1000.0,0.0,0.9,0.01,0.01,val2)
      call sigtot(1,10000.0,0.0,0.9,0.01,0.01,val3)
      sigera(620)=val1
      sigera(621)=val2
      sigera(622)=val3

*.... For CC ds/dq2
      do i=601,609
         if(q2f(i).gt.1.0) then
            sigera(i)=dsdq2(1,q2f(i),0.9,0.1,0.1,0.0,0.0)
         else
            sigera(i)=0.0
         endif
c         write(6,*) 'f2vec: CC dsdq2 = ',q2f(i),sigera(i)
      enddo

*.... For NC in the xF3 binning
      do i=631,680
         if(yf(i).lt.1.0.and.yf(i).gt.0.0.AND.xf(i).lt.0.9) then
            sigera(i)=signc(xf(i),q2f(i),1,1,0.0,0.0)
         else
            sigera(i)=0.0
         endif
      enddo

      return
      end

*.....****************
      subroutine flfac
*.....****************

      vector flcor(800),xf,yf,q2f

*.... For NC first
      do i=456,490
         if (yf(i).lt.1.0.and.yf(i).gt.0.0) then
            flcor(i)=sigflcor(xf(i),q2f(i),1)
         else
            flcor(i)=0.0
         endif
      enddo
      return
      end




