Parallelization of Monte Carlo numerical integration

Question

I'm trying to parallelize some integral. What I do: I define integrand:

Integrand105[p0_, p1_, t1_] =  2 (2*Pi)^4 (2/Pi)^4*(p0*
     p1 (ArcTan[p0/2]/(4*Pi)) (ArcTan[
         p1/2]/(4*Pi))/((p0^2 - 2*p0*p1*t1 + p1^2 + 1)^4*(p0^2 + 
          1)*(p1^2 + 1)));

Give a definition to all the kernels:

LaunchKernels[];
DistributeDefinitions[Integrand105];

The integral itself:

res = Total[
   ParallelTable[
    NIntegrate[
     Integrand105[p0, p1, t1], {p0, 0, Infinity}, {p1, 0, 
      Infinity}, {t1, -1 + i, i}, Method -> "MonteCarlo", 
     MaxPoints -> 1000000, PrecisionGoal -> 5], {i, 0, 1}]] // 
  AbsoluteTiming

That means that all parallelization consists in splitting the domain of integration. Ok. What I get: - from first kernel:

NIntegrate::maxp: The integral failed to converge after 1000100 integrand evaluations. NIntegrate obtained 0.0024072378420660233` and 5.0413357069178085`*^-6 for the integral and error estimates.

-from second one:

NIntegrate::maxp: The integral failed to converge after 1000100 integrand evaluations. NIntegrate obtained 0.05468151281718856` and 0.0007526057830408755` for the integral and error estimates.

And in the end:

{3.569,0.05708}

So, my first question is: How can I get common estimation for error? Of course, I can sum these two error bars, but it will be the highest border. Any ideas? But even with this way, unfortunately, I have troubles:

I found hear: NIntegrate error bound how I can obtain errors:

Unprotect[Message];
Message[NIntegrate::maxp, its_, int_, err_] := Sow[err]

If I apply it just to integral:

NIntegrate[
  Integrand105[p0, p1, t1], {p0, 0, Infinity}, {p1, 0, 
   Infinity}, {t1, -1, 1}, Method -> "MonteCarlo", 
  MaxPoints -> 2000000, PrecisionGoal -> 10] // Reap

I get:

{0.0583558, {{0.0007278}}}

Excellent, I thought! I did all kinds of things to apply it for my case:

Total[ParallelTable[
   NIntegrate[
     Integrand105[p0, p1, t1], {p0, 0, Infinity}, {p1, 0, 
      Infinity}, {t1, -1 + i, i}, Method -> "MonteCarlo", 
     MaxPoints -> 1000000, PrecisionGoal -> 5] // Reap, {i, 0, 
    1}]] // AbsoluteTiming 

And... {5.217, {0.05672, {}}}

Ok, may be without Total - no - {5.69, {{0.002407, {}}, {0.05574, {}}}}.

Also I tried other ways, but I didn't find smth that would help me.

I believe in the wisdom of the community. Thank you very much in advance.

p.s. In common, I'm looking for a way to calculate a integral(not like in the example above, integral as in the example one can calculate by means of, for instance:

Method -> {"GlobalAdaptive", Method -> "CartesianRule"},PrecisionGoal -> 15

or with:

"SingularityHandler" -> "DuffyCoordinates"}, PrecisionGoal -> 15

Within 5 minutes on a laptop. But applying this method, say, to 7-dimensional integral you 100% fail (If I'm wrong, please, make me the happiest man and correct me). I have an access to a cluster and I try to find an approach to calculate this(7-dim) integral with highest as far as possible precision.