1
$\begingroup$

I have the following question. I want to generate three vectors p1Unrotated,p2Unrotated,p3Unrotated at some reference frame, then rotate them by some random angles $\theta,\phi$ around x- and z-axes using matrices phRotMatrix, thRotMatrix, and finally output the resulting vectors p1Rotated,p2Rotated,p3Rotated. My current implementation (the block BlockThreeBodyPhaseSpaceNew2) is

(*Preliminary definitions*)

(*Momenta,angles at initial frame*)
pPar[En_, mas_] = Sqrt[En^2 - mas^2];
E2[m_, E1_, E3_] = m - E1 - E3;
costh12[E1_, E3_, m_, m1_, m2_, 
   m3_] = (pPar[E3, m3]^2 - pPar[E1, m1]^2 - 
     pPar[E2[m, E1, E3], m2]^2)/(2 pPar[E1, m1]*
     pPar[E2[m, E1, E3], m2]);
costh13[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]);
th12comp = 
  Compile[{E1, E3, m, m1, m2, m3}, 
   ArcCos[costh12[E1, E3, m, m1, m2, m3]]];
th13comp = 
  Compile[{E1, E3, m, m1, m2, m3}, 
   ArcCos[costh13[E1, E3, m, m1, m2, m3]]];
(*Generator of two random values of E1,E3 available for the given \
values of m,m1,m2,m3. Please see \
https://mathematica.stackexchange.com/questions/259054/how-to-compute-\
randomreal-faster*)
BlockRandomEnergies = 
  Compile[{m, m1, m2, m3}, 
   Module[{E1 = 0., E2v = 0., E3 = 0.}, 
    While[E1 = RandomReal[{m1, (m^2 + m1^2 - (m3 + m2)^2)/(2*m)}];
     E3 = RandomReal[{m3, (m^2 + m3^2 - (m1 + m2)^2)/(2*m)}];
     E2v = m - E1 - E3;
     E2v <= m2 || (E2v^2 - m2^2 - (E1^2 - m1^2) - (E3^2 - m3^2))^2 >= 
       4*(E1^2 - m1^2)*(E3^2 - m3^2)];
    {E1, E3}], CompilationTarget -> "C", RuntimeOptions -> "Speed"];
(*Rotation matrices parametrized by theta,phi*)
phRotMatrix[
   ph_] = {{Cos[ph], Sin[ph], 0}, {-Sin[ph], Cos[ph], 0}, {0, 0, 1}};
thRotMatrix[
   th_] = {{1, 0, 0}, {0, Cos[th], Sin[th]}, {0, -Sin[th], Cos[th]}};
(*Block which rotates vectors*)
BlockThreeBodyPhaseSpaceNew2[E1_, E3_, m_, m1_, m2_, m3_] := Block[{},
  th12val = th12comp[E1, E3, m, m1, m2, m3];
  th13val = th13comp[E1, E3, m, m1, m2, m3];
  (*Momenta at frame where Subscript[p, 1] is aligned along z axis*)
  p1Unrotated = {0, 0, pPar[E1, m1]};
  p2Unrotated = {0, pPar[E2[m, E1, E3], m2]*Sin[th12val], 
    pPar[E2[m, E1, E3], m2]*Cos[th12val]};
  p3Unrotated = {0, -pPar[E3, m3]*Sin[th13val], 
    pPar[E3, m3]*Cos[th13val]};
  thrand = RandomReal[{0, Pi}];
  phrand = RandomReal[{-Pi, Pi}];
  p1Rotated = phRotMatrix[phrand].thRotMatrix[thrand].p1Unrotated;
  p2Rotated = phRotMatrix[phrand].thRotMatrix[thrand].p2Unrotated;
  p3Rotated = phRotMatrix[phrand].thRotMatrix[thrand].p3Unrotated;
  {{0., 0., 0., m, m}, {p1Rotated[[1]], p1Rotated[[2]], 
    p1Rotated[[3]], E1, m1}, {p2Rotated[[1]], p2Rotated[[2]], 
    p2Rotated[[3]], E2[m, E1, E3], m2}, {p3Rotated[[1]], 
    p3Rotated[[2]], p3Rotated[[3]], E3, m3}}]

BlockThreeBodyPhaseSpaceNew2 is slow (I need $10^{6}-10^{7}$ evaluations):

MASS = 5.3;
MASS1 = 1.;
MASS2 = 1.87;
MASS3 = 0.1;
ERAND = BlockRandomEnergies[MASS, MASS1, MASS2, MASS3];
BlockThreeBodyPhaseSpaceNew2[ERAND[[1]], ERAND[[2]], MASS, MASS1, 
   MASS2, MASS3]; // AbsoluteTiming

{0.0001341,Null}

I want to make it faster using Compile, but I do something wrong:

ThreeBodyPhaseSpaceNew2Compiled = Compile[{E1, E3, m, m1, m2, m3},
  th12val = th12comp[E1, E3, m, m1, m2, m3];
  th13val = th13comp[E1, E3, m, m1, m2, m3];
  (*Momenta at frame where Subscript[p, 1] is aligned along z axis*)
  p1Unrotated = {0, 0, pPar[E1, m1]};
  p2Unrotated = {0, pPar[E2[m, E1, E3], m2]*Sin[th12val], 
    pPar[E2[m, E1, E3], m2]*Cos[th12val]};
  p3Unrotated = {0, -pPar[E3, m3]*Sin[th13val], 
    pPar[E3, m3]*Cos[th13val]};
  thrand = RandomReal[{0, Pi}];
  phrand = RandomReal[{-Pi, Pi}];
  p1Rotated = phRotMatrix[phrand].thRotMatrix[thrand].p1Unrotated;
  p2Rotated = phRotMatrix[phrand].thRotMatrix[thrand].p2Unrotated;
  p3Rotated = phRotMatrix[phrand].thRotMatrix[thrand].p3Unrotated;
  {{0., 0., 0., m, m}, {p1Rotated[[1]], p1Rotated[[2]], 
    p1Rotated[[3]], E1, m1}, {p2Rotated[[1]], p2Rotated[[2]], 
    p2Rotated[[3]], E2[m, E1, E3], m2}, {p3Rotated[[1]], 
    p3Rotated[[2]], p3Rotated[[3]], E3, m3}}, 
  CompilationTarget -> "C", RuntimeOptions -> "Speed"]

Namely, the code

ERAND = BlockRandomEnergies[MASS, MASS1, MASS2, MASS3];
ThreeBodyPhaseSpaceNew2Compiled[ERAND[[1]], ERAND[[2]], MASS, MASS1,MASS2,MASS3]; // AbsoluteTiming

gives errors like

Compile::cplist: phRotMatrix(phrand) should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function.

Compile::part: Part specification p1Rotated[[1]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function.

Compile::part: Part specification p1Rotated[[1]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function.

It seems that I should define some tensor structures for Compile, but it is unclear to me how. This is because tensors and vectors used inside are not input, they are rather computed in-flight.

Could you please tell me how to compile the code evaluation in the correct way?

$\endgroup$
1
  • $\begingroup$ @HighPerformanceMark : thanks, that was accidentally, I have removed that question. $\endgroup$ Dec 8 '21 at 11:12

Your Answer

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

Browse other questions tagged or ask your own question.