**The problem**

I would like to compute the Lorentz transformation $\mathbf{p}'$ of an arbitrary 3-momentum vector $\mathbf{p}$ given a boost characterized by some other momentum vector $\mathbf{p}\_{N}$, and then calculate polar and azimuthal angles $\theta',\phi'$. The procedure is given below:
$$
\mathbf{p}' = \mathbf{p} + \gamma\_{N}\mathbf{v}\_{N}E + \Gamma \mathbf{v}\_{N}(\mathbf{v}\_{N}\cdot \mathbf{p})
$$
Here,
$$
\quad |\mathbf{p}| = \sqrt{E^{2}-m^{2}}, \quad |\mathbf{p}\_{N}| = \sqrt{E\_{N}^{2}-m\_{N}^{2}},\quad \mathbf{v}\_{N} = \frac{\mathbf{p}\_{N}}{E\_{N}}, \quad \gamma\_{N} = \frac{E\_{N}}{m\_{N}}, \quad \Gamma = \frac{\gamma\_{N}-1}{v\_{N}^{2}}
$$
where in spherical coordinates
$$
\mathbf{p} = |\mathbf{p}|(s(\theta)c(\phi),s(\theta)s(\phi),c(\theta)), \quad \mathbf{p}\_{N} = |\mathbf{p}\_{N}|(s(\theta\_{N})c(\phi\_{N}),s(\theta\_{N})s(\phi\_{N}),c(\theta\_{N}))
$$
Finally,
$$
\theta' = \arccos\left( \frac{p'\_{z}}{|\mathbf{p}'|}\right), \quad \phi' = \begin{cases}\arccos\left(\frac{p\_{x}'}{\sqrt{p\_{x}^{'2}+p\_{y}^{'2}}} \right), \quad p\_{y}' > 0, \\ -\arccos\left(\frac{p\_{x}'}{\sqrt{p\_{x}^{'2}+p\_{y}^{'2}}} \right),\quad p\_{y}' < 0 \end{cases}
$$

The initial parameters here are
$$
\theta\_{N},\phi\_{N},E\_{N},m\_{N},m,E,\theta,\phi
$$

The output must be

$$
\mathbf{p}',\quad \theta', \quad \phi\
$$

**My attempt to resolve.**

I have written a piece of code that computes all these parameters. Please find it below:

```
    (*Position of the point in coordinate space in terms of spherical \
    coordinates*)
    θVal[x_, y_, z_] = ArcCos[z/Sqrt[x^2 + y^2 + z^2]];
    ϕVal[x_, y_] = 
      If[y > 0, ArcCos[x/Sqrt[x^2 + y^2]], -ArcCos[x/Sqrt[x^2 + y^2]]];
    (*N momentum, velocity and γ,Γ factors at its \
    lab frame*)
    pNvec[EN_, mN_, θN_, ϕN_] = 
      Sqrt[EN^2 - mN^2] {Sin[θN] Cos[ϕN], 
        Sin[θN] Sin[ϕN], Cos[θN]};
    vvec[EN_, mN_, θN_, ϕN_] = +(
       pNvec[EN, mN, θN, ϕN]/EN);
    γFactor[EN_, mN_] = EN/mN;
    Γfactor[EN_, mN_] = 
      Simplify[(γFactor[EN, mN] - 
        1)/((v /. 
          Solve[γFactor[EN, mN] == 1/Sqrt[1 - v^2], v])^2)[[1]]];
    (*Momentum of X at the rest frame of N*)
    pproductRestVec[EXrest_, mX_, θXrest_, ϕXrest_] = 
      Sqrt[EXrest^2 - mX^2] {Sin[θXrest] Cos[ϕXrest], 
        Sin[θXrest] Sin[ϕXrest], Cos[θXrest]};
    (*Decay product's momentum and energy at N's lab frame*)
    pproductLabVec[EN_, mN_, θN_, ϕN_, EXrest_, 
       mX_, θXrest_, ϕXrest_] = 
      Simplify[pproductRestVec[EXrest, 
         mX, θXrest, ϕXrest] + γFactor[EN, mN]*
         vvec[EN, mN, θN, ϕN]*
         EXrest + Γfactor[EN, mN]*
         vvec[EN, 
          mN, θN, ϕN] (vvec[EN, 
            mN, θN, ϕN].pproductRestVec[EXrest, 
            mX, θXrest, ϕXrest])];
    EproductLabVec[EN_, mN_, θN_, ϕN_, EXrest_, 
       mX_, θXrest_, ϕXrest_] = 
      Simplify[γFactor[EN, 
         mN] (EXrest + 
          vvec[EN, mN, θN, ϕN].pproductRestVec[EXrest, 
            mX, θXrest, ϕXrest])];
    (*X product's θ and ϕ at lab frame*)
    pproductLabVec3[EN_, mN_, θN_, ϕN_, EXrest_, 
       mX_, θXrest_, ϕXrest_] = 
      pproductLabVec[EN, mN, θN, ϕN, EXrest, 
        mX, θXrest, ϕXrest][[3]];
    pproductLabVec2[EN_, mN_, θN_, ϕN_, EXrest_, 
       mX_, θXrest_, ϕXrest_] = 
      pproductLabVec[EN, mN, θN, ϕN, EXrest, 
        mX, θXrest, ϕXrest][[2]];
    pproductLabVec1[EN_, mN_, θN_, ϕN_, EXrest_, 
       mX_, θXrest_, ϕXrest_] = 
      pproductLabVec[EN, mN, θN, ϕN, EXrest, 
        mX, θXrest, ϕXrest][[1]];
    θproductLabVec = 
 Compile[{{EN, _Real}, {mN, _Real}, {θN, _Real}, {ϕN, \
_Real}, {EXrest, _Real}, {mX, _Real}, {θXrest, _Real}, \
{ϕXrest, _Real}}, 
  ArcCos[pproductLabVec3[EN, mN, θN, ϕN, EXrest, 
     mX, θXrest, ϕXrest]/((*Norm[pproductLabVec[EN,
     mN,θN,ϕN,EXrest,
     mX,θXrest,ϕXrest]]*)√(pproductLabVec1[EN, 
         mN, θN, ϕN, EXrest, 
         mX, θXrest, ϕXrest]^2 + 
        pproductLabVec2[EN, mN, θN, ϕN, EXrest, 
         mX, θXrest, ϕXrest]^2 + 
        pproductLabVec3[EN, mN, θN, ϕN, EXrest, 
         mX, θXrest, ϕXrest]^2))]]
    ϕproductLabVec = 
     Compile[{{EN, _Real}, {mN, _Real}, {θN, _Real}, {ϕN, \
    _Real}, {EXrest, _Real}, {mX, _Real}, {θXrest, _Real}, \
    {ϕXrest, _Real}}, ϕVal[
       pproductLabVec1[EN, mN, θN, ϕN, EXrest, 
        mX, θXrest, ϕXrest], 
       pproductLabVec2[EN, mN, θN, ϕN, EXrest, 
        mX, θXrest, ϕXrest]]]

```

**Question**

However, it turns out that it works quite slow:

```
ENvalue = 20;
mNvalue = 1.5;
θNvalue = 0.03;
ϕNvalue = Pi/3;
EXrestValue = mNvalue/2;
mXvalue = 0.01;
θXrestValue = Pi/8;
ϕXrestValue = 0.1;
AbsoluteTiming[ϕproductLabVec[ENvalue, 
   mNvalue, θNvalue, ϕNvalue, EXrestValue, 
   mXvalue, θXrestValue, ϕXrestValue]]*10^6
AbsoluteTiming[θproductLabVec[ENvalue, 
   mNvalue, θNvalue, ϕNvalue, EXrestValue, 
   mXvalue, θXrestValue, ϕXrestValue]]*10^6
AbsoluteTiming[
  pproductLabVec[ENvalue, mNvalue, θNvalue, ϕNvalue, 
   EXrestValue, mXvalue, θXrestValue, ϕXrestValue]]*10^6

```

> 
> {75.1, 736183.}
> 
> 
> 

> 
> {74.7, 39506.9}
> 
> 
> 

> 
> {265.5, {564138., 511174., 1.92596\*10^7}}
> 
> 
> 

In situations when there are a lot of points (say, $10^7$), it is very crucial. Could you please point out what parts of the code may be improved in order to speed up it?