How to compile the following evaluation?

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[], p1Rotated[],
p1Rotated[], E1, m1}, {p2Rotated[], p2Rotated[],
p2Rotated[], E2[m, E1, E3], m2}, {p3Rotated[],
p3Rotated[], p3Rotated[], 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[], ERAND[], 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[], p1Rotated[],
p1Rotated[], E1, m1}, {p2Rotated[], p2Rotated[],
p2Rotated[], E2[m, E1, E3], m2}, {p3Rotated[],
p3Rotated[], p3Rotated[], E3, m3}},
CompilationTarget -> "C", RuntimeOptions -> "Speed"]

Namely, the code

ERAND = BlockRandomEnergies[MASS, MASS1, MASS2, MASS3];
ThreeBodyPhaseSpaceNew2Compiled[ERAND[], ERAND[], 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[] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function.

Compile::part: Part specification p1Rotated[] 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?

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