1
$\begingroup$

The task

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}$. The procedure is given below: $$ \tag 1 \mathbf{p}' = \mathbf{p} + \gamma_{N}\mathbf{v}_{N}E + \Gamma \mathbf{v}_{N}(\mathbf{v}_{N}\cdot \mathbf{p}) $$ Here, $$ |\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})) $$ The initial parameters here are $$ 0 < \theta_{N} < \pi,\quad -\pi < \phi_{N} < \pi, \quad E_{N}>m_{N}>0, $$ $$ E>m>0, \quad 0<\theta<\pi, \quad -\pi <\phi < \pi $$ If fixing $\theta_{N},\phi_{N},m_{N},m$, the output is $$ E_{N},\theta,\phi,\mathbf{p}_{N},\mathbf{p}' $$

The code

The algorithm is as follows:

First, I compute $(1)$ (see also this question):

(*HNL 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 decay product at the rest frame of an HNL*)
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 HNL'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])];

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]; 

Then, I compute the table of values $\{E_{N},\theta,\phi\}$. This step is needed since in the full version of the code I make some manipulations with this table which require a lot of time, so it is preferable to generate it preliminarly:

mXvalue = 0.1;
mNvalue = 0.5;
EXrestValue = (mNvalue^2 + mXvalue^2)/(2*mNvalue);
ENvalue = 10;
θNvalue = 0.3;
ϕNvalue = 0.6;
ENθXϕXtable = 
  Flatten[ParallelTable[{EN, θX, ϕX}, {EN, 1.5, 20.5, 
     1}, {θX, 0.01, Pi, 0.01}, {ϕX, -Pi/2, Pi/2, 
     0.01}], {1, 2, 3}];

Finally, I compute the table:

tab = ParallelTable[{ENθXϕXtable[[i]][[1]], 
    ENθXϕXtable[[i]][[2]], 
    ENθXϕXtable[[i]][[3]], 
    pNvec[ENθXϕXtable[[i]][[1]], 
     mNvalue, θNvalue, ϕNvalue], 
    pproductLabVecCompiled[ENθXϕXtable[[i]][[1]], 
     mNvalue, θNvalue, ϕNvalue, EXrestValue, mXvalue, 
     ENθXϕXtable[[i]][[2]], 
     ENθXϕXtable[[i]][[3]]]}, {i, 1, 
    Length[ENθXϕXtable], 1}][[1]] // AbsoluteTiming 

The problem

In my machine, the time required for generating tab is 75 s. If, however, I generate is while omitting the step 2 of my algorithm I get time which is almost three times shorter:

Flatten[ParallelTable[{EN, θX, ϕX, 
     pNvec[EN, mNvalue, θNvalue, ϕNvalue], 
     pproductLabVecCompiled[EN, mNvalue, θNvalue, ϕNvalue,
       EXrestValue, mXvalue, θX, ϕX]}, {EN, 1.5, 20.5, 
     1}, {θX, 0.01, Pi, 0.01}, {ϕX, -Pi/2, Pi/2, 
     0.01}], {1, 2, 3}][[1]] // AbsoluteTiming

Maybe this is because in the algorithm with step 2 the table ENθXϕXtable is called many times.

Given that step 2 is unavoidable in my algorithm, could you please tell me how to shorten the duration of the table building with step 2?

P.S. If increasing the number of rows of tab (say, if decreasing the step in $E_{N}$ down to $0.1$), another problem occurs: Mathematica eats too much RAM (up to 10 Gb and maybe larger in my machine). This also suggests that my algorithm with step 2 included is very problematic.

$\endgroup$
3
  • $\begingroup$ It should be possible to vectorize at least one level of your loops since your arithmetic is fairly simple. Another thing: can you explain why your Γfactor contains a Solve? I cannot see why that is needed. $\endgroup$ Mar 29 at 14:28
  • $\begingroup$ @MariusLadegårdMeyer : I just need to express $v$ in terms of $\gamma = 1/\sqrt{1-v^{2}}$, so I inserted there Solve. $\endgroup$ Mar 29 at 14:36
  • $\begingroup$ You should do the algebra and solve it yourself, calling Solve many, many times would normally be a bottleneck for sure. But ah, I see you are Seting and not SetDelaying, so then it is fine. $\endgroup$ Mar 30 at 8:16
1
$\begingroup$

…This step is needed since in the full version of the code I make some manipulations with this table which require a lot of time, so it is preferable to generate it preliminarly.

I don't think so, probably there's a design miss somewhere. Anyway, here's my solution. The key idea is: Parallel* only speeds up calculation when code inside it is really slow, in other words, simple arithmetical calculation won't benefit (that much, at least) from Parallel*. Quick test:

ENθXϕXtable = 
   Flatten[ParallelTable[{EN, θX, ϕX}, {EN, 1.5, 20.5, 1}, {θX, 0.01, 
      Pi, 0.01}, {ϕX, -Pi/2, Pi/2, 0.01}], {1, 2, 3}]; // AbsoluteTiming
(* {9.59866, Null} *)

tabletst = 
   Flatten[Table[{EN, θX, ϕX}, {EN, 1.5, 20.5, 1}, {θX, 0.01, Pi, 
      0.01}, {ϕX, -Pi/2, Pi/2, 0.01}], 2]; // AbsoluteTiming
(* {0.530388, Null} *)

tabletst2 = 
   Flatten[Compile[{}, 
      Table[{EN, θX, ϕX}, {EN, 1.5, 20.5, 1}, {θX, 0.01, Pi, 
        0.01}, {ϕX, -Pi/2, Pi/2, 0.01}]][], 2]; // AbsoluteTiming *)
(* {0.237478, Null} *)

tabletst == tabletst2 == ENθXϕXtable
(* True *)

So, forget ParallelTable, once again, let's Compile. Quick fix:

cf = Hold@Compile[{{ENθXϕXtable, _Real, 2}, 
        mNvalue, θNvalue, ϕNvalue, EXrestValue, mXvalue}, Table[{
              pNvec[ENθXϕXtable[[i]][[1]], 
                mNvalue, θNvalue, ϕNvalue], 
              pproductLabVecCompiled[ENθXϕXtable[[i]][[1]], 
                mNvalue, θNvalue, ϕNvalue, EXrestValue, mXvalue, 
                ENθXϕXtable[[i]][[2]], 
                ENθXϕXtable[[i]][[3]]]}, {i, 1, 
              Length[ENθXϕXtable], 1}][[1]], 
       CompilationOptions -> "InlineCompiledFunctions" -> True] /. DownValues@pNvec /. 
    OwnValues@pproductLabVecCompiled // ReleaseHold;    

tabc = 
   cf[ENθXϕXtable, mNvalue, θNvalue, ϕNvalue, EXrestValue, 
    mXvalue]; // AbsoluteTiming
(* {2.14798, Null} *)

ENθXϕXtable[[1]]~Join~tabc
(* {1.5, 0.01, -1.5708, {0.344931, 0.23598, 1.35105}, {0.291008, 0.196689, 
  1.37983}} *)

If you have a C compiler installed, adding CompilationTarget -> C to cf will speed it up even further. The calculation of your tab takes 116.235 seconds on my laptop, BTW.

$\endgroup$
3
  • $\begingroup$ Thanks! Am I right that it is straightforward to generalize this table (cf) on the case of different input (e.g., several different tables with data in addition to ENθXϕXtable, different functions in addition to vectors...)? $\endgroup$ Mar 29 at 20:08
  • $\begingroup$ Also, a stupid question: why do we need DownValues@pNvec and OwnValues@pproductLabVecCompiled ? $\endgroup$ Mar 29 at 20:39
  • 1
    $\begingroup$ @John It depends on how the functions and tables are defined, but usually arithmetical functions are easy to compile. We need DownValues@pNvec, etc. because they're not compilable, but the functions inside them are. (If you don't know what e.g. DownValues@pNvec means, execute it separately and observe the output. ) For more info, please refer to the posts linked in the answer under your previous question. (the question isn't stupid at all BTW, it's one of the most advanced topic about the core language. ) $\endgroup$
    – xzczd
    Mar 30 at 2:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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