subroutine xsect(Ebeam_,Mtarg_,Mmeson_)
      real Ebeam_,Mtarg_,Mmeson_
      include 'xsect.inc'
      Ebeam=Ebeam_
      Mtarg=Mtarg_
      Mmeson=Mmeson_
      print *, 'xsect initialized with beam energy',Ebeam
      print *, '                       target mass',Mtarg
      print *, '                        meson mass',Mmeson
      print *, '|tmin| is',getQ2(0),' GeV^2'
      if (Mmeson.lt.0.2) then
        GGwidth=GGwidthpi0
      else
        GGwidth=GGwidtheta
      endif
      end

c returns the value of -t (GeV^2) for a given lab angle theta (radians)
      real function getQ2(theta)
      real theta
      include 'xsect.inc'
      double precision Q2
      double precision Q2was
      double precision dtheta
      dtheta=theta
      Q2was=-1
      Q2=0
      niter=0
      do while (abs(Q2was-Q2).gt.1e-12)
        Q2was=Q2
        Emeson=Ebeam-Q2/(2*Mtarg)
        Pmeson=sqrt(Emeson**2-Mmeson**2)
        Q2=2*Ebeam*(Emeson-Pmeson*cos(dtheta))-Mmeson**2
        niter=niter+1
      enddo
c     print *, 'getQ2 converged after',niter,' interations'
      getQ2=Q2
      end

c returns the value of theta (radians) for a given value of -t (GeV^2)
      real function getTheta(Q2)
      real Q2
      include 'xsect.inc'
      double precision costheta
      Emeson=Ebeam-Q2/(2*Mtarg)
      Pmeson=sqrt(Emeson**2-Mmeson**2)
      costheta=Emeson/Pmeson-(Q2+Mmeson**2)/(2*Ebeam*Pmeson)
      getTheta=acos(costheta)
      end

c returns the electromagnetic form factor of the target
      real function getFFem(Q2)
      real Q2
      include 'xsect.inc'
      getFFem=(betaFF**2/(Q2+betaFF**2))**2
      end

c returns the Primakov differential cross section (microbarns/sr)
      real function diffxs(theta)
      real theta
      include 'xsect.inc'
      real Q2,FF
      Q2=getQ2(theta)
      FF=getFFem(Q2)
      diffxs=GGwidth*8*alphaQED*Ztarg**2/Mmeson**3
     +      *(Pmeson/Emeson)**3
     +      *(Ebeam**2/Q2)**2
     +      *(FF*sin(theta))**2 
     +      *hbarc**2*1e4
      end

c returns the missing mass at angle theta, meson energy E
      real function recoilM2(theta,E)
      real theta,E
      include 'xsect.inc'
      double precision p,dtheta
      dtheta=theta
      p=sqrt(E**2-Mmeson**2)
      recoilM2=Mtarg**2+Mmeson**2
     +        +2*Mtarg*(Ebeam-E)
     +        -2*Ebeam*(E-p*cos(dtheta))
      end