Not sure what is the strange behavior of Plot and NIntegrate. But if you rescale your function to simplify:

 res[y_, g_] = f/(40000000 \[Pi]^2 ) /. {x -> y/40000000, \[Gamma] -> g/40000000} // FullSimplify;
 res[y, g] // TraditionalForm

and find

y0[g_] = Limit[res[y, g], y -> 0] // FullSimplify

it looks like your function is not divergent at y = 0 in the whole range of g:

Plot[y0[g], {g, 0, 100}, PlotTheme -> "Detailed"]

I could be wrong though, did not have much time to dig in.

I think this is due to some error / precision loss. You have a complicated function and need to first analyze it a bit analytically, maybe take it apart on simpler components to understand larger behavior. Then you could approach numerical integration. But nearest future I will not have time to play with this.