c sig_factorized_2007.f
      
* p(e,e'K+)Lambda and p(e,e'K+)Sigma cross sections 
* via factorized model described in Tanja Horn's December 2007 report
* model is based on Q^2, W dependence of Bebek and Brauel data, using
* Koltenuk fits which are updated to improve epsilon, etc. dependence

* Cross section (at theta_cm=0, i.e. t=tmin) is F(Q^2)*G(W)
* The t-dependence at fixed Q^2,W is h(t)=D*exp(-b*t),
*     so the cross section is F(Q^2)*G(W)*(h(t)/h(tmin))=F(Q^2)*G(W)*h(t-tmin)

c gh 21.06.14

      implicit none

      real*8 mkg,mpg,mLg,mSg,mpig,mng,pi,hbarc
      real*8 mrec,epsilon
      real*8 q2g,wg,nug,pgamg,betacm,gammacm
      real*8 pcmLg,pcmSg,ekcmLg,ekcmSg,pkaoncm
      real*8 pgamcmg,nuparcmg,tminLg,tminSg,tmin,tt
      real*8 A,X,factQ2,factW,D,hb,factt,jac
      real*8 dsig,dsigT,dsigL,ratio
      
      integer itt,ifs

c stuff for linux_suppl.inc
      real*8 sind, cosd, tand
      real*8 asind, acosd, atand, atan2d
      external sind, cosd, tand
      external asind, acosd, atand, atan2d

      mkg=0.49368
      mpg=0.93827
      mLg=1.11568
      mSg=1.19264
      mpig=0.13957
      mng=0.93957
      
      pi=3.14159
      hbarc=0.197  ! (GeV-fm)
      epsilon=0.5  ! lazy and hard code it for now

 90   write(6,91)
 91   format(' 1:Lambda or 2:Sigma final state?')
      read(5,*)ifs
      if (ifs.eq.1) then
         mrec=mLg
      elseif (ifs.eq.2) then
         mrec=mSg
      else
         goto 90
      endif

      write(6,100)
 100  format(' Enter Q^2, W in GeV')
      read(5,*)q2g,wg

      nug = (wg**2 + q2g - mpg**2)/(2.*mpg)
      pgamg = sqrt( nug**2 + q2g)

c find speed of virtual photon+proton c.m. frame
      betacm  = pgamg/(nug+mpg)
      gammacm = (nug+mpg)/wg

      pcmLg = sqrt( (wg**2 + mLg**2 - mkg**2)**2 - 4.*(wg*mLg)**2 )
     1        /(2.*wg)
      ekcmLg= sqrt(pcmLg**2 + mkg**2)

      pcmSg  = sqrt( (wg**2 + mSg**2 - mkg**2)**2 - 4.*(wg*mSg)**2 )
     1        /(2.*wg)
      ekcmSg= sqrt(pcmSg**2 + mkg**2)

      if (ifs.eq.1) then
         pkaoncm = pcmLg
      else
         pkaoncm = pcmSg
      endif

! calculate t_min
      pgamcmg = (pgamg-betacm*nug)*gammacm
      nuparcmg= nug/gammacm-betacm*pgamcmg                     !p.59 of notes

      tminLg   = q2g-mkg**2+2.*(nuparcmg*ekcmLg-pgamcmg*pcmLg) !p.105
      tminSg   = q2g-mkg**2+2.*(nuparcmg*ekcmSg-pgamcmg*pcmSg)

      if (ifs.eq.1) then
         tmin = -tminLg
      else
         tmin = -tminSg
      endif

! the parameterization factorizes the cross section into three parts

! Q2 parameterization
      if (ifs.eq.1) then        ! Lambda
         A=4.50                 ! ub/sr-GeV^4
         X=2.67                 ! GeV^2
      elseif (ifs.eq.2) then    ! Sigma
         A=0.1587                 
         X=0.785                
      endif
      factQ2=A/(q2g+X)**2   ! ub/sr

! W-dependent part scaled by phase space factor from W=2.15, pK*=0.68234
      if (ifs.eq.1) then        ! Lambda

         if (Wg.ne.0) then   !=0 for central event!!!
            factW = 4.1959*pkaoncm/(sqrt(wg**2)*(wg**2-mpg**2))
            if (Wg.lt.2.) then ! resonance region correction
               factW = 0.959*factW + (0.18*1.72**2*0.10**2) /
     >              ( (wg**2-1.72**2)**2 + 1.72**2*0.10**2 )
            endif
         endif

c old W-dependent part scaled by phase space factor from W=2.15, pK*=0.68234
c      factW=(pkaoncm/(w**2-mp**2)/w)/
c     *     ( 0.682352458165202/(2.15**2-mp**2)/2.15)  ! dimensionless

      elseif (ifs.eq.2) then    ! Sigma
         
         if (Wg.ne.0) then   !=0 for central event!!!
            factW = 1.3176*pkaoncm/(sqrt(wg**2)*(wg**2-mpg**2))
         endif
      endif
         
! t-dependent part, actually a fit to 2pi(d2sig/dtdphi)
      if (ifs.eq.1) then     ! Lambda
         D=1.10                 ! dimensionless
         hb=2.10                ! 1/GeV^2
      elseif (ifs.eq.2) then
         D=1.00                 ! dimensionless
         hb=1.00                ! 1/GeV^2
      endif

      write(6,10)
 10   format(/,4x'Q2',7x,'W',6x'-t',9x,'sigma',6x,'T/L')

      do itt=1,30
         tt= tmin - float(itt-1)*0.05

         factt=D*exp(-hb*abs(tt-tmin)) ! dimensionless

! the product of these functions yields dsigma/dOmegaK* in ub/sr
         dsig=factQ2*factW*factt*(2.+epsilon)/3.

! to convert to dsigma/dt in ub/GeV2, need a Jacobian
! Jacobian (p. 106 of logbook) is dt/dOmega*, so actually need 1/J
         jac = pkaoncm*pgamcmg / 3.1415

c Ratio=T/L
c take T to be epsilon indepenent part above

         dsigT = (factQ2*factW*factt*2./3.)/jac
         dsigL = (factQ2*factW*factt/3.)/jac

         ratio = dsigT / dsigL

         write(6,110)q2g,wg,tt,dsig/jac,ratio
 110     format(3f8.3,3x,2f10.4)
      
      enddo
      end


      include 'linux_suppl.inc'