c Maps 3n-4 dimensional unit cube into n-particle phase space,
c and computes jacobian.
c Phase space is normalized in standard way (as (2.8) in fno2006)
c
c n: number of particles
c et: CM energy
c xm(n): masses
c xpar(3n-4): the 3n-4 parameters between zero and 1 that parametrize the phase space
c p(4,n): output 5 vectors: px py pz energy
c wt: output weight, wt d xpar_1 ... d xpar_{3n-4} is the phase space element
      subroutine bambo(n,et,xm,xpar,p,wt)
      implicit none
      integer n
      real * 8 et,xm(n),p(4,n),xpar(3*n-4),wt
      integer m,ipar,j,nn
      real * 8 esys,beta,vec(3),recm(4),mrec2,mrec,mrmax,mrmin,mr,
     #   mrecmin,x,xmin,xmax,phi,cth,sth,km,em,xjac,tmp
      real * 8 pi
      parameter (pi=3.141592653589793d0)

      if(n.lt.2) then
         write(*,*) ' error: cannot do 1 body phase space'
         stop
      endif

      xjac=1
c initial energy of subsystem
      esys=et
c initial boost parameters of subsystem
      beta=0
      vec(1)=1
      vec(2)=0
      vec(3)=0

c phasespace parametrized by xpar(ipar)
      ipar=1
c minimum mass of subsystem
      mrecmin=0
      do j=1,n
         mrecmin=mrecmin+xm(j)
      enddo

      do m=1,n-1
c mass of recoil system
         mrecmin=mrecmin - xm(m)
c     Recoil sistem must have mass >= sum j=2,m-1 m_j
         if(et.lt.xm(m)+mrecmin) then
            wt=0
            return
         endif
c     generate uniform 2-body phase space as if recoil system was a single
c     particle of mass from mrec to maximum, scaling as if massless
c     n-m phase space: mrec^(2*nn-2)/(2*nn-2) - mrec^(2*nn) Esys^(-2) /(2*nn)
         if(m.lt.n-1) then
            nn=(n-m)
c mr is mrec^2/esys^2
            mrmax = ((esys-xm(m))/esys)**2

c 2 body massless phase space is proportional to
c (E^2-mrec^2)*mrec^(2nn-4) d mrec^2
c      \propto Esys^2 mrec^(2*nn-2)/(2*nn-2) - mrec^(2*nn) /(2*nn)

            xmax = mrmax**(nn-1)/(nn-1)-mrmax**nn/nn
            mrmin = (mrecmin/esys)**2
            xmin = mrmin**(nn-1)/(nn-1)-mrmin**nn/nn
            x = xmax*xpar(ipar)+(1-xpar(ipar))*xmin
            xjac=xjac*(xmax-xmin)
c solves m^(n-1)/(n-1)-m^(n) / n =x
            call solvespec(nn,x,mr)
c jacobian from dx to d mr:
            xjac = xjac/(mr**(nn-2) - mr**(nn-1))
            mrec2=mr*esys**2
            xjac=xjac*esys**2
c divide by 2 pi, when introducing 2 pi delta(prec^2-mrec^2) d mrec^2/(2 pi)
            xjac=xjac/(2 * pi)
            mrec=sqrt(mrec2)
            ipar=ipar+1
         else
            mrec=mrecmin
         endif
         phi=xpar(ipar)*2*pi
         xjac=xjac*2*pi
         ipar=ipar+1
         cth=1-xpar(ipar)*2
         ipar=ipar+1
         xjac=xjac*2
         sth=sqrt(1-cth**2)
c energy of mth particle
         em=(esys**2-mrec**2+xm(m)**2)/(2*esys)
c its momentum
         km=sqrt(em**2-xm(m)**2)
c     include two body phase space factor
c     
         xjac=xjac*km/esys /(16*pi**2)
c     
         p(3,m)=km*cth
         p(1,m)=km*sth*cos(phi)
         p(2,m)=km*sth*sin(phi)
         p(4,m)=em
c     recoil momentum
         recm(1)=-p(1,m)
         recm(2)=-p(2,m)
         recm(3)=-p(3,m)
         recm(4)=esys-p(4,m)
c     Boost p_m and recm according to beta of current system
         call mboost(1,vec,beta,p(1,m),p(1,m))
         call mboost(1,vec,beta,recm,recm)
         if(m.eq.n-1) then
c nothing else to be done! last momentum is recoil momentum
            p(1,n)=recm(1)
            p(2,n)=recm(2)
            p(3,n)=recm(3)
            p(4,n)=recm(4)
         else
c     Compute velocity of recoil system; at next iteration
c     we compute the phase space of recoil system, so esys=mrec.
            tmp=sqrt(recm(1)**2+recm(2)**2+recm(3)**2)
            vec(1)=recm(1)/tmp
            vec(2)=recm(2)/tmp
            vec(3)=recm(3)/tmp
            beta=tmp/recm(4)
            esys=mrec
         endif
      enddo
      if(ipar.ne.3*n-3) then
         write(*,*) '?'
      endif
      wt=xjac
      end


      subroutine solvespec(n,x,m)
      implicit none
      real * 8 x,m,m0,dm0,vm0
      integer n
c solves the equation
c     m^(n-1)/(n-1)-m^(n) / n =x
c for integer n
      if(x.lt.0.or.x.gt.1d0/(n-1)-1d0/n) then
         m=-1
         return
      endif
      if(n.eq.2) then
         m=1-sqrt(1-2*x)
      else
         m0=0.5d0
 1       continue
         dm0=m0**(n-2)-m0**(n-1)
         vm0=m0**(n-1)/(n-1)-m0**n/n
c solve vm0+(m-m0)*dm0=x
         m=(x-vm0)/dm0+m0
         if(m.gt.1) then
            m=1
         elseif(m.lt.0) then
            m=0
         endif
         if(abs(m-m0)/(abs(m)+abs(m0)).lt.1d-13) then
            return
         else
            m0=m
            goto 1
         endif
      endif
      end