c This routine performs the inverse mapping from barred and radiation
c variables to the n+1 momenta, as in Sec. 5.2.1 in fno2006.
c All particle can have masses, except for the n-th and (n+1)-th.
c conventions: vector(4)=(x,y,z,t)
c Input:
c n           : number of final state barred momenta
c q0          : CM energy
c barredk(4,n): the n barred-k 4-vectors
c xi,y,phi   : the radiation variables
c Output:
c xk(4,n+1)   : the n+1 real final state momenta
c jac         : jacobian factor on phirad
      subroutine barradmap(n,q0,barredk,xi,y,phi,xk,jac)
      implicit none
c parameters
      real * 8 pi
      parameter (pi=3.141592653589793d0)
      integer n
      real * 8 q0,barredk(4,n),xi,y,phi,xk(4,n+1),jac
C Local variables
      real * 8 q2,mrec2,k0np1,uknp1,ukn,uk,cpsi,cpsi1,ubkn,vec(3),
     #     norm,k0rec,ukrec,beta,k2
      integer i
c     according to fno2006: by k0 we mean the 0 component in the CM, by
c     uk (underlined k) we mean the modulus of its 3-momentum n and np1
c     in a variable name suggests n and n+1, etc.
      q2=q0**2
c (5.42) of fnw2006
      k0np1=xi*q0/2
      uknp1=k0np1
c compute Mrec^2 (5.45)
      mrec2=(q0-barredk(4,n))**2
     #     -barredk(1,n)**2-barredk(2,n)**2-barredk(3,n)**2
      ukn=(q2-mrec2-2*q0*uknp1)/(2*(q0-uknp1*(1-y)))
c compute the length of k (5.44)
      uk=sqrt(ukn**2+uknp1**2+2*ukn*uknp1*y)
c compute cos psi (angle between knp1 and k)
      cpsi=(uk**2+uknp1**2-ukn**2)/(2*uk*uknp1)
c get the cosine of the angle between kn and k
      cpsi1=(uk**2+ukn**2-uknp1**2)/(2*uk*ukn)
c Set k_n and k_n+1 parallel to kbar_n
      ubkn=barredk(4,n)
      do i=1,4
         xk(i,n)=ukn*barredk(i,n)/ubkn
         xk(i,n+1)=uknp1*barredk(i,n)/ubkn
      enddo
c Set up a unit vector orthogonal to kbar_n and to the z axis
      vec(3)=0
      norm=sqrt(barredk(1,n)**2+barredk(2,n)**2)
      vec(1)=barredk(2,n)/norm
      vec(2)=-barredk(1,n)/norm
c Rotate k_n+1 around vec of an amount psi
      call mrotate(vec,sqrt(1-cpsi**2),cpsi,xk(1,n+1))
c Rotate k_n around vec of an amount psi1 in opposite direction
      call mrotate(vec,-sqrt(1-cpsi1**2),cpsi1,xk(1,n))
c set up a unit vector parallel to kbar_n
      do i=1,3
         vec(i)=barredk(i,n)/ubkn
      enddo
c Rotate k_n and k_n+1 around this vector of an amount phi
      call mrotate(vec,sin(phi),cos(phi),xk(1,n+1))
      call mrotate(vec,sin(phi),cos(phi),xk(1,n))
c compute the boost velocity
      k0rec=q0-ukn-uknp1
      ukrec=sqrt(k0rec**2-mrec2)
      beta=(q2-(k0rec+ukrec)**2)/(q2+(k0rec+ukrec)**2)
c     Boost all other barred k (i.e. 1 to n-1) along vec with velocity
c     beta in the k direction (same as barred k_n)
      do i=1,3
         vec(i)=barredk(i,n)/ubkn
      enddo
      call mboost(n-1,vec,beta,barredk(1,1),xk(1,1))
      k2=2*ukn*uknp1*(1-y)
      jac=q2*xi/(4*pi)**3*ukn**2/ubkn/(ukn-k2/(2*q0))
      end