* this routine is calculate the two body reaction m1+m2--->m3+m4
* with target m2 at rest

*************************************
*      input: m1, E1, m2, m3, m4  in lab frame
*
*      output: theta3,P3,theta4,P4   
*    
*       all angles are defined in lab frame with z axis along initial m1 direction
*       event are generated in phase space only
*        Note: all E are total energy
**********************************
      subroutine body2(m1,E1,m2,m3,m4,theta3,P3,theta4,P4)
      double precision m1,E1,m2,m3,m4,theta3,E3,theta4,E4
       double precision angcos,P,PI,E_cm,V_c,P3_c,E3_c,alpha,P3,P4
        real*4 ran_true
       
      PARAMETER (PI=3.141592653)

        call ranlux(ran_true,1)

      angcos=2.0d0*(dble(ran_true)-0.5d0)
*  energy for center of mass
      E_cm=dsqrt(m1**2+m2**2+2.0d0*E1*m2)
      
* velocity for center of mass
      P=dsqrt(E1**2-m1**2)
      V_c=P/(E1+m2)

      E3_c=(E_cm**2+m3**2-m4**2)/2.0d0/E_cm
      P3_c=dsqrt(E3_c**2-m3**2)
c      P3_c=E_cm/2.0d0*((1.0d0-((m3+m4)/E_cm)**2)*
c     >    (1.0d0-((m3-m4)/E_cm)**2))**0.5

* calculate variables in Lab

      E3=0.5d0*(E1+m2)*(1.0d0+(m3**2-m4**2)/E_cm**2)
     >   +P/2.0d0*((1.0d0-((m3+m4)/E_cm)**2)*
     >   (1.0d0-((m3-m4)/E_cm)**2))**0.5
     >   *angcos

      E4=0.5d0*(E1+m2)*(1.0d0+(m4**2-m3**2)/E_cm**2)
     >   -P/2.0d0*((1.0d0-((m3+m4)/E_cm)**2)*
     >   (1.0d0-((m3-m4)/E_cm)**2))**0.5
     >   *angcos
      
      alpha=angcos+V_c*E3_c/P3_c
      if(alpha.ne.0.0d0)then
        theta3=datan(E_cm*dsqrt(1.0d0-angcos**2)/(E1+m2)/alpha)
      else
        theta3=PI/2.0d0
      endif
      if(theta3.lt.0.0d0) theta3=theta3+PI
      
      P3=dsqrt(E3**2-m3**2)
      P4=dsqrt(E4**2-m4**2)
      if(P4.ne.0.0d0)then
         theta4=dacos((P-P3*dcos(theta3))/P4)
         else
            call ranlux(ran_true,1)
            theta4=PI*dble(ran_true)
      endif     
      
      return
      end