Numerical integration of MeijerG function with variables

The root cause of the difficulties encountered by the OP is that MeijerG, as implemented in Mathematica, often is not accurate when evaluated at machine precision. For instance, with the parameters,

params = {p -> 1, q1 -> 1, q2 -> 1, q13 -> 1, qi1 -> 1, mu -> 3/2, 
    L -> 3, gammak1 -> 1, BI -> 9, C1 -> 4, Bsd -> 4, M -> 2};

merely plotting the integrand in the Question over a portion of its domain

mg = MeijerG[{{-(q13 + q1 + q2 + (2 mu) + p - 1), 
    1 - (q13 + q1 + qi1 + mu + p + (L mu)), 
    1 - (q13 + q1 + mu + p)}, {}}, {{0}, {-(q13 + q1 + mu + p)}}, 
    C1/(BI*gammak1*(Bsd + ((Sin[Pi/M]^2)/(Sin[theta]^2))))] /. params;
LogPlot[mg, {theta, 1/2, Pi/2}, PlotRange -> {{0, Pi/2}, All}, 
   ImageSize -> Large, AxesLabel -> {theta, "mg"}, LabelStyle -> {15, Bold, Black}]

is very slow, noisy, and highly inaccurate.

Not surprisingly, trying to NIntegrate this fails. In fact, obtaining a smooth curve requires a very large WorkingPrecision, especially near theta = 0.

lp = LogPlot[mg, {theta, 1/100, Pi/2}, WorkingPrecision -> 2000, 
    PlotRange -> {{0, Pi/2}, All}, ImageSize -> Large, AxesLabel -> {theta, "mg"}, 
    LabelStyle -> {15, Bold, Black}])
Cases[lp, Line[z_] -> z, Infinity][[1, 2, 1]]

(* 0.022446 *)

Even with WorkingPrecision -> 2000, Plot can reach only down to 0.022446, well above the requested lower bound, 0.01. About 70 seconds was required for this computation. (Note that the first value obtained by Plot has been discarded, because it visibly is inaccurate.) Performing the desired integral requires the integrand to be evaluated all the way to theta = 0. To obtain a reasonable estimate of the function there, Fit the very left of the curve above to a biquadratic in theta.

Rest[First[Cases[lp, Line[z_] -> z, Infinity]]];
fi = Interpolation[%];
val = %%[[1, 1]];
Fit[%%%[[;; 10]], {1, theta^2, theta^4}, theta]
s = Piecewise[{{fi[theta], theta >= val}}, %];
(* 17.877 - 25.3885 theta^2 + 159.428 theta^4 *)

mg at theta = 1/100 can be computed from

Block[{$MaxExtraPrecision = 100000}, N[mg /. theta -> 10^-2, 6]]
(* 5.79133*10^7 *)

in about 70 seconds and agrees well with Exp[s /. theta -> 1/100]. Plotting Exp[s] has a credible appearance as well. The integral now is trivial to perform.

NIntegrate[Exp[s], {theta, 0, Pi/2}]
(* 1.62895*10^7 *)

In general, the same procedure should work for other parameters.

Addendum: Faster, simpler code

Most of the code above can be replaced by

fb = Interpolation@Block[{$MaxExtraPrecision = 100000}, 
    Table[{theta, N[mg, 6]}, {theta, 2/100, Ceiling[Pi/2, 1/100], 1/100}]];
Quiet@NIntegrate[fb[theta], {theta, 0, Pi/2}]
(* 1.6289*10^7 *)

as before. This computation requires about 18 seconds, rather than 70 seconds for the previous approach.

Applying the present approach to a different set of parameters, say,

SeedRandom[1812];
params = Join[Thread[{p, q1, q2, q13, qi1} -> Rationalize[RandomReal[2, 5], 0]], 
    {mu -> 3/2, L -> 3, gammak1 -> 1, BI -> 9, C1 -> 4, Bsd -> 4, M -> 2}];

and defining mg and fb as before, but with the new params, then gives for the integral

(* 2.0744*10^6 *)

in about 30 seconds. For completeness, the correspond plot of the integrand is

Quiet@LogPlot[fb[theta], {theta, 0, Pi/2}, ImageSize -> Large, 
    AxesLabel -> {theta, "mg"}, LabelStyle -> {15, Bold, Black}]