!  GeDecay2Body.f90
!
!  Takes care of the two body decay. Parent particle and decay type are
!  specified in old_tree. Kinematical relations:
!  energy daughter 1: e1 = (mp**2 + m1**2 - m2**2)/(2mp)
!  energy daughter 2: e2 = (mp**2 + m2**2 - m1**2)/(2mp)
!  mom.   daughter 1: p1 = sqrt(e1**2 - m1**2)
!
!  Detlef Irmscher, Thomas Ullrich, Uni Heidelberg, September 1993
!  Last update: 18.7.97 ad
!                           implement 1+cos**2 distribution ...
!
subroutine GeDecay2Body (oldindex)
    use genesis_m, self=>GeDecay2Body
    use f90_kind
    implicit none
    logical, save        :: init = .true.
    integer, save        :: angleid
    integer, intent(in)  :: oldindex
    type(tree), pointer  :: old_tree, child1, child2
    real(kind=double)    :: mp, md1, md2, e1, e2, p1, p2
    real(kind=double)    :: costheta, sintheta, cosphi, sinphi, phi
    real                 :: rand(2)
    real, external       :: GeCos2
!
! generate 1+cos**2 angular distribution
!
    if (init) then
        angleid  = workingID
        angleid  = GeFreeHbookId(angleid )
        call hbfun1(angleid,'1+cos(t)**2 2-body',1000,-1.,1.,GeCos2)
        init = .false.
    endif

    ! check index
    if (oldindex < 1 .or. oldindex > lasttree) then
        write (errmess, *) 'funny index: ', oldindex
        call GeErrorMessage ("GeDecay2Body", errmess)
        return
    endif

    ! set up trees
    old_tree => treeList(oldindex)
    ! check indices of children
    if (old_tree%child(1) < 1 .or. old_tree%child(1) > lasttree) then
        write (errmess, *) 'funny index for child 1: ', old_tree%child(1)
        call GeErrorMessage ("GeDecay2Body", errmess)
        return
    endif
    if (old_tree%child(2) < 1 .or. old_tree%child(2) > lasttree) then
        write (errmess, *) 'funny index for child 2: ', old_tree%child(2)
        call GeErrorMessage ("GeDecay2Body", errmess)
        return
    endif

    child1   => treeList(old_tree%child(1))
    child2   => treeList(old_tree%child(2))

    ! check if really 2-body decay
    if (old_tree%decaytype%nBody /= 2) then
                write (errmess, *) 'Incorrect decay type:',         &
                old_tree%decaytype%nBody,                           &
                '-body decay ', old_tree%decaytype%name
        call GeErrorMessage ("GeDecay2Body", errmess, fatal = .true.)
    endif

    ! get the masses
    mp = -1.                    ! has to be smaller than zero for pi0->2 photon
    md1 = old_tree%decaytype%child(1)%mass
    md2 = old_tree%decaytype%child(2)%mass
    do while (mp < md1 + md2)
       mp  = GeMassDistribution(oldindex)
    enddo

    ! calculate CM energies and momenta of children
    e1  = (mp**2 + md1**2 - md2**2)/(2.*mp)
    e2  = mp - e1
    p1  = sqrt((e1 + md1)*(e1 - md1))
    p2  = sqrt((mp**2-(md1+md2)**2)*(mp**2-(md1-md2)**2))/(2.*mp)

    ! orientation in decaying particle rest frame
    call random_number(rand)
    costheta = 2.*rand(1) - 1.
!    costheta = hrndm1 (angleid)

    sintheta = sqrt((1. + costheta)*(1. - costheta))
    phi      = TWOPI * rand(2)

    ! fill 4-vectors of children
    child1%p4%E    = e1
    child1%p4%p%px = p1*sintheta*cos(phi)
    child1%p4%p%py = p1*sintheta*sin(phi)
    child1%p4%p%pz = p1*costheta
    child2%p4%E    = e2
    child2%p4%p    = -child1%p4%p

    ! rotate lepton vectors into coordinate system of virtual photon
    costheta = cos(thetaof(old_tree%p4%p))
    sintheta = sin(thetaof(old_tree%p4%p))
    cosphi   = cos(phiof(old_tree%p4%p))
    sinphi   = sin(phiof(old_tree%p4%p))
    child1%p4%p = GeP3rot(child1%p4%p, costheta, sintheta, cosphi, sinphi)
    child2%p4%p = GeP3rot(child2%p4%p, costheta, sintheta, cosphi, sinphi)

    ! boost secondary particles into the lab
    child1%p4  = GeBoost (child1%p4, old_tree%p4)
    child2%p4  = GeBoost (child2%p4, old_tree%p4)

contains

   function GeMassDistribution (ipart)
       implicit none
       real                        :: GeMassDistribution
       integer, intent(in)         :: ipart
       integer                     :: GeFreeHbookId
       integer, save               :: id_rho, id_omega, id_phi
       type(tree), pointer         :: ptree
       real, external              :: hrndm1, GeGounarisSakurai, GeRndmBreitWigner
       real                        :: width, mass
       logical, save               :: init = .true.
       integer, parameter          :: nbins = 600
       real, parameter             :: upperX = 3.
       real, save                  :: e_mass

       !
       ! initialize the mass distribution for the rho,omega,phi -> e+ e- decay
       !
       if (init) then
           e_mass = particleList(GeParticleIndex("electron"))%mass
           vmass  = particleList(GeParticleIndex("rho"))%mass
	   vwidth = particleList(GeParticleIndex("rho"))%width 
           id_rho = GeFreeHbookId(workingID)
           call hbfun1(id_rho, 'rho mass', nbins, 2*e_mass, upperX, GeGounarisSakurai)
           vmass  = particleList(GeParticleIndex("omega"))%mass
	   vwidth = particleList(GeParticleIndex("omega"))%width 
           id_omega = GeFreeHbookId(workingID)
           call hbfun1(id_omega, 'omega mass', nbins, 2*e_mass, upperX, GeGounarisSakurai)
           vmass  = particleList(GeParticleIndex("phi"))%mass
	   vwidth = particleList(GeParticleIndex("phi"))%width 
           id_phi = GeFreeHbookId(workingID)
           call hbfun1(id_phi, 'phi mass', nbins, 2*e_mass, upperX, GeGounarisSakurai)
	   
	   ! copy histogram to keep original distribution before integration 
	   call HCOPY(id_rho, GeFreeHbookId(originalID),' ')
	   call HCOPY(id_omega, GeFreeHbookId(originalID),' ')
	   call HCOPY(id_phi, GeFreeHbookId(originalID),' ')
	   
           init = .false.
       endif

       ! pointer to particle
       ptree => treeList(ipart)

       ! get width and mass of decay
       width = ptree%parttype%width
       mass  = ptree%parttype%mass

       ! return pure mass if no width given
       if (width == 0) then
           GeMassDistribution = mass
           return
       endif

       !
       ! generate mass, special treatment for the rho, omega, phi
       ! Here we have to check for the particle Id which may cause
       ! trouble if the user changed it.
       !
       if (ptree%parttype%id == 111) then
          GeMassDistribution = hrndm1(id_rho)
       else if (ptree%parttype%id == 221) then
          GeMassDistribution = hrndm1(id_omega)
       else if (ptree%parttype%id == 331) then
          GeMassDistribution = hrndm1(id_phi)
       else
          GeMassDistribution = GeRndmBreitWigner(mass, width)
       endif

   end function GeMassDistribution

end subroutine GeDecay2Body

!
! Mass distribution according to a Breit-Wigner distribution
! Method taken from:
! Review of Particle Properties, Phys. Rev. D 45 part 2 (1992) III.42
!
real function GeRndmBreitWigner(mass, width)
   use genesis_m
   implicit none
   real :: mass, width, u(2), v(2)
   do
      call random_number(u)
      v = 2.*u-1.
      if ( v(1)**2+v(2)**2 <= 1 ) exit
   enddo
   GeRndmBreitWigner = v(1)/v(2)*width/2.+mass
end function GeRndmBreitWigner

!
! Mass distribution for the V -> e+ e- decay (p-wave Breit-Wigner)
! Reference:
! G.J. Gounaris and J.J Sakurai, Phys. Rev. Lett. 21 (1968) 244
! The width and the mass of the vector mesons are passed
! to the function via module genesis_m. This is necessary since
! 'hbfun1' assumes one argument only.
!
real function GeGounarisSakurai(xm)
   use genesis_m
   implicit none
   real          :: xm, eps, corr
   real, save    :: e_mass
   logical, save :: init = .true.
   
   if (init) then
      e_mass = particleList(GeParticleIndex("electron"))%mass
      init = .false.
   endif
   GeGounarisSakurai = 0.
   corr = vwidth*vmass/xm* &
          ((0.25*xm**2-e_mass**2)/(0.25*vmass**2-e_mass**2))**1.5
   eps = (e_mass/xm)**2
   if( 1.-4.*eps < 0.) return
   GeGounarisSakurai = sqrt(1.-4.*eps)*(1.+2.*eps)/ &
              ((vmass**2-xm**2)**2+(vmass*corr)**2)
end function GeGounarisSakurai