Determining the range of parameters that yield real values for a certain NIntegrate form

Question

I have specified just one set of $s$ and $g$ values that yields a real value for the NIntegrate below. It is possible that some $s,g$ combination can give rise to an imaginary output. I would like to know the method of determining the whole range of $s,g$ values that give rise to real value outcome for the integration. A region plot might be useful, but nevertheless a neat way of doing this will be appreciated.

s = 1.1; 
g = 0.1; 

NIntegrate[
  (p^2* Sqrt[1 - (2*p)/(-g + p*(1 + p))])/(-1 + (p/g)*(p - 1)) - 
  p^2/((1 + p/g)*(1 + p))*Sqrt[1 - (2*p)/(g + p*(1 + p))]), 
  {p, s, 1000}]

Answer 1

Amplifying on Chenminqi's answer

g = 1/10;

func = (p^2*Sqrt[1 - (2*p)/(-g + p*(1 + p))])/(-1 + (p/g)*(p - 1)) - 
    p^2/((1 + p/g)*(1 + p))*Sqrt[1 - (2*p)/(g + p*(1 + p))] // Simplify;

fd = FunctionDomain[func, p]

This is equivalent to requiring that the arguments of Sqrt be positive

fd == Reduce[
  Thread[
   Cases[func, Sqrt[x_] -> x, Infinity] > 0],
  p]

True

fd // N

p < -1.09161 || -0.091608 < p < 0.091608 || p > 1.09161

Plot[func, {p, -5, 2.5},
 Frame -> True,
 Axes -> False,
 PlotPoints -> 101]