* this routine is calculate the three body decay M--->m1+m2+m3
c
c     initial particle M has energy E0, Mass M0, and momentum P0 in lab frame 
c     m(3), P(3), theta(3), phi(3) gives the mass, momentum, theta, and phi
c     for three decay particles in lab frame.  
c  
c     input: (1) yamx, the max value of P1_cc*P3_c, P1_cc is momentum of m1 in 
c                CMS of m1+m2, and P3_c is momentum of m3 in CMS of M
c            (2) P0, momentum of decay particle M
c            (3) M0, mass of   decay particle M
c             
c    output: (1) P(3), momentum for m1, m2, m3 in lab frame
c            (2) theta(3), theta angles for m1, m2, m3 in lab frame 
c            (3) phi(3), phi angles for m1, m2, m3 in lab frame
c
c    Important: theta(3) and phi(3) are respecting
c     to the direction of initial direction of decay particle M.
c 
c
c     Sept 12, 2003    by Liping Gan
c
c
c***********************************************************************






      subroutine body3(ymax,P0,M0,m,P,theta,phi)
     
      double precision P0,M0,m(3),P(3),theta(3),phi(3),E(3)
      double precision cos3_c,cos1_cc,m12,V_c,r_c,PI
      double precision P3_c,E0,P1_cc,E1_cc,E3_c,alpha1,alpha2
c  note: E12,P12, theta12,phi_12 are energy, momentum and angle for the
c        center of mass of m1+m2 in lab frame, 
c        the12(2), phi12(2) are the angles for m1 and m2 in lab with 
c       z direction along theta12
      double precision E12,P12, theta12,phi_12,the12(2),phi12(2)
      double precision px_12(2), py_12(2), pz_12(2), px(2), 
     >           py(2), pz(2),p_xy(2)
      double precision xl,xh,ymax,ranf
        real*4 ran_true,x
        

      PARAMETER (PI=3.141592653d0)
	  
          E0=dsqrt(P0**2+M0**2)
      xl=m(1)+m(2)
      xh=M0-m(3)
      
      m12=ranf(xl,xh,ymax,m,M0)


c in the CMS of decay particle M      
        call ranlux(ran_true,1)
      cos3_c=(dble(ran_true)-0.5d0)*2.0d0
      
      P3_c=((M0**2-(m12+m(3))**2)*(M0**2-(m12-m(3))**2))**0.5/2.d0/M0
 
      V_c=P0/E0
      r_c=E0/M0

* calculate in the center mass of m1 and m2
      P1_cc=dsqrt((m12**2-(m(1)+m(2))**2)*
     >      (m12**2-(m(1)-m(2))**2))/2.d0/m12
      E1_cc=(m12**2-m(2)**2+m(1)**2)/2.d0/m12
        call ranlux(ran_true,1)
      cos1_cc=(dble(ran_true)-0.5d0)*2.0d0

* calculate m(3) and m12 variables in lab

      E(3)=E0/2.d0*(1.d0+(m(3)**2-m12**2)/M0**2)
     >   +P0/2.d0*((1.d0-((m12+m(3))/M0)**2)*
     >    (1.d0-((m12-m(3))/M0)**2))**0.5
     >    *cos3_c

      E12=E0/2.d0*(1.d0+(m12**2-m(3)**2)/M0**2)
     >   -P0/2.d0*((1.d0-((m12+m(3))/M0)**2)*
     >    (1.d0-((m12-m(3))/M0)**2))**0.5
     >    *cos3_c
      
      
      E3_c=(M0**2-m12**2+m(3)**2)/2.d0/M0
      alpha1=p3_c*cos3_c+V_c*E3_c
      if(alpha1.ne.0.0d0)then
         theta(3)=datan(P3_c*dsqrt(1.d0-cos3_c**2)/r_c/alpha1)
	else 
	 theta(3)=PI/2.d0
      endif
      if(theta(3).lt.0.0d0) theta(3)=theta(3)+PI

      P12=dsqrt(E12**2-m12**2)
      P(3)=dsqrt(E(3)**2-m(3)**2)
      if(P12.ne.0.0d0)then
         theta12=dacos((P0-P(3)*dcos(theta(3)))/P12)
         else
        call ranlux(ran_true,1)
            theta12=PI*dble(ran_true)
      endif
      
      E(1)=E12/2.d0*(1.d0+(m(1)**2-m(2)**2)/m12**2)
     >   +P12/2.d0*((1.d0-((m(1)+m(2))/m12)**2)
     >   *(1.d0-((m(1)-m(2))/m12)
     >   **2))**0.5*cos1_cc

      E(2)=E12/2.d0*(1.d0+(m(2)**2-m(1)**2)/m12**2)
     >   -P12/2.d0*((1.d0-((m(1)+m(2))/m12)**2)*
     >    (1.d0-((m(1)-m(2))/m12)**2
     >   ))**0.5*cos1_cc

c calculate angles of m1 and m2 in lab, with z along direction of m12
      
      alpha2= P1_cc*cos1_cc+P12/E12*E1_cc
      if(alpha2.ne.0.0d0)then     
      the12(1)=datan(P1_cc*dsqrt(1.d0-cos1_cc**2)/(E12/m12)/alpha2)
      else
         the12(1)=PI/2.d0
      endif
      if(the12(1).lt.0.0d0) the12(1)=the12(1)+PI
      P(1)=dsqrt(E(1)**2-m(1)**2)
      P(2)=dsqrt(E(2)**2-m(2)**2)
      if(P(2).ne.0.0d0)then
          the12(2) =dacos((P12-P(1)*dcos(the12(1)))/P(2))
          else
             call ranlux(ran_true,1)
            the12(2)  =PI*dble(ran_true)
       endif
 
c*** phi angles of m1 and m2 in the lab, z along diretion of P12   

        call ranlux(ran_true,1)
        phi12(1)=dble(ran_true)*2.d0*PI
        if(phi12(1).le.PI)then
           phi12(2)=phi12(1)+PI
           else
              phi12(2)=phi12(1)-PI
        endif     


c***   phi angles of m3 and m1+m2 in the lab,z along P0 of 
c***   initial decay particle momentum

          call ranlux(ran_true,1)
      phi(3)=dble(ran_true)*2.d0*PI

         if(phi(3).le.PI)then
         phi_12=phi(3)+PI
         else
            phi_12=phi(3)-PI
         endif



c convert the angles the12,phi12 in the lab frame where z along the
c direction of decay particle M


        do i=1,2
             Px_12(i)=P(i)*dsin(the12(i))*dcos(phi12(i))
             Py_12(i)=P(i)*dsin(the12(i))*dsin(phi12(i))
             Pz_12(i)=P(i)*dcos(the12(i))

             Px(i)=Pz_12(i)*dsin(theta12)*dcos(phi_12)+
     >             Px_12(i)*dcos(theta12)*dcos(phi_12)-
     >             Py_12(i)*dsin(phi_12)

             Py(i)=Pz_12(i)*dsin(theta12)*dsin(phi_12)+
     >             Px_12(i)*dcos(theta12)*dsin(phi_12)+
     >             Py_12(i)*dcos(phi_12)

             Pz(i)=Pz_12(i)*dcos(theta12)-Px_12(i)*dsin(theta12)

             if(p(i).gt.0.d0)then
                theta(i)=dacos(Pz(i)/p(i))
                p_xy(i)=dsqrt(Px(i)**2+Py(i)**2)
                if(p_xy(i).gt.0.0d0)then
                   phi(i)=dacos(px(i)/p_xy(i))
                   if(py(i).lt.0.0d0)phi(i)=2.d0*PI-phi(i)
                else
                   call ranlux(ran_true,1)
                   phi(i)=2.d0*PI*dble(ran_true)
                endif
            else
                   call ranlux(ran_true,1)  
                   theta(i)= PI*dble(ran_true)
                   call ranlux(ran_true,1)

                   phi(i)=2.d0*PI*dble(ran_true)

           endif
                

       enddo

       return
       end