C
C   12/22/09  sross section with Lage/Sergey and 3 variable parameters
C
C     f77 laget_eta.f -o laget -L $CERN/2004/lib -lpacklib -lmathlib -lkernlib
C
      SUBROUTINE CRSEC6(egama,thpi0,sumall,t)
C
C      program laget_eta ! dsigma/dt(microbarn/Gev^2)
c	Labels: o-omega contribution,r-rho meson,c-coulomb
      implicit none
      double precision meta, mn, ebeam
      parameter(meta = 0.54775d0, mn = 0.938272d0)
      real dgamma
      real egama,sumall
      double precision thpi0
      double precision s, s0, t,ao,ar,alpha,phi
      double precision ao_neg,ar_neg
c     common/spar_com/ s, s0, t,ao,ar,alpha,phi
      double precision pi,co,cr,cc,w,wo,wr,wc
      double precision wdomeg
      double precision go,gr,fo,fr,fc
      double precision par(3),dwidth
C
	parameter(pi = 3.1416d0)
C                    PDG after 2004
Ct	parameter(dwidth = 0.51d0)
C                    PDG before 2004
	parameter(dwidth = 0.46d0)
C
       integer n
      double precision tmin, theta

      go=0
      gr=0
C
C       !I put here the phase from Cornell, Sergey
Ct      phi=pi/3.d0     ! not correct, Ashot
C
      phi=0.0d0
C
	par(1) = dwidth   ! magn. of Primakoff process
	par(2) = 1.0d0    ! magn. of hadr. amplitude
	par(3) = phi      ! diff. of interf. phase angle
C
CCCCC	ebeam = 11.0d0
C
	ebeam = DBLE(egama)
C
C       do n=0,100
C      theta=n*0.03d0
*      thetal = theta  / sqrt((mn+2.d0*ebeam)/mn)
C
      theta = thpi0
C
      s0 = 1.0d0
      s  = 2.d0*mn*ebeam + mn**2
      
      co=0.063d0     !0.29*sqrt(6.44/137)
      cr=0.0664d0    !0.81*sqrt(0.92/137)
	 
      tmin = meta**4/(4.d0*ebeam**2)
      t = -4.d0*ebeam**2*dsin(0.5d0*theta*pi/180.d0)**2-tmin
c      print*,'t=', t
      ao=0.44d0+0.9d0*t
      ar=0.55d0+0.8d0*t 
c      print*,'ao=', ao, 'ar=',ar

c      if (ao.lt.1e-16 .or. ar.lt.1e-16) then
c         sumall=0.
c         return
c      endif

     
      

c      if (abs(go) .lt. 1e-16 .or. abs(gr) .lt. 1e-16) then
c         sumall=0.
c         return
c      endif


C                           the omega amplitude
      wo=0.
      fo=0.;
      if (ao.gt.0.) then
         go=dgamma(ao)          !Euler Gamma functions
      else if (ao.lt.0) then
         ao_neg=-ao
         go=-pi/(ao_neg*dgamma(ao_neg)*dsin(pi*ao_neg));
      endif

c      if (ao.le.0.) then 
c         print *,'ao',dsin(pi*ao),go
c         endif     
c         if (go.ne.0.) then
      fo= par(2)*20.0d0*co*ebeam**2*s**(ao-1.d0)*dsin(theta*pi/180.d0)*
     &(pi*0.9d0/(go*dsin(pi*ao)))/(1.414d0*meta)
c      print *,go, go*dsin(pi*ao)

      wo=fo**2
   
      wr=0.
      fr=0.
      if (ar.gt.0.) then 
         gr=dgamma(ar)          !Euler Gamma functions
      else if (ar.lt.0.) then
         ar_neg=-ar
         gr=-pi/(ar_neg*dgamma(ar_neg)*dsin(pi*ar_neg));
      endif

c      if (ar.le.0.) then 
c         print *,'ar',dsin(pi*ar),gr
c         endif
c      if (gr.ne.0) then
      fr= par(2)*20.0d0*cr*ebeam**2*s**(ar-1.d0)*dsin(theta*pi/180.d0)*
     &(pi*0.8d0/(gr*dsin(pi*ar)))/(1.414d0*meta)
      wr=fr**2
     
	fc=-dsqrt((8.0d0*par(1))/(137.0d0*meta**3))*ebeam**2*
     &dsin(theta*pi/180.0d0)/(t*50.0d0)
	wc=fc**2
	w=pi*(wo+wr+wc+2.0d0*dcos(pi*(ar-ao))*fo*fr+
     &2*fc*(fo*dcos(pi*ao+par(3))+fr*dcos(pi*ar+par(3))))/ebeam**2

c        print *,'go',go,'gr',gr
c        print *,'wo',wo,'wr',wr,'wc',wc
C
C         !to change to 
c      dsigma/domega drop out pi at front and ebeam**2 at end of 
C      this expression.
C
*       from Ashot, calculation of d(sigma)/d(omega)
C
C 122809	wdomeg = (ebeam*ebeam*w)/pi
C 122809        in this program the cross sections are in MICROBARNs
C 122809        for the GEANT programs we need it in MILIBARNs
C 122809       Therfore we need to divide to 1000.
C
c	wdomeg = ((ebeam*ebeam*w)/pi)/(1.0d3)
c this is for dsig/domega
	wdomeg = (ebeam*ebeam*w)/pi
c        print *,w,wdomeg
C
C
C      
C for dsig/dt 
c      	wdomeg = w
C
        sumall = SNGL(wdomeg)
C
      return
      end