*------------------------------------------------------------------------*
*    This program calculates the process \mu^- e^- -> \mu^- e^-  as      *
*    described in pages 90-93 of ITC.                                    *
*                                                                        *
*    The purpose is to show how to do numerical integrations. We use 5   *
*    methods: Gauss, Vegas (Monte Carlo), Suave (Monte Carlo), Cuhre     *
*    and Divonne (Monte Carlo). The first uses a routine that I wrote.   *
*    The other 4 are described in the preprint hep-ph/0404043 by T. Hahn *
*    and are part of a free package CUBA. The conventions of my package  *
*    are the same for compatibility. Note that Divonne and Cuhre only    *
*    work for ndim >=2 and my Gauss routine only for ndim <=6            *
*                                                                        *
*    All the methods (including gauss) use integration in an hypercube   *
*    of unit length. So the program as to map the intervals into [0,1].  *
*    The limits of the intervals are stored in the arrays xmin and xmax. *
*                                                                        *
*------------------------------------------------------------------------*

      program muemue
      implicit none

********************************************************************
*                                                                  *      
***   Start variables for Integration                              *
*                                                                  *
********************************************************************


*
***   Variables for Gauss
*

      integer npoints,method
      parameter (npoints = 64)
      parameter (method = 1)

*
***   Variables for all the methods
*

	integer ndim, ncomp, nvec, last, seed, mineval, maxeval
	double precision epsrel, epsabs, userdata
	parameter (ndim = 1)
	parameter (ncomp = 2)
	parameter (userdata = 0)
	parameter (nvec = 1)
	parameter (epsrel = 1D-3)
	parameter (epsabs = 1D-12)
	parameter (last = 4)
	parameter (seed = 0)
	parameter (mineval = 0)
	parameter (maxeval = 50000)

	integer nstart, nincrease, nbatch, gridno
	character*(*) statefile	
	parameter (nstart = 1000)
	parameter (nincrease = 500)
	parameter (nbatch = 1000)
	parameter (gridno = 0)
	parameter (statefile = "")

	integer nnew
	double precision flatness
	parameter (nnew = 1000)
	parameter (flatness = 25D0)

	integer key1, key2, key3, maxpass
	double precision border, maxchisq, mindeviation
	integer ngiven, ldxgiven, nextra
	parameter (key1 = 47)
	parameter (key2 = 1)
	parameter (key3 = 1)
	parameter (maxpass = 5)
	parameter (border = 0D0)
	parameter (maxchisq = 10D0)
	parameter (mindeviation = .25D0)
	parameter (ngiven = 0)
	parameter (ldxgiven = ndim)
	parameter (nextra = 0)

	integer key
	parameter (key = 0)


	double precision integral(ncomp), error(ncomp), prob(ncomp)
	integer  verbose, nregions, neval, fail

********************************************************************
*                                                                  *      
***   End of variables for Integration                             *
*                                                                  *
********************************************************************

      integer i,j
      double precision p1L(4),p2L(4),p2Lab,p1Lab
      double precision m1,m2,m3,m4,pi,alfa,factor
      double precision xmin(10),xmax(10)
      double precision xs_exact,cosxmin,dot,rs
      double precision xs_gauss_CM,xs_vegas_CM
      double precision teta_min_Lab,teta_min_CM
      double precision teta_CM,teta_Lab
      double precision a,b,res1,res2

**************  Standard Model Physical Constants ******
*
      INCLUDE 'smpc.h'
*
********************************************************

      common/boundary/xmin,xmax
      common/data/pi,alfa,factor,m1,m2,m3,m4,rs
      
      external integrand_CM,integrand_Lab

      verbose = 0

      m1=SM_MUMASS      
      m2=SM_EMASS

      alfa=1.d0/128.9d0


*
***  1GeV^(-2)= 3.893795 10^(8) pb (Conversion factor)
*

      factor=3.893795d+08
      pi=acos(-1.d0)

*
**    The cross section will be evaluated in the CM frame and Lab frame.
**    We take in the Lab the electron at rest and cycle the muon momentum
**    in the Lab.
*
*
***   The theta angle in the CM frame is the only integration variable, 
***   The integral will blow up if xmin(1)=0.d0.
***
*

*
***  Loop the CM energy
*

      rs=0.d0

*
**    Cut: xmin = 10 degrees  in CM
*
      teta_min_CM=10.d0*(pi/180.d0)   

      do i=1,100


      xmin(1)=teta_min_CM
      xmax(1)=pi
      cosxmin=cos(xmin(1))


*
**    Calculate the CM energy
*
      rs=rs+1.d0

      call TetaLab(rs,m1,m2,m3,m4,teta_min_Lab,teta_min_CM)
      
      if(i.eq.49.AND..FALSE.)then
         teta_CM=pi/180.d0*175.d0
         do j=1,51
            call TetaLab(rs,m1,m2,m3,m4,teta_Lab,teta_CM)
            write(20,*)teta_CM/pi*180.d0,teta_Lab/pi*180.d0
            teta_CM=teta_CM+pi/180.d0/50.*5.
         enddo
      endif


********************************************************************
*                                                                  *      
***   Call methods with Cuba-3.3                                   *
*                                                                  *
********************************************************************

*
**   Modify the parameter method to choose the method of integration. 
*

      if(method.eq.1)then
         call gauss(ndim, ncomp, npoints, integrand_CM, integral, fail)
         if(fail.eq.1)then
            write(*,*)'Gauss method only works for ndim <= 6'
            stop
         endif

         if(fail.eq.2)then
            write(*,*)'Gauss method only works for
     &           npoints=4,6,8,10,16,20,24,32,64,128'
            stop
         endif
         
      endif
      
      if(method.eq.2)then
	call vegas(ndim, ncomp, integrand_CM, userdata, nvec,
     &    epsrel, epsabs, verbose, seed,
     &    mineval, maxeval, nstart, nincrease, nbatch,
     &    gridno, statefile,
     &    neval, fail, integral, error, prob)
      endif


      if(method.eq.3)then
	call suave(ndim, ncomp, integrand_CM, userdata, nvec,
     &    epsrel, epsabs, verbose + last, seed,
     &    mineval, maxeval, nnew, flatness,
     &    statefile,
     &    nregions, neval, fail, integral, error, prob)
      endif
      
      if(method.eq.4)then
	call divonne(ndim, ncomp, integrand_CM, userdata, nvec,
     &    epsrel, epsabs, verbose, seed,
     &    mineval, maxeval, key1, key2, key3, maxpass,
     &    border, maxchisq, mindeviation,
     &    ngiven, ldxgiven, 0, nextra, 0,
     &    statefile,
     &    nregions, neval, fail, integral, error, prob)
      endif
      
      if(method.eq.5)then
	call cuhre(ndim, ncomp, integrand_CM, userdata, nvec,
     &    epsrel, epsabs, verbose + last,
     &    mineval, maxeval, key,
     &    statefile,
     &    nregions, neval, fail, integral, error, prob)
      endif


********************************************************************
*                                                                  *      
***   End of  methods with Cuba-3.3                                *
*                                                                  *
********************************************************************
        

*
***   Write the results boths for the numerical integration, integral(1)
***   that comes from the called subroutine and for the analytical 
***   integration, from the function xs_exact(cos(teta_min)). 
*
      xs_gauss_CM=integral(1)

	call vegas(ndim, ncomp, integrand_CM, userdata, nvec,
     &    epsrel, epsabs, verbose, seed,
     &    mineval, maxeval, nstart, nincrease, nbatch,
     &    gridno, statefile,
     &    neval, fail, integral, error, prob)
 
      xs_vegas_CM=integral(1)

      xmin(1)=teta_min_Lab
      if(m1.gt.m2)then
         xmax(1)=asin(m2/m1)
      else
         xmax(1)=pi
      endif

      a=xmin(1)
      b=xmax(1)

      xmin(1)=a
      xmax(1)=a+(b-a)/(10.d0*(log(rs)+1.d0))
      call gauss(ndim, ncomp, npoints, integrand_Lab, integral, fail)
      res1=integral(1)
      xmin(1)=xmax(1)
      xmax(1)=b
      call gauss(ndim, ncomp, npoints, integrand_Lab, integral, fail)
      res1=res1+integral(1)
      xmin(1)=a
      xmax(1)=b
      call gauss(ndim, ncomp, npoints, integrand_Lab, integral, fail)
      res2=integral(1)

      write(*,99)rs,xs_exact(cosxmin),xs_gauss_CM,xs_vegas_CM,res1,res2


      end do



 98   format(f8.2,4x,3(e11.5,2x))
 99   format(f8.2,2x,5(e11.5,2x))

      end

*
*** The integration has to be mapped from the interval [0,1] into
*** [xmin,xmax] 
*

      subroutine integrand_CM(ndim, xx, ncomp, ff)
      implicit none
      integer i, ndim, ncomp
      double precision xx(*), ff(*)
      double precision x(ndim),xmin(10),xmax(10),Jacob
      double precision dsig_CM

      common/boundary/xmin,xmax

      Jacob=1.d0

      do i=1,ndim
         Jacob=Jacob*(xmax(i)-xmin(i))
         x(i)=xmin(i)+(xmax(i)-xmin(i))*xx(i)
      enddo
      
      ff(1) = dsig_CM(x(1))*Jacob

      return
      end
*
*** The integration has to be mapped from the interval [0,1] into
*** [xmin,xmax] 
*

      subroutine integrand_Lab(ndim, xx, ncomp, ff)
      implicit none
      integer i, ndim, ncomp
      double precision xx(*), ff(*)
      double precision x(ndim),xmin(10),xmax(10),Jacob
      double precision dsig_Lab

      common/boundary/xmin,xmax

      Jacob=1.d0

      do i=1,ndim
         Jacob=Jacob*(xmax(i)-xmin(i))
         x(i)=xmin(i)+(xmax(i)-xmin(i))*xx(i)
      enddo
      
      ff(1) = dsig_Lab(x(1))*Jacob

      return
      end

*
*** This is the function to integrate. See my ITC notes.
*

      double precision function dsig_CM(teta)
      implicit none
      double precision teta,cst
      double precision pi,alfa,factor
      double precision m1,m2,m3,m4
      double precision p1(4),p2(4),p3(4),p4(4)
      double precision rs,dot,ksq

      common/data/pi,alfa,factor,m1,m2,m3,m4,rs

      call CM_Kinematics(m1,m2,m3,m4,rs,p1,p2,p3,p4,teta)

      cst=factor/(64.d0*pi**2*rs**2)*(4.d0*pi*alfa)**2*8.d0
      ksq=dot(p1,p1)+dot(p3,p3)-2.d0*dot(p1,p3)

      dsig_CM=cst/ksq**2*(dot(p3,p4)*dot(p1,p2)+dot(p3,p2)*dot(p1,p4)
     &     -m2**2*dot(p1,p3)-m1**2*dot(p2,p4)+2.d0*m2**2*m1**2)
      dsig_CM=dsig_CM*sin(teta)*2.d0*pi

      return
      end

*
*** This is the function to integrate. See my ITC notes.
*

      double precision function dsig_Lab(tetaL)
      implicit none
      integer flag,opt
      double precision tetaL,cst
      double precision pi,alfa,factor
      double precision m1,m2,m3,m4
      double precision p1(4),p2(4),p3(4),p4(4)
      double precision rs,dot,ksq,p3Lab,aux1,aux2

      common/data/pi,alfa,factor,m1,m2,m3,m4,rs


      cst=factor/(64.d0*pi**2)*(4.d0*pi*alfa)**2*8.d0


      opt=1
      call Lab_Kinematics(m1,m2,m3,m4,rs,p1,p2,p3,p4,tetaL,opt)
      ksq=dot(p1,p1)+dot(p3,p3)-2.d0*dot(p1,p3)
      aux1=cst/ksq**2*(dot(p3,p4)*dot(p1,p2)+dot(p3,p2)*dot(p1,p4)
     &     -m2**2*dot(p1,p3)-m1**2*dot(p2,p4)+2.d0*m2**2*m1**2)
      p3Lab=sqrt(p3(4)**2-m3**2)

      aux1=aux1/m2*p3Lab/p1(3)/
     &     abs(m2+p1(4)-p1(3)*p3(4)/p3Lab*cos(tetaL))
      aux1=aux1*sin(tetaL)*2.d0*pi

      dsig_Lab=aux1

      if(m1.gt.m2)then
      opt=2
      call Lab_Kinematics(m1,m2,m3,m4,rs,p1,p2,p3,p4,tetaL,opt)
      ksq=dot(p1,p1)+dot(p3,p3)-2.d0*dot(p1,p3)
      aux2=cst/ksq**2*(dot(p3,p4)*dot(p1,p2)+dot(p3,p2)*dot(p1,p4)
     &     -m2**2*dot(p1,p3)-m1**2*dot(p2,p4)+2.d0*m2**2*m1**2)
      p3Lab=sqrt(p3(4)**2-m3**2)

      aux2=aux2/m2*p3Lab/p1(3)/
     &     abs(m2+p1(4)-p1(3)*p3(4)/p3Lab*cos(tetaL))
      aux2=aux2*sin(tetaL)*2.d0*pi

      dsig_Lab=dsig_Lab + aux2
      endif

      return
      end

*
***  Exact cross section
*

      double precision function xs_exact(xmax)
      implicit none
      double precision cst
      double precision pi,alfa,factor
      double precision m1,m2,m3,m4
      double precision p1(4),p2(4),p3(4),p4(4)
      double precision rs,dot,pCM,a0,a1,a2,xmax

      common/data/pi,alfa,factor,m1,m2,m3,m4,rs

*
*** CM_Kinematics is called with teta=0.5d0 to use the same subroutine,
*** but the needed variables do not depend on teta.
*
      call CM_Kinematics(m1,m2,m3,m4,rs,p1,p2,p3,p4,0.5d0)
      cst=factor/(64.d0*pi**2*rs**2)*(4.d0*pi*alfa)**2*8.d0
      
      pCM=sqrt(p3(4)**2-m1**2)


      A0= (m1**8 + m2**8 - 2.d0*m2**6*rs**2 + 8.d0*pCM**4*rs**4 
     &     - 2.d0*m2**2*rs**6 + 4.d0*pCM**2*rs**6 + rs**8 - 
     &     2.d0*m1**6*(2.d0*m2**2 + rs**2) + 2.d0*m1**4*(3.d0*m2**4 
     &     + m2**2*rs**2 -  2.d0*pCM**2*rs**2 + rs**4) + 
     &     m2**4*(-4*pCM**2*rs**2 +  2.d0*rs**4) - 
     &     2.d0*m1**2*(2.d0*m2**6 - m2**4*rs**2 + rs**6 - 2.d0*m2**2*
     &     (2.d0*pCM**2*rs**2 + 3.d0*rs**4)))/(32.d0*pCM**4*rs**4)


      A1=(-m1**4 - m2**4 + 2.d0*m2**2*rs**2 + rs**4 
     &   + 2.d0*m1**2*(m2**2 + rs**2))/ (8.d0*pCM**2*rs**2)

      A2=1.d0/4.d0

      xs_exact=cst*(A0*(1.d0+xmax)/2.d0/(1.d0-xmax)
     &     + A1*(1.d0+xmax+2.d0*(1.d0-xmax)*log((1.d0-xmax)/2.d0))/
     &     2.d0/(1.d0-xmax)
     &     + A2*(2.d0*xmax**2-xmax-3.d0-4.d0*(1.d0-xmax)*
     &     log((1.d0-xmax)/2.d0))/2.d0/(xmax-1.d0))

      xs_exact=xs_exact*2.d0*pi

      return
      end


*
***  Useful function to evaluate scalar products.
*

      double precision function dot(p,q)
      implicit none
      double precision p(4),q(4)
      dot=p(4)*q(4)-p(1)*q(1)-p(2)*q(2)-p(3)*q(3)
      return
      end



*
***  Useful function for pCM momenta.
*

      double precision function lambda(x,y,z)
      implicit none
      double precision x,y,z
      lambda=sqrt((x**2-y**2-z**2)**2-4.d0*y**2*z**2)
      return
      end

*
***  Subroutine that evaluates the CM momenta
*

      subroutine CM_Kinematics(m1,m2,m3,m4,rs,p1,p2,p3,p4,teta)
      implicit none
      double precision p1(4),p2(4),p3(4),p4(4)
      double precision rs,beta,gama,E1CM,E2CM,E3CM,E4CM,p1CM,p3CM
      double precision teta
      double precision m1,m2,m3,m4
      double precision dot,lambda

      p1CM=lambda(rs,m1,m2)/2.d0/rs
      p3CM=lambda(rs,m3,m4)/2.d0/rs

      E1CM=(rs**2+m1**2-m2**2)/(2.d0*rs)
      E2CM=(rs**2+m2**2-m1**2)/(2.d0*rs)
      E3CM=(rs**2+m3**2-m4**2)/(2.d0*rs)
      E4CM=(rs**2+m4**2-m3**2)/(2.d0*rs)

      p1(4)=E1CM
      p1(1)=0.d0
      p1(2)=0.d0
      p1(3)=p1CM

      p2(4)=E2CM
      p2(1)=0.d0
      p2(2)=0.d0
      p2(3)=-p1CM
  
      p3(4)=E3CM
      p3(1)=0.d0
      p3(2)=p3CM*sin(teta)
      p3(3)=p3CM*cos(teta)

      p4(4)=E4CM
      p4(1)=0.d0
      p4(2)=-p3CM*sin(teta)
      p4(3)=-p3CM*cos(teta)

      return
      end


*
***  Subroutine that evaluates the Lab (p2 at rest) momenta
*

      subroutine Lab_Kinematics(m1,m2,m3,m4,rs,p1,p2,p3,p4,teta,opt)
      implicit none
      integer flag,opt
      double precision p1(4),p2(4),p3(4),p4(4)
      double precision rs,beta,gamma,E1CM,E2CM,p1CM
      double precision teta
      double precision m1,m2,m3,m4
      double precision dot,lambda,p3Lab
      double precision A,B,C,Ap,Bp,Cp,x1,x2,delta,E3Lab

      p1CM=lambda(rs,m1,m2)/2.d0/rs
      E1CM=(rs**2+m1**2-m2**2)/(2.d0*rs)
      E2CM=(rs**2+m2**2-m1**2)/(2.d0*rs)
      beta=p1CM/E2CM
      gamma=1.d0/sqrt(1.d0-beta**2)

      if(m1.gt.m2)then
         flag=1
      else
         flag=0
      endif

      p1(4)=gamma*(E1CM+beta*p1CM)
      p1(1)=0.d0
      p1(2)=0.d0
      p1(3)=gamma*(p1CM+beta*E1CM)

      p2(4)=gamma*(E2CM-beta*p1CM)
      p2(1)=0.d0
      p2(2)=0.d0
      p2(3)=gamma*(-p1CM+beta*E2CM)


*
**    To find p3Lab use p4=(E1+m2-E3,0,-p3Lab*sin(teta),p1lab-p3Lab*cos(teta)
**    with p4.p4 = m4**2
*

      A=2.d0*(p1(4)+p2(4))
      B=(p1(4)+p2(4))**2-m4**2+m3**2-p1(3)**2
      C=2.d0*p1(3)*cos(teta)

      Ap=A**2-C**2
      Bp=B*C
      Cp=A**2*m3**2-B**2
      delta=Bp**2-Ap*Cp

      x1=(Bp+sqrt(delta))/Ap
      x2=(Bp-sqrt(delta))/Ap
      
      if(flag.eq.1)then
         if(opt.eq.1)then
            p3Lab=x1
            E3Lab=sqrt(p3Lab**2+m3**2)
         else
            p3Lab=x2
            E3Lab=sqrt(p3Lab**2+m3**2)
         endif
      else
         p3Lab=x1
         E3Lab=sqrt(p3Lab**2+m3**2)
      endif


      p3(4)=E3Lab
      p3(1)=0.d0
      p3(2)=p3Lab*sin(teta)
      p3(3)=p3Lab*cos(teta)

      p4(4)=p1(4)+p2(4)-E3Lab
      p4(1)=0.d0
      p4(2)=-p3Lab*sin(teta)
      p4(3)=p1(3)-p3Lab*cos(teta)


      return
      end

      subroutine TetaLab(rs,m1,m2,m3,m4,teta_Lab,teta_CM)
      implicit none
      double precision rs,m1,m2,m3,m4,teta_Lab,teta_CM
      double precision beta,gamma,lambda
      double precision p1CM,E1CM,E2CM,p3CM,E3CM,num,den

      
      p1CM=lambda(rs,m1,m2)/2.d0/rs
      E1CM=(rs**2+m1**2-m2**2)/(2.d0*rs)
      E2CM=(rs**2+m2**2-m1**2)/(2.d0*rs)
      beta=p1CM/E2CM
      gamma=1.d0/sqrt(1.d0-beta**2)
      
      E3CM=(rs**2+m3**2-m4**2)/(2.d0*rs)
      
      num=sin(teta_CM)
      den=gamma*(cos(teta_CM)+E3CM/E2CM)

      if((num/den).ge.0.d0)then
         teta_Lab=atan(num/den)
      else
         teta_Lab=atan(num/den)+acos(-1.d0)
      endif

      return
      end


