"Preparations for Mathematica";

SetDirectory[NotebookDirectory[]]

"Evaluating the imaginary part of the potential in lattice perturbation theory according to [arXiv:0707.2458] Laine & Tassler (3.23)";

"Go over to new spherical coordinates for (OverVector[y]) to explicitly regularize the denominator 1/|(Overscript[y, ∼])|. We need to consider 8 different cubes with length 1/2 as the denominator 1/|(Overscript[y, ∼])| vanishes on all corners of the cube due to the Sin[(Subscript[πy, i])] nature.

In each cube we intorduce new spherical coordinates for which the measure goes proprotional to ry^2 around
a shifted origin. The shift in the origin shows up as parameter s in the definition of the coordinates.
";

y1[ry_,θy_,ϕy_,s_]:=If[s≠0,-1,1]*ry*Sin[θy]*Cos[ϕy]+If[s≠0,1,0];

y2[ry_,θy_,ϕy_,s_]:=If[s≠0,-1,1]*ry*Sin[θy]*Sin[ϕy]+If[s≠0,1,0];

y3[ry_,θy_,s_]:=If[s≠0,-1,1]*ry*Cos[θy]+If[s≠0,1,0];

MYy1[ry_,θy_,ϕy_,c_]:=Switch[c,1,y1[ry,θy,ϕy,0],2,y1[ry,θy,ϕy,1],3,y1[ry,θy,ϕy,0],4,y1[ry,θy,ϕy,0],5,y1[ry,θy,ϕy,1],6,y1[ry,θy,ϕy,1],7,y1[ry,θy,ϕy,0],8,y1[ry,θy,ϕy,1]];

MYy2[ry_,θy_,ϕy_,c_]:=Switch[c,1,y2[ry,θy,ϕy,0],2,y2[ry,θy,ϕy,0],3,y2[ry,θy,ϕy,1],4,y2[ry,θy,ϕy,0],5,y2[ry,θy,ϕy,1],6,y2[ry,θy,ϕy,0],7,y2[ry,θy,ϕy,1],8,y2[ry,θy,ϕy,1]];

MYy3[ry_,θy_,ϕy_,c_]:=Switch[c,1,y3[ry,θy,ϕy,0],2,y3[ry,θy,ϕy,0],3,y3[ry,θy,ϕy,0],4,y3[ry,θy,ϕy,1],5,y3[ry,θy,ϕy,0],6,y3[ry,θy,ϕy,1],7,y3[ry,θy,ϕy,1],8,y3[ry,θy,ϕy,1]];

JacDet[s1_,s2_,s3_]:=FullSimplify[Det[(({{D[y1[ry,θy,ϕy,s1],ry], D[y2[ry,θy,ϕy,s2],ry], D[y3[ry,θy,s3],ry]}, {D[y1[ry,θy,ϕy,s1],θy], D[y2[ry,θy,ϕy,s2],θy], D[y3[ry,θy,s3],θy]}, {D[y1[ry,θy,ϕy,s1],ϕy], D[y2[ry,θy,ϕy,s2],ϕy], D[y3[ry,θy,s3],ϕy]}}))]]

JacDet[1,1,1]

"In each cube we then need to replace the argument of the delta function. In the paper it is mentioned that one can just replace x3=x3[x1,x2,y1,y2,y3] however this leads to a pathological situation when (Overscript[(Subscript[y, 3]), ∘])=0. Instead we can choose different replacements depending on the direction of the (Overscript[y, ∘]) vector.

We choose to use the maximum's norm:

			1) max[(Overscript[(Subscript[y, 1]), ∘]),(Overscript[(Subscript[y, 2]), ∘]),(Overscript[(Subscript[y, 3]), ∘])]=(Overscript[(Subscript[y, 3]), ∘]) -> replace x3=x3[x1,x2,y1,y2,y3]
			2) max[(Overscript[(Subscript[y, 1]), ∘]),(Overscript[(Subscript[y, 2]), ∘]),(Overscript[(Subscript[y, 3]), ∘])]=(Overscript[(Subscript[y, 2]), ∘]) -> replace x2=x2[x1,x3,y1,y2,y3]
			3) max[(Overscript[(Subscript[y, 1]), ∘]),(Overscript[(Subscript[y, 2]), ∘]),(Overscript[(Subscript[y, 3]), ∘])]=(Overscript[(Subscript[y, 1]), ∘]) -> replace x1=x1[x2,x3,y1,y2,y3]
";

AccPrecGoal=10;

"Case 1:";

gFuncMaxy3[x1_,x2_,x3_,ry_,θy_,ϕy_]:=2*Sin[Pi*x1/2]*Sin[Pi*ry*Sin[θy]*Cos[ϕy]]+2*Sin[Pi*x2/2]*Sin[Pi*ry*Sin[θy]*Sin[ϕy]]-2*Sin[Pi*x3/2]*Sin[Pi*ry*Cos[θy]];

JacobianMaxy3=(({{1, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {D[gFuncMaxy3[x1,x2,x3,ry,θy,ϕy],x1], D[gFuncMaxy3[x1,x2,x3,ry,θy,ϕy],x2], D[gFuncMaxy3[x1,x2,x3,ry,θy,ϕy],x3], D[gFuncMaxy3[x1,x2,x3,ry,θy,ϕy],ry], D[gFuncMaxy3[x1,x2,x3,ry,θy,ϕy],θy], D[gFuncMaxy3[x1,x2,x3,ry,θy,ϕy],ϕy]}, {0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}}));

Abs[FullSimplify[Det[JacobianMaxy3]]]

NormalizationMaxy3[x3_?NumericQ,ry_?NumericQ,θy_?NumericQ]:=π Abs[Cos[((π x3)/2)] Sin[π ry Cos[θy]]];

AFuncMaxy3[x1_?NumericQ,x2_?NumericQ,ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ]:=2*Sin[Pi*x1/2]*Sin[Pi*ry*Sin[θy]*Cos[ϕy]]+2*Sin[Pi*x2/2]*Sin[Pi*ry*Sin[θy]*Sin[ϕy]];

x3funcMaxy3[x1_?NumericQ,x2_?NumericQ,ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ]:=(2/Pi)*ArcSin[AFuncMaxy3[x1,x2,ry,θy,ϕy]/(2*Sin[Pi*ry Cos[θy]])];

x3ArgumentMaxy3[x1_?NumericQ,x2_?NumericQ,ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ]:=AFuncMaxy3[x1,x2,ry,θy,ϕy]/(2*Sin[Pi*ry Cos[θy]]);

NormalizationMaxy3[x3_,ry_,θy_]:=π Abs[Cos[((π x3)/2)] Sin[π ry Cos[θy]]];

AFuncMaxy3[x1_,x2_,ry_,θy_,ϕy_]:=2*Sin[Pi*x1/2]*Sin[Pi*ry*Sin[θy]*Cos[ϕy]]+2*Sin[Pi*x2/2]*Sin[Pi*ry*Sin[θy]*Sin[ϕy]];

x3funcMaxy3[x1_,x2_,ry_,θy_,ϕy_]:=(2/Pi)*ArcSin[AFuncMaxy3[x1,x2,ry,θy,ϕy]/(2*Sin[Pi*ry Cos[θy]])];

x3ArgumentMaxy3[x1_,x2_,ry_,θy_,ϕy_]:=AFuncMaxy3[x1,x2,ry,θy,ϕy]/(2*Sin[Pi*ry Cos[θy]]);

NormalizationMaxy3Cartesian[x3_?NumericQ,y3_?NumericQ]:=π Abs[Cos[((π x3)/2)] Sin[π y3]];

AFuncMaxy3Cartesian[x1_?NumericQ,x2_?NumericQ,y1_?NumericQ,y2_?NumericQ]:=2*Sin[Pi*x1/2]*Sin[Pi*y1]+2*Sin[Pi*x2/2]*Sin[Pi*y2];

x3funcMaxy3Cartesian[x1_?NumericQ,x2_?NumericQ,y1_?NumericQ,y2_?NumericQ,y3_?NumericQ]:=(2/Pi)*ArcSin[AFuncMaxy3Cartesian[x1,x2,y1,y2]/(2*Sin[Pi*y3])];

x3ArgumentMaxy3Cartesian[x1_?NumericQ,x2_?NumericQ,y1_?NumericQ,y2_?NumericQ,y3_?NumericQ]:=AFuncMaxy3Cartesian[x1,x2,y1,y2]/(2*Sin[Pi*y3]);

CaseIxIntegral[ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ,r_?NumericQ]:=NIntegrate[HeavisideTheta[1-4(x3ArgumentMaxy3[x1,x2,ry,θy,ϕy]-0.5)^2]*((((3-Cos[Pi*x1*r/a]-Cos[Pi*x2*r/a]-Cos[Pi*x3funcMaxy3[x1,x2,ry,θy,ϕy]*r/a]))/((4*Sin[Pi*x1/2]^2+4*Sin[Pi*x2/2]^2+4*Sin[Pi*x3funcMaxy3[x1,x2,ry,θy,ϕy]/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)))*(1/(NormalizationMaxy3[x3funcMaxy3[x1,x2,ry,θy,ϕy],ry,θy]))*(HeavisideTheta[0.0001-ry]*(((Sin[θy] ry)/π)+(1/96) π Sin[θy] (4 ((Cos[θy])^4)+(3+Cos[4 ϕy]) ((Sin[θy])^4)) (ry^3))+HeavisideTheta[-0.0001+ry]*((ry^2*Sin[θy])/(Sqrt[(4*Sin[Pi*ry*Sin[θy]*Cos[ϕy]/2]^2+4*Sin[Pi*ry*Sin[θy]*Sin[ϕy]/2]^2+4*Sin[Pi*ry *Cos[θy]/2]^2)]))),{x1,0,1},{x2,0,1},Method->{"GlobalAdaptive"},PrecisionGoal->AccPrecGoal,AccuracyGoal->AccPrecGoal,WorkingPrecision->15,MaxRecursion->20]

(*Timing[CaseIxIntegral[0.1,0.1,0.1,a]]*)

"Case 2:";

gFuncMaxy2[x1_,x2_,x3_,ry_,θy_,ϕy_]:=2*Sin[Pi*x1/2]*Sin[Pi*ry*Sin[θy]*Cos[ϕy]]+2*Sin[Pi*x2/2]*Sin[Pi*ry*Sin[θy]*Sin[ϕy]]-2*Sin[Pi*x3/2]*Sin[Pi*ry*Cos[θy]];

JacobianMaxy2=(({{1, 0, 0, 0, 0, 0}, {D[gFuncMaxy2[x1,x2,x3,ry,θy,ϕy],x1], D[gFuncMaxy2[x1,x2,x3,ry,θy,ϕy],x2], D[gFuncMaxy2[x1,x2,x3,ry,θy,ϕy],x3], D[gFuncMaxy2[x1,x2,x3,ry,θy,ϕy],ry], D[gFuncMaxy2[x1,x2,x3,ry,θy,ϕy],θy], D[gFuncMaxy2[x1,x2,x3,ry,θy,ϕy],ϕy]}, {0, 0, 1, 0, 0, 0}, {0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}}));

Abs[FullSimplify[Det[JacobianMaxy2]]]

NormalizationMaxy2[x2_?NumericQ,ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ]:=π Abs[Cos[((π x2)/2)] Sin[π ry Sin[θy] Sin[ϕy]]];

AFuncMaxy2[x1_?NumericQ,x3_?NumericQ,ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ]:=2*Sin[Pi*x1/2]*Sin[Pi*ry*Sin[θy]*Cos[ϕy]]-2*Sin[Pi*x3/2]*Sin[Pi*ry*Cos[θy]];

x2funcMaxy2[x1_?NumericQ,x3_?NumericQ,ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ]:=(2/Pi)*ArcSin[-AFuncMaxy2[x1,x3,ry,θy,ϕy]/(2*Sin[Pi*ry*Sin[θy]*Sin[ϕy]])];

x2ArgumentMaxy2[x1_?NumericQ,x3_?NumericQ,ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ]:=-AFuncMaxy2[x1,x3,ry,θy,ϕy]/(2*Sin[Pi*ry*Sin[θy]*Sin[ϕy]]);

NormalizationMaxy2Cartesian[x2_?NumericQ,y2_?NumericQ]:=π Abs[Cos[((π x2)/2)] Sin[π y2]];

AFuncMaxy2Cartesian[x1_?NumericQ,x3_?NumericQ,y1_?NumericQ,y3_?NumericQ]:=2*Sin[Pi*x1/2]*Sin[Pi*y1]-2*Sin[Pi*x3/2]*Sin[Pi*y3];

x2funcMaxy2Cartesian[x1_?NumericQ,x3_?NumericQ,y1_?NumericQ,y2_?NumericQ,y3_?NumericQ]:=(2/Pi)*ArcSin[-AFuncMaxy2Cartesian[x1,x3,y1,y3]/(2*Sin[Pi*y2])];

x2ArgumentMaxy2Cartesian[x1_?NumericQ,x3_?NumericQ,y1_?NumericQ,y2_?NumericQ,y3_?NumericQ]:=-AFuncMaxy2Cartesian[x1,x3,y1,y3]/(2*Sin[Pi*y2]);

CaseIIxIntegral[ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ,r_?NumericQ]:=NIntegrate[HeavisideTheta[1-4( x2ArgumentMaxy2[x1,x3,ry,θy,ϕy]-0.5)^2]*((((3-Cos[Pi*x1*r/a]-Cos[Pi*x2funcMaxy2[x1,x3,ry,θy,ϕy]*r/a]-Cos[Pi*x3*r/a]))/((4*Sin[Pi*x1/2]^2+4*Sin[Pi*x2funcMaxy2[x1,x3,ry,θy,ϕy]/2]^2+4*Sin[Pi*x3/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)))*(1/(NormalizationMaxy2[x2funcMaxy2[x1,x3,ry,θy,ϕy],ry,θy,ϕy]))*If[ry<0.0001,((Sin[θy] ry)/π)+(1/96) π Sin[θy] (4 ((Cos[θy])^4)+(3+Cos[4 ϕy]) ((Sin[θy])^4)) (ry^3),((ry^2*Sin[θy])/(Sqrt[(4*Sin[Pi*ry*Sin[θy]*Cos[ϕy]/2]^2+4*Sin[Pi*ry*Sin[θy]*Sin[ϕy]/2]^2+4*Sin[Pi*ry *Cos[θy]/2]^2)]))],{x1,0,1},{x3,0,1},Method->{"GlobalAdaptive"},PrecisionGoal->AccPrecGoal,AccuracyGoal->AccPrecGoal,WorkingPrecision->15,MaxRecursion->20]

(*Timing[CaseIIxIntegral[0.1,0.1,0.1,a]]*)

"Case 3:";

gFuncMaxy1[x1_,x2_,x3_,ry_,θy_,ϕy_]:=2*Sin[Pi*x1/2]*Sin[Pi*ry*Sin[θy]*Cos[ϕy]]+2*Sin[Pi*x2/2]*Sin[Pi*ry*Sin[θy]*Sin[ϕy]]-2*Sin[Pi*x3/2]*Sin[Pi*ry*Cos[θy]];

JacobianMaxy1=(({{D[gFuncMaxy1[x1,x2,x3,ry,θy,ϕy],x1], D[gFuncMaxy1[x1,x2,x3,ry,θy,ϕy],x2], D[gFuncMaxy1[x1,x2,x3,ry,θy,ϕy],x3], D[gFuncMaxy1[x1,x2,x3,ry,θy,ϕy],ry], D[gFuncMaxy1[x1,x2,x3,ry,θy,ϕy],θy], D[gFuncMaxy1[x1,x2,x3,ry,θy,ϕy],ϕy]}, {0, 1, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0}, {0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}}));

Abs[FullSimplify[Det[JacobianMaxy1]]]

NormalizationMaxy1[x1_?NumericQ,ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ]:=π Abs[Cos[((π x1)/2)] Sin[π ry Cos[ϕy] Sin[θy]]];

AFuncMaxy1[x2_?NumericQ,x3_?NumericQ,ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ]:=2*Sin[Pi*x2/2]*Sin[Pi*ry*Sin[θy]*Sin[ϕy]]-2*Sin[Pi*x3/2]*Sin[Pi*ry*Cos[θy]];

x1funcMaxy1[x2_?NumericQ,x3_?NumericQ,ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ]:=(2/Pi)*ArcSin[-AFuncMaxy1[x2,x3,ry,θy,ϕy]/(2*Sin[Pi*ry*Sin[θy]*Cos[ϕy]])];

x1ArgumentMaxy1[x2_?NumericQ,x3_?NumericQ,ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ]:=-AFuncMaxy1[x2,x3,ry,θy,ϕy]/(2*Sin[Pi*ry*Sin[θy]*Sin[ϕy]]);

NormalizationMaxy1Cartesian[x1_?NumericQ,y2_?NumericQ]:=π Abs[Cos[((π x1)/2)] Sin[π y2]];

AFuncMaxy1Cartesian[x2_?NumericQ,x3_?NumericQ,y2_?NumericQ,y3_?NumericQ]:=2*Sin[Pi*x2/2]*Sin[Pi*y2]-2*Sin[Pi*x3/2]*Sin[Pi*y3];

x1funcMaxy1Cartesian[x2_?NumericQ,x3_?NumericQ,y1_?NumericQ,y2_?NumericQ,y3_?NumericQ]:=(2/Pi)*ArcSin[-AFuncMaxy1Cartesian[x2,x3,y2,y3]/(2*Sin[Pi*y1])];

x1ArgumentMaxy1Cartesian[x2_?NumericQ,x3_?NumericQ,y1_?NumericQ,y2_?NumericQ,y3_?NumericQ]:=-AFuncMaxy1Cartesian[x2,x3,y2,y3]/(2*Sin[Pi*y1]);

CaseIIIxIntegral[ry_?NumericQ,θy_?NumericQ,ϕy_?NumericQ,r_?NumericQ]:=NIntegrate[HeavisideTheta[1-4( x1ArgumentMaxy1[x2,x3,ry,θy,ϕy]-0.5)^2]*((((3-Cos[Pi*x1funcMaxy1[x2,x3,ry,θy,ϕy]*r/a]-Cos[Pi*x2*r/a]-Cos[Pi*x3*r/a]))/((4*Sin[Pi*x1funcMaxy1[x2,x3,ry,θy,ϕy]/2]^2+4*Sin[Pi*x2/2]^2+4*Sin[Pi*x3/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)))*(1/(NormalizationMaxy1[x1funcMaxy1[x2,x3,ry,θy,ϕy],ry,θy,ϕy]))*If[ry<0.0001,((Sin[θy] ry)/π)+(1/96) π Sin[θy] (4 ((Cos[θy])^4)+(3+Cos[4 ϕy]) ((Sin[θy])^4)) (ry^3),((ry^2*Sin[θy])/(Sqrt[(4*Sin[Pi*ry*Sin[θy]*Cos[ϕy]/2]^2+4*Sin[Pi*ry*Sin[θy]*Sin[ϕy]/2]^2+4*Sin[Pi*ry *Cos[θy]/2]^2)]))],{x2,0,1},{x3,0,1},Method->{"GlobalAdaptive"},PrecisionGoal->AccPrecGoal,AccuracyGoal->AccPrecGoal,WorkingPrecision->15,MaxRecursion->20]

(*Timing[CaseIIIxIntegral[0.1,0.1,0.1,a]]*)

"Numerical implementation of the integral";

SigmaLat=Gamma[1/24]^2*Gamma[11/24]^2*(Sqrt[3]-1)/(48*Pi^2);

Nc=3;

CA=Nc;

CF=((Nc^2-1)/(2 Nc));

Clear[a]

a=1;

"
A) Since y integral is over the 3d unit cube, let the speherical integration go over the ball with radius Sqrt[3].
B) Set integrand zero outside of the unit cube.
C) Che3ck which case is fullfilled and provide the corresponding integrand with the replaced integration variables
D) Check that the argument of the delta function vanishes by making sure that the argument of the ArcSin lies within the correct range of valued 0 <= arg <=1 so that the ArcSin remains real valued.
";

FastMaximum[x_]:=#&@@Ordering[-x,1]

"Cube 1: 0<y1<0.5  0<y2<0.5 0<y3<0.5 ";

LatImVCube1QMC[betaLat_,r_,ag_]:=-((2*Pi*CF*CA^2)/betaLat)*NIntegrate[ 
If[y1[ry,θy,ϕy,0]>0.5||y2[ry,θy,ϕy,0]>0.5||y3[ry,θy,0]>0.5,0,

Switch[ FastMaximum[{Sin[Pi*y3[ry,θy,0]],Sin[Pi*y2[ry,θy,ϕy,0]],Sin[Pi*y1[ry,θy,ϕy,0]]}],
1,

If[ Between[x3ArgumentMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]],{0,1}],

((3-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a]-Cos[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]]*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+4*Sin[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]]/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]],y3[ry,θy,0]]),0 ],

2,

If[Between[x2ArgumentMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]],{0,1}],((3-Cos[Pi*xA*r/a]-Cos[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]]*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]]/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]],y2[ry,θy,ϕy,0]])),0],

3,

If[Between[ x1ArgumentMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]],{0,1}],((3-Cos[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]]*r/a]-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]]/2]^2+4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]],y1[ry,θy,ϕy,0]])),0]

]*((ry^2) Sin[θy])*(1/(Sqrt[(4*Sin[Pi*y1[ry,θy,ϕy,0]/2]^2+4*Sin[Pi*y2[ry,θy,ϕy,0]/2]^2+4*Sin[Pi*y3[ry,θy,0]/2]^2)]))], {ry,0,1},{ϕy,0,Pi/2},{θy,0,Pi/2},{xA,0,1},{xB,0,1},Method->{"AdaptiveQuasiMonteCarlo","MaxPoints"->10^8},PrecisionGoal->10,AccuracyGoal->ag,Exclusions->{NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]],y3[ry,θy,0]]==0,NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]],y2[ry,θy,ϕy,0]]==0,NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,0]],y1[ry,θy,ϕy,0]]==0}]

(*LatImVCube1QMC[a,3]*)

"Cube 2: 0.5<y1<1  0<y2<0.5 0<y3<0.5 ";

LatImVCube2QMC[betaLat_,r_,ag_]:=-((2*Pi*CF*CA^2)/betaLat)*NIntegrate[ 
If[y1[ry,θy,ϕy,1]<0.5||y1[ry,θy,ϕy,1]>1||y2[ry,θy,ϕy,0]>0.5||y3[ry,θy,0]>0.5,0,

Switch[ FastMaximum[{Sin[Pi*y3[ry,θy,0]],Sin[Pi*y2[ry,θy,ϕy,0]],Sin[Pi*y1[ry,θy,ϕy,1]]}],
1,

If[ Between[x3ArgumentMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]],{0,1}],

((3-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a]-Cos[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]]*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+4*Sin[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]]/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]],y3[ry,θy,0]]),0 ],

2,

If[Between[x2ArgumentMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]],{0,1}],((3-Cos[Pi*xA*r/a]-Cos[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]]*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]]/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]],y2[ry,θy,ϕy,0]])),0],

3,

If[Between[ x1ArgumentMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]],{0,1}],((3-Cos[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]]*r/a]-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]]/2]^2+4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]],y1[ry,θy,ϕy,1]])),0]

]*((ry^2) Sin[θy])*(1/(Sqrt[(4*Sin[Pi*y1[ry,θy,ϕy,1]/2]^2+4*Sin[Pi*y2[ry,θy,ϕy,0]/2]^2+4*Sin[Pi*y3[ry,θy,0]/2]^2)]))], {ry,0,1},{ϕy,0,Pi/2},{θy,0,Pi/2},{xA,0,1},{xB,0,1},Method->{"AdaptiveQuasiMonteCarlo","MaxPoints"->10^8},PrecisionGoal->10,AccuracyGoal->ag,Exclusions->{NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]],y3[ry,θy,0]]==0,NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]],y2[ry,θy,ϕy,0]]==0,NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,0]],y1[ry,θy,ϕy,1]]==0}]

(*LatImVCube2QMC[a,3]*)

"Cube 3: 0<y1<0.5  0.5<y2<1 0<y3<0.5 ";

LatImVCube3QMC[betaLat_,r_,ag_]:=-((2*Pi*CF*CA^2)/betaLat)*NIntegrate[ 
If[y1[ry,θy,ϕy,0]>0.5||y2[ry,θy,ϕy,1]<0.5||y2[ry,θy,ϕy,1]>1||y3[ry,θy,0]>0.5,0,

Switch[ FastMaximum[{Sin[Pi*y3[ry,θy,0]],Sin[Pi*y2[ry,θy,ϕy,1]],Sin[Pi*y1[ry,θy,ϕy,0]]}],
1,

If[ Between[x3ArgumentMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]],{0,1}],

((3-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a]-Cos[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]]*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+4*Sin[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]]/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]],y3[ry,θy,0]]),0 ],

2,

If[Between[x2ArgumentMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]],{0,1}],((3-Cos[Pi*xA*r/a]-Cos[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]]*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]]/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]],y2[ry,θy,ϕy,1]])),0],

3,

If[Between[ x1ArgumentMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]],{0,1}],((3-Cos[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]]*r/a]-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]]/2]^2+4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]],y1[ry,θy,ϕy,0]])),0]

]*((ry^2) Sin[θy])*(1/(Sqrt[(4*Sin[Pi*y1[ry,θy,ϕy,0]/2]^2+4*Sin[Pi*y2[ry,θy,ϕy,1]/2]^2+4*Sin[Pi*y3[ry,θy,0]/2]^2)]))], {ry,0,1},{ϕy,0,Pi/2},{θy,0,Pi/2},{xA,0,1},{xB,0,1},Method->{"AdaptiveQuasiMonteCarlo","MaxPoints"->10^8},PrecisionGoal->10,AccuracyGoal->ag,Exclusions->{NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]],y3[ry,θy,0]]==0,NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]],y2[ry,θy,ϕy,1]]==0,NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,0]],y1[ry,θy,ϕy,0]]==0}]

(*LatImVCube3QMC[a,3]*)

"Cube 4: 0<y1<0.5  0<y2<0.5 0.5<y3<1 ";

LatImVCube4QMC[betaLat_,r_,ag_]:=-((2*Pi*CF*CA^2)/betaLat)*NIntegrate[ 
If[y1[ry,θy,ϕy,0]>0.5||y2[ry,θy,ϕy,0]>0.5||y3[ry,θy,1]<0.5||y3[ry,θy,1]>1,0,

Switch[ FastMaximum[{Sin[Pi*y3[ry,θy,1]],Sin[Pi*y2[ry,θy,ϕy,0]],Sin[Pi*y1[ry,θy,ϕy,0]]}],
1,

If[ Between[x3ArgumentMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]],{0,1}],

((3-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a]-Cos[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]]*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+4*Sin[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]]/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]],y3[ry,θy,1]]),0 ],

2,

If[Between[x2ArgumentMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]],{0,1}],((3-Cos[Pi*xA*r/a]-Cos[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]]*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]]/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]],y2[ry,θy,ϕy,0]])),0],

3,

If[Between[ x1ArgumentMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]],{0,1}],((3-Cos[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]]*r/a]-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]]/2]^2+4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]],y1[ry,θy,ϕy,0]])),0]

]*((ry^2) Sin[θy])*(1/(Sqrt[(4*Sin[Pi*y1[ry,θy,ϕy,0]/2]^2+4*Sin[Pi*y2[ry,θy,ϕy,0]/2]^2+4*Sin[Pi*y3[ry,θy,1]/2]^2)]))], {ry,0,1},{ϕy,0,Pi/2},{θy,0,Pi/2},{xA,0,1},{xB,0,1},Method->{"AdaptiveQuasiMonteCarlo","MaxPoints"->10^8},PrecisionGoal->10,AccuracyGoal->ag,Exclusions->{NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]],y3[ry,θy,1]]==0,NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]],y2[ry,θy,ϕy,0]]==0,NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,0],y3[ry,θy,1]],y1[ry,θy,ϕy,0]]==0}]

(*LatImVCube4QMC[a,3]*)

"Cube 5: 0.5<y1<1  0.5<y2<1 0<y3<0.5 ";

LatImVCube5QMC[betaLat_,r_,ag_]:=-((2*Pi*CF*CA^2)/betaLat)*NIntegrate[ 
If[y1[ry,θy,ϕy,1]<0.5||y1[ry,θy,ϕy,1]>1||y2[ry,θy,ϕy,1]<0.5||y2[ry,θy,ϕy,1]>1||y3[ry,θy,0]>0.5,0,

Switch[ FastMaximum[{Sin[Pi*y3[ry,θy,0]],Sin[Pi*y2[ry,θy,ϕy,1]],Sin[Pi*y1[ry,θy,ϕy,1]]}],
1,

If[ Between[x3ArgumentMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]],{0,1}],

((3-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a]-Cos[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]]*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+4*Sin[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]]/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]],y3[ry,θy,0]]),0 ],

2,

If[Between[x2ArgumentMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]],{0,1}],((3-Cos[Pi*xA*r/a]-Cos[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]]*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]]/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]],y2[ry,θy,ϕy,1]])),0],

3,

If[Between[ x1ArgumentMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]],{0,1}],((3-Cos[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]]*r/a]-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]]/2]^2+4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]],y1[ry,θy,ϕy,1]])),0]

]*((ry^2) Sin[θy])*(1/(Sqrt[(4*Sin[Pi*y1[ry,θy,ϕy,1]/2]^2+4*Sin[Pi*y2[ry,θy,ϕy,1]/2]^2+4*Sin[Pi*y3[ry,θy,0]/2]^2)]))], {ry,0,1},{ϕy,0,Pi/2},{θy,0,Pi/2},{xA,0,1},{xB,0,1},Method->{"AdaptiveQuasiMonteCarlo","MaxPoints"->10^8},PrecisionGoal->10,AccuracyGoal->ag,Exclusions->{NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]],y3[ry,θy,0]]==0,NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]],y2[ry,θy,ϕy,1]]==0,NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,0]],y1[ry,θy,ϕy,1]]==0}]

(*LatImVCube5QMC[a,3]*)

"Cube 6: 0.5<y1<1  0<y2<0.5 0.5<y3<1 ";

LatImVCube6QMC[betaLat_,r_,ag_]:=-((2*Pi*CF*CA^2)/betaLat)*NIntegrate[ 
If[y1[ry,θy,ϕy,1]<0.5||y1[ry,θy,ϕy,1]>1||y2[ry,θy,ϕy,0]>0.5||y3[ry,θy,1]<0.5||y3[ry,θy,0]>1,0,

Switch[ FastMaximum[{Sin[Pi*y3[ry,θy,1]],Sin[Pi*y2[ry,θy,ϕy,0]],Sin[Pi*y1[ry,θy,ϕy,1]]}],
1,

If[ Between[x3ArgumentMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]],{0,1}],

((3-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a]-Cos[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]]*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+4*Sin[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]]/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]],y3[ry,θy,1]]),0 ],

2,

If[Between[x2ArgumentMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]],{0,1}],((3-Cos[Pi*xA*r/a]-Cos[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]]*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]]/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]],y2[ry,θy,ϕy,0]])),0],

3,

If[Between[ x1ArgumentMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]],{0,1}],((3-Cos[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]]*r/a]-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]]/2]^2+4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]],y1[ry,θy,ϕy,1]])),0]

]*((ry^2) Sin[θy])*(1/(Sqrt[(4*Sin[Pi*y1[ry,θy,ϕy,1]/2]^2+4*Sin[Pi*y2[ry,θy,ϕy,0]/2]^2+4*Sin[Pi*y3[ry,θy,1]/2]^2)]))], {ry,0,1},{ϕy,0,Pi/2},{θy,0,Pi/2},{xA,0,1},{xB,0,1},Method->{"AdaptiveQuasiMonteCarlo","MaxPoints"->10^8},PrecisionGoal->10,AccuracyGoal->ag,Exclusions->{NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]],y3[ry,θy,1]]==0,NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]],y2[ry,θy,ϕy,0]]==0,NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,0],y3[ry,θy,1]],y1[ry,θy,ϕy,1]]==0}]

(*LatImVCube6QMC[a,3]*)

"Cube 7: 0<y1<0.5  0.5<y2<1 0.5<y3<1 ";

LatImVCube7QMC[betaLat_,r_,ag_]:=-((2*Pi*CF*CA^2)/betaLat)*NIntegrate[ 
If[y1[ry,θy,ϕy,0]>0.5||y2[ry,θy,ϕy,1]<0.5||y2[ry,θy,ϕy,1]>1||y3[ry,θy,1]<0.5||y3[ry,θy,1]>1,0,

Switch[ FastMaximum[{Sin[Pi*y3[ry,θy,1]],Sin[Pi*y2[ry,θy,ϕy,1]],Sin[Pi*y1[ry,θy,ϕy,0]]}],
1,

If[ Between[x3ArgumentMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]],{0,1}],

((3-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a]-Cos[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]]*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+4*Sin[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]]/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]],y3[ry,θy,1]]),0 ],

2,

If[Between[x2ArgumentMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]],{0,1}],((3-Cos[Pi*xA*r/a]-Cos[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]]*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]]/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]],y2[ry,θy,ϕy,1]])),0],

3,

If[Between[ x1ArgumentMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]],{0,1}],((3-Cos[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]]*r/a]-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]]/2]^2+4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]],y1[ry,θy,ϕy,0]])),0]

]*((ry^2) Sin[θy])*(1/(Sqrt[(4*Sin[Pi*y1[ry,θy,ϕy,0]/2]^2+4*Sin[Pi*y2[ry,θy,ϕy,1]/2]^2+4*Sin[Pi*y3[ry,θy,1]/2]^2)]))], {ry,0,1},{ϕy,0,Pi/2},{θy,0,Pi/2},{xA,0,1},{xB,0,1},Method->{"AdaptiveQuasiMonteCarlo","MaxPoints"->10^8},PrecisionGoal->10,AccuracyGoal->ag,Exclusions->{NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]],y3[ry,θy,1]]==0,NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]],y2[ry,θy,ϕy,1]]==0,NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,0],y2[ry,θy,ϕy,1],y3[ry,θy,1]],y1[ry,θy,ϕy,0]]==0}]

(*LatImVCube7QMC[a,3]*)

"Cube 8: 0.5<y1<1  0.5<y2<1 0.5<y3<1 ";

LatImVCube8QMC[betaLat_,r_,ag_]:=-((2*Pi*CF*CA^2)/betaLat)*NIntegrate[ 
If[y1[ry,θy,ϕy,1]<0.5||y1[ry,θy,ϕy,1]>1||y2[ry,θy,ϕy,1]<0.5||y2[ry,θy,ϕy,1]>1.||y3[ry,θy,1]<0.5||y3[ry,θy,1]>1,0,

Switch[ FastMaximum[{Sin[Pi*y3[ry,θy,1]],Sin[Pi*y2[ry,θy,ϕy,1]],Sin[Pi*y1[ry,θy,ϕy,1]]}],
1,

If[ Between[x3ArgumentMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]],{0,1}],

((3-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a]-Cos[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]]*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+4*Sin[Pi*x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]]/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]],y3[ry,θy,1]]),0 ],

2,

If[Between[x2ArgumentMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]],{0,1}],((3-Cos[Pi*xA*r/a]-Cos[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]]*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*xA/2]^2+4*Sin[Pi*x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]]/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]],y2[ry,θy,ϕy,1]])),0],

3,

If[Between[ x1ArgumentMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]],{0,1}],((3-Cos[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]]*r/a]-Cos[Pi*xA*r/a]-Cos[Pi*xB*r/a])/(4*Sin[Pi*x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]]/2]^2+4*Sin[Pi*xA/2]^2+4*Sin[Pi*xB/2]^2+CA^2*SigmaLat/(Pi*betaLat))^2)*(1/(NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]],y1[ry,θy,ϕy,1]])),0]

]*((ry^2) Sin[θy])*(1/(Sqrt[(4*Sin[Pi*y1[ry,θy,ϕy,1]/2]^2+4*Sin[Pi*y2[ry,θy,ϕy,1]/2]^2+4*Sin[Pi*y3[ry,θy,1]/2]^2)]))], {ry,0,1},{ϕy,0,Pi/2},{θy,0,Pi/2},{xA,0,1},{xB,0,1},Method->{"AdaptiveQuasiMonteCarlo","MaxPoints"->10^8},PrecisionGoal->10,AccuracyGoal->ag,Exclusions->{NormalizationMaxy3Cartesian[x3funcMaxy3Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]],y3[ry,θy,1]]==0,NormalizationMaxy2Cartesian[x2funcMaxy2Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]],y2[ry,θy,ϕy,1]]==0,NormalizationMaxy1Cartesian[x1funcMaxy1Cartesian[xA,xB,y1[ry,θy,ϕy,1],y2[ry,θy,ϕy,1],y3[ry,θy,1]],y1[ry,θy,ϕy,1]]==0}]

(*LatImVCube8QMC[a,3]*)

"Now combine all eight cubes and compute the final result:";

FullImVAllCorners[betaLat_,r_,ag_]:=LatImVCube1QMC[betaLat,r,ag]+LatImVCube2QMC[betaLat,r,ag]+LatImVCube3QMC[betaLat,r,ag]+LatImVCube4QMC[betaLat,r,ag]+LatImVCube5QMC[betaLat,r,ag]+LatImVCube6QMC[betaLat,r,ag]+LatImVCube7QMC[betaLat,r,ag]+LatImVCube8QMC[betaLat,r,ag]