How to evaluate a multidimensional NIntegral in a range where certain conditions must be satisfied by the variables

I have a function f which looks similar to this :
f[a_, b_, c_] := NIntegrate[sigma[a, b, c], {a, -valueA, valueA}, {b, -valueB, valueB}, {c, -valueC, valueC}]

and wish to integrate it only for valuesa, b, c where another function g fulfills a certain condition : g[a,b,c] >= threshold.

I did try using a boolean in this way:
f[a_, b_, c_] := NIntegrate[Boole[g[a, b, c] >= threshold]*sigma[a, b, c], {a, -valueA, valueA}, {b, -valueB, valueB}, {c, -valueC, valueC}] but I do not get the desired result.

I have also tried to define a Piecewise function for g this way and include it in the integral instead of the Boole:
Piecewise[{{g[a,b,c] , g[a,b,c]>= threshold}}]

However, I’m afraid that when using the Piecewise it gets integrated as well, which is not what I wish for. This is just a basic example and in reality I need to pass at least 3 different conditions before I integrate. Looking forward for any tips and help, it’s gonna be much appreciated.

tl;dr Trying to numerically integrate a multidmensional integral, and only pass certain values for the variables where conditions a-priori to the integration are fulfilled.

Here’s the full integral with prerequisites and values:

(*Transferred energy*)

Tmaxc12[vx_, vy_, vz_, U_, phi_, theta_] := 
 0.5*MC12 (vx^2. + vy^2. + vz^2.) + (1 - 
     Cos[theta])*(Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2] + MC12*vz)*
   Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2]/MC12 - 
  Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2]*
   Sin[theta]*(vx*Cos[phi] + vy*Sin[phi])

(*CONSTANTS DEFINITION*)

Te[U_?NumericQ] := U*e;
[Beta][U_?NumericQ] := Sqrt[1. - 1./((U/m1) + 1.)^2.];
pe[U_] := Sqrt[Te[U]*(Te[U] + 2.*m*c^2.)/c^2.];
c = 299792458.; (*speed of light*)

m = 9.10938356*10^(-31.); 

m1 = 510998.;(*electron mass in eV*)

MC12 = 12.011*1.660539040*10^(-27.); 

e = 1.60217662*10^(-19.); (*elementary charge*)
[HBar] = 
  1.054571800*10^(-34.); (*reduced Planck constant*)

Zc12 = 6.;

eps = 8.85418*10^(-12. );(*vacuum permittivity*)
(*Velocity 
distributions*)

Pvel[v_?NumericQ, Vfit_?NumericQ] := 
 1./Sqrt[2.*Pi*Vfit]*Exp[-v^2./(2.*Vfit)]
(*mean squared velocities for C12*)
VfitxyC12 = 1146080.;
VfitxC12 = VfitxyC12/2.; VfityC12 = VfitxyC12/2.; VfitzC12 = 317000.;
vxvalC12 = Sqrt[VfitxC12]; vyvalC12 = Sqrt[VfityC12]; vzvalC12 = 
 Sqrt[VfitzC12];

(*cross section*)

k1C12 = ((Zc12 e^2.)/(4. [Pi] eps 2. m c^2.))^2.;
k2C12 = [Pi] Zc12 e^2. /([HBar] c);
sigmaC12[theta_, U_] := 
  k1C12* (1. - [Beta][U]^2.) /[Beta][
    U]^4.*(Csc[theta/2.])^4.*(1. - [Beta][U]^2.*Sin[theta/2.]^2. + 
     k2C12*[Beta][U]*Sin[theta/2.] (1. - Sin[theta/2.]))*10.^28.;

This is how I defined my region of interest, where Tmax>= 21.14:

region = ImplicitRegion[
   Tmaxc12[vx, vy, vz, U, phi, theta]/e >= 
    21.14, {{vx, -vxvalC12, vxvalC12}, {vy, -vyvalC12, 
     vyvalC12}, {vz, -vzvalC12, vzvalC12}, {phi, 0, 2 Pi}, {theta, 0, 
     Pi}}];

and now the integral I was trying to solve :

sigma5D[U_] := 
 NIntegrate[ 
  sigmaC12[theta, U]*Sin[theta]*Pvel[vx, VfitxC12]*Pvel[vy, VfityC12]*
   Pvel[vz, VfitzC12], {vx, vy, vz, theta, phi} [Element] region, 
  Method -> "GlobalAdaptive"]
sigma5D[100000] // Timing

error msg:

The region given at position 1 in DiscretizeRegion[ImplicitRegion[...]] is in dimension 5. DiscretizeRegion only supports dimensions 1 through 3.

after which mathematica crashes and quits the kernel.

region[U_?NumericQ] := ImplicitRegion[
   (0.5*MC12 (vx^2. + vy^2. + vz^2.) + (1 - 
          Cos[theta])*(Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2] + MC12*vz)*
        Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2]/MC12 - 
       Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2]*
        Sin[theta]*(vx*Cos[phi] + vy*Sin[phi]))/e >= 
    21.14, {{vx, -vxvalC12, vxvalC12}, {vy, -vyvalC12, 
     vyvalC12}, {vz, -vzvalC12, vzvalC12}, {phi, 0, 2 Pi}, {theta, 0, Pi}}];

pts = RandomPoint[region[100000.], 50000];

integrand[U_, {vx_, vy_, vz_, theta_, phi_}] := 
 sigmaC12[theta, U]*Sin[theta]*Pvel[vx, VfitxC12]*Pvel[vy, VfityC12]*Pvel[vz, VfitzC12]

ListPlot[Sort[integrand[100000., #] & /@ pts], PlotRange -> All]

Mean[integrand[100000.,#]& /@ RandomPoint[region[100000.],50000]]

With[{reg = region[100000.]},
 NIntegrate[
  If[RegionMember[reg, {vx, vy, vz, theta, phi}], 
   integrand[100000., {vx, vy, vz, theta, phi}], 0], {vx, -vxvalC12, 
   vxvalC12}, {vy, -vyvalC12, vyvalC12}, {vz, -vzvalC12, 
   vzvalC12}, {phi, 0, 2 Pi}, {theta, 0, Pi}
  ]]

(* NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 1.20423050211285083223861747561433368647454170854808161214061758389`65.954589770191*^645 and 4.35609789552659774486067653532170671114285705699384650588785747247`65.954589770191*^643 for the integral and error estimates. *)

(* 1.204230502112851*10^645 *)