8
$\begingroup$

I have variables $E_{1},E_{3}$ that are bounded in the following way: $$ \tag 1 E_{1}>m_{1}, \quad E_{3}>m_{3}, E_{2}\equiv m-E_{1}-E_{3}>m_{2}, \\ -1 < \frac{E_{2}^{2}-m_{2}^{2}-(E_{1}^{2}-m_{1}^{2})-(E_{3}^{2}-m_{3}^{2})}{2\sqrt{E_{1}^{2}-m_{1}^{2}}\sqrt{E_{3}^{2}-m_{3}^{2}}}< 1 $$ Here, $m>m_{1}+m_{2}+m_{3}$, and all $m,m_{1},m_{2},m_{3}$ are positive real numbers.

I would like to generate random values of $E_{1},E_{3}$ that satisfy $(1)$. However, my code is slow. Could you please tell me how to make it faster?

(*Generating the domain (1)*)
pPar[En_, mass_] = Sqrt[En^2 - mass^2];
E2[m_, E1_, E3_] = m - E1 - E3;
costheta13[E1_, E3_, m_, m1_, m2_, 
   m3_] = (pPar[E2[m, E1, E3], m2]^2 - pPar[E1, m1]^2 - 
     pPar[E3, m3]^2)/(2*pPar[E1, m1]*pPar[E3, m3]);
E3domain[E1_, m_, m1_, m2_, m3_] = 
 E3 /. Assuming[m > 0 && E1 > m1 > 0, 
   Simplify[Solve[costheta13[E1, E3, m, m1, m2, m3] == 1, E3]]]
E1domain[m_, m1_, m2_, 
  m3_] = {E1 /. 
    Solve[E3domain[E1, m, m1, m2, m3][[1]] - 
        E3domain[E1, m, m1, m2, m3][[2]] == 0, E1][[2]], 
   E1 /. Solve[
      E3domain[E1, m, m1, m2, m3][[1]] - 
        E3domain[E1, m, m1, m2, m3][[2]] == 0, E1][[3]]} // Simplify
(*Generating random reals*)
BlockRandomVectors[m_, m1_, m2_, m3_] := 
 Block[{mass1 = m}, 
  E1rand = RandomReal[{E1domain[m, m1, m2, m3][[1]], 
     E1domain[m, m1, m2, m3][[2]]}];
  E3rand = 
   RandomReal[{E3domain[E1rand, m, m1, m2, m3][[1]], 
     E3domain[E1rand, m, m1, m2, m3][[2]]}];
  {E1rand, E3rand}]
Table[BlockRandomVectors[10, 0.5, 1, 0.2], {i, 1, 10^6, 
    1}]; // AbsoluteTiming

{84.4381,Null}

$\endgroup$
8
  • $\begingroup$ Do you want $E_1$ and $E_3$ to be distributed uniformly within these constraints? $\endgroup$
    – Roman
    Dec 2, 2021 at 14:39
  • $\begingroup$ Thanks! No, initially I wanted to generate with some weights, but them decided to simplify the task via separating. $\endgroup$ Dec 2, 2021 at 17:45
  • $\begingroup$ It looks like you're trying to generate a 4 point kinematical configuration in special relativity? I think it would be easier to generate the momenta and calculate the energies from there, why did you decide to generate the energies? $\endgroup$
    – Jojo
    Dec 2, 2021 at 20:31
  • $\begingroup$ @Joe : this is because energies specify all the relevant kinematics. $\endgroup$ Dec 2, 2021 at 20:32
  • $\begingroup$ Ok sure but there's no harm in generating the momenta even if this includes extra info that you don't need. How many dimensions are you working in? And what does the last inequality refer to physically? It looks like $\frac{|p_4|}{2|p_1||p_3|}\in (-1,1)$, which is strange because the mass dimension doesn't match up $\endgroup$
    – Jojo
    Dec 2, 2021 at 20:36

1 Answer 1

13
$\begingroup$

Use a rejection step, and compilation:

r = Compile[{m, m1, m2, m3},
      Module[{E1 = 0., E2 = 0., E3 = 0.},
        While[
          E1 = RandomReal[{m1, m - m3 - m2}];
          E3 = RandomReal[{m3, m - m1 - m2}];
          E2 = m - E1 - E3;
          E2 <= m2 || (Sqrt[E2^2-m2^2]-Sqrt[E1^2-m1^2]-Sqrt[E3^2-m3^2])^2 >= 4*(E1^2-m1^2)*(E3^2-m3^2)];
        {E1, E3}],
      CompilationTarget -> "C", RuntimeOptions -> "Speed"];

T = Table[r[10, 0.5, 1, 0.2], 10^6]; // AbsoluteTiming
(*    {0.714959, Null}    *)
$\endgroup$

Your Answer

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

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