cosϕSol[θS_, sinθh_, cosθh_, a_] =
  cosϕ /. Solve[cosϕ*Sin[θS] sinθh + Cos[θS]*cosθh == a, cosϕ][[1]];
θhSol[θS_, a_] =
  θh /. Solve[cosϕSol[θS, Sin[θh], Cos[θh], a] == 1, θh] /. {C[1] -> 0};
A[ES_, mS_] = Sqrt[4*ES^2 - 125^2]/Sqrt[4*ES^2 - 4*mS^2];
integral1[ES_, θS_, mS_] :=
 NIntegrate[1,
  {θh,
   Max[Min[θhSol[θS, A[ES, mS]][[1]], θhSol[θS, A[ES, mS]][[2]]], 0],
   Min[Max[θhSol[θS, A[ES, mS]][[2]], θhSol[θS, A[ES, mS]][[1]]], Pi]},
  {ϕh, 0,
   Min[Pi, ArcCos[cosϕSol[θS, Sin[θh], Cos[θh], A[ES, mS]]]]}]

Clear[cosϕSol, θhSol, A, integral1]

cosϕSol[θS_, sinθh_, cosθh_, a_] :=
  cosϕ /.
   Solve[cosϕ*Sin[θS] sinθh +
       Cos[θS]*cosθh == a, cosϕ][[1]];

θhSol[θS_,
  a_] := θh /.
   Solve[cosϕSol[θS, Sin[θh], Cos[θh], a] ==
      1, θh] /. {C[1] -> 0}

A[ES_, mS_] := Sqrt[4*ES^2 - 125^2]/Sqrt[4*ES^2 - 4*mS^2];

integral1[ES_, θS_, mS_] :=
 NIntegrate[1,
  {θh,
   Max[Min[θhSol[θS,
       A[ES, mS]][[1]], θhSol[θS, A[ES, mS]][[2]]], 0],
   Min[Max[θhSol[θS,
       A[ES, mS]][[2]], θhSol[θS, A[ES, mS]][[1]]], Pi]},
  {ϕh, 0,
   Min[Pi, ArcCos[
     cosϕSol[θS, Sin[θh], Cos[θh],
      A[ES, mS]]]]},
  Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}]

integral1[1500, 0.01, 40]

(* During evaluation of In[132]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[132]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[132]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[132]:= General::stop: Further output of Solve::ifun will be suppressed during this calculation.

During evaluation of In[132]:= NIntegrate::nlim: ϕh = Min[3.14159,3.14159 -0.545349 I] is not a valid limit of integration. *)

(* NIntegrate[1, {θh,
  Max[Min[θhSol[0.01, A[1500, 40]][[
     1]], θhSol[0.01, A[1500, 40]][[2]]], 0],
  Min[Max[θhSol[0.01, A[1500, 40]][[
     2]], θhSol[0.01, A[1500, 40]][[1]]], π]}, {ϕh, 0,
   Min[π,
   ArcCos[cosϕSol[0.01, Sin[θh], Cos[θh],
     A[1500, 40]]]]},
 Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}] *)

integral1[ES_, θS_, mS_] :=
 NIntegrate[1,
  {θh,
   Max[Min[θhSol[θS,
       A[ES, mS]][[1]], θhSol[θS, A[ES, mS]][[2]]], 0],
   Min[Max[θhSol[θS,
       A[ES, mS]][[2]], θhSol[θS, A[ES, mS]][[1]]], Pi]},
  {ϕh, 0,
   Min[Pi, Re@
     ArcCos[cosϕSol[θS, Sin[θh], Cos[θh],
       A[ES, mS]]]]},
  Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}]

integral1[1500, 0.01, 40]

(* During evaluation of In[134]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[134]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[134]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[134]:= General::stop: Further output of Solve::ifun will be suppressed during this calculation. *)

(* 0.0981353 *)

integral1[ES_, θS_, mS_] :=
  NIntegrate[1,
   {θh,
    Max[Min[θhSol[θS, A[ES, mS]][[1]], θhSol[θS, A[ES, mS]][[2]]], 0],
    Min[Max[θhSol[θS, A[ES, mS]][[2]], θhSol[θS, A[ES, mS]][[1]]], Pi]},
   {ϕh, 0,
    Min[Pi, ArcCos[Piecewise[{{#, -1 <= # <= 1}}, 1.] &@
       cosϕSol[θS, Sin[θh], Cos[θh], A[ES, mS]]]
     ]}];

integral1[1500, 0.01, 40]

(*  0.0289182  *)

integral1[ES_, θS_, mS_] :=
  With[{θh0 = Sort@Clip[θhSol[θS, A[ES, mS]], {0, Pi}]},
   NIntegrate[1,
    {θh, θh0[[1]], θh0[[2]]},
    {ϕh,0,
     ArcCos[Piecewise[{{#, -1 <= # <= 1}}, 1.] &@
       cosϕSol[θS, Sin[θh], Cos[θh], A[ES, mS]]
      ]},
    MinRecursion -> 2]
   ];

integral1[1500, 0.01, 40]
(*  0.0289182  *)

ClearAll[realACos];  (* another way to code a real arc cosine *)
realACos[x_?NumericQ /; -1 <= x <= 1] := ArcCos[x];
realACos[_?NumericQ] := 0;

integral2[ES_, θS_, mS_] :=
  With[{θh0 = Sort@Clip[θhSol[θS, A[ES, mS]], {0, Pi}]},
   NIntegrate[
    realACos[cosϕSol[θS, Sin[θh], Cos[θh], A[ES, mS]]],
    {θh, θh0[[1]], θh0[[2]]},
    MinRecursion -> 5, MaxRecursion -> 20,
    Method -> "GaussKronrodRule"]
   ];

integral2[1500, 0.01, 40]
(*  0.0289182  *)