1
$\begingroup$

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?

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ You might find it helpful to simply use $\phi'=\arctan(p'_y/p'_x)$, which accounts for both the positive and negative case! Or, in mathematica, ArcTan[px, py], which accounts also for the $p'_x=0$ case. $\endgroup$
    – thorimur
    Mar 27, 2021 at 0:53
  • $\begingroup$ @John Taylor You made a typo in formula for $\bf {v}_N$ in Latex. $\endgroup$
    – Alex Trounev
    Mar 27, 2021 at 12:59
  • $\begingroup$ @AlexTrounev : thanks! $\endgroup$
    – John Taylor
    Mar 27, 2021 at 14:19

2 Answers 2

Reset to default
4
$\begingroup$

You're using Compile in wrong way. I'd suggest reading this and this post as a start. The following is a quick fix for your code:

rule = Flatten[
   DownValues /@ {pproductLabVec3, pproductLabVec2, pproductLabVec1, ϕVal}];
θproductLabVec = 
  Hold@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]/(√(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))]] //. rule // ReleaseHold;
ϕproductLabVec = 
  Hold@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]]] //. rule // ReleaseHold;

pproductLabVecCompiled = 
  Compile[{{EN, _Real}, {mN, _Real}, {θN, _Real}, {ϕN, _Real}, {EXrest, _Real}, 
           {mX, _Real}, {θXrest, _Real}, {ϕXrest, _Real}}, #] &@
   pproductLabVec[EN, mN, θN, ϕN, EXrest, mX, θXrest, ϕXrest];

Speed test:

ENvalue = 20;
mNvalue = 1.5;
θNvalue = 0.03;
ϕNvalue = Pi/3;
EXrestValue = mNvalue/2;
mXvalue = 0.01;
θXrestValue = Pi/8;
ϕXrestValue = 0.1;
AbsoluteTiming[
 Do[ϕproductLabVec[ENvalue, mNvalue, θNvalue, ϕNvalue, EXrestValue, 
   mXvalue, θXrestValue, ϕXrestValue], {10^6}]]
(* {2.217114, Null} *)

AbsoluteTiming[
 Do[θproductLabVec[ENvalue, mNvalue, θNvalue, ϕNvalue, EXrestValue, 
   mXvalue, θXrestValue, ϕXrestValue], {10^6}]]
(* {1.812695, Null} *)

AbsoluteTiming[
 Do[pproductLabVecCompiled[ENvalue, mNvalue, θNvalue, ϕNvalue, EXrestValue, 
   mXvalue, θXrestValue, ϕXrestValue], {10^6}]]
(* {1.947791, Null} *)
Improve this answer
$\endgroup$
3
$\begingroup$

The following takes approx. 4 seconds for 10^5 evaluations. A sped up could be obtained if the data is not fed one by one, but in bunches. Then instead of of calculating scalars, one could work with vectors.

boost[thn_, phn_, en_, mn_, m_, e_, ph_, th_] := 
  Module[{mn2 = mn^2, vn, gamn, gam, pt, tht, pht, pn},
   p = Sqrt[e^2 - m^2] {Sin[th] Cos[ph], Sin[th] Sin[ph], Cos[th]};
   pn = Sqrt[en^2 - mn^2] {Sin[thn] Cos[phn], Sin[thn] Sin[phn], 
      Cos[thn]};
   vn = pn/en;
   gamn = en/mn;
   gam = (gamn - 1)/vn.vn;
    
   pt = p + gamn e vn + gam (vn.p) vn;
   tht = ArcCos[pt[[3]]/Norm[pt]];
   pht = ArcTan[pt[[1]], pt[[2]]];
   
   {pt, tht, pht}
   ];

Here is a test case:

{m, mn} = RandomReal[1, {2}];
{e, en} = RandomReal[{1, 2}, {2}];
{th, thn} = RandomReal[Pi, {2}];
{ph, phn} = RandomReal[2 Pi, {2}];
p = RandomReal[1];

Do[boost[thn, phn, en, mn, m, e, ph, th], 10^5]; // Timing

(** {3.95313, Null} *)
Improve this answer
$\endgroup$
8
  • $\begingroup$ Is the convention in your code {EN, mN, θN, ϕN, EXrest, mX, θXrest, ϕXrest} -> {en, mn, thn, phn, e, m, th, ph} // Thread? If so, the output of your code doesn't match that of OP's. $\endgroup$
    – xzczd
    Mar 27, 2021 at 11:12
  • $\begingroup$ You have a typo in definition gamn, it should be gamn = en/mn. Also code by John Taylor has a typo in definition of $\Gamma$. $\endgroup$
    – Alex Trounev
    Mar 27, 2021 at 11:21
  • $\begingroup$ @Alex Trounev Thank's, nobody is perfect. $\endgroup$
    – Daniel Huber
    Mar 27, 2021 at 12:16
  • $\begingroup$ @Alex Perhaps $\mathbf{v}_{N}$ and $v_{N}$ are different? This should be clarified by OP, of course. $\endgroup$
    – xzczd
    Mar 27, 2021 at 12:34
  • $\begingroup$ @xzczd Sorry, his code has no typo, but he made typo in Latex with vn definition, it should be vn = pn/en. $\endgroup$
    – Alex Trounev
    Mar 27, 2021 at 12:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.