Loop correction to the 3-point vertex

In chapter 18 in Srednicki, the loop correction to the 3-point vertex in $\phi^3$ theory in six dimensions is calculated. In this post, I give comments on its numerical calculation in Mathematica.
The calculation of the 3-point vertex $V_3$ involves a two-dimensional integral over Feynman parameters [1] $$\begin{equation*} V_3(s) = g - g\ \alpha \int_0^1\hspace{-0.5em} dx_1 \int_0^{1 - x_1}\hspace{-2em} dx_2\ \log \left( D_3(s) / m^2 \right) \end{equation*}$$ with $$\begin{equation*} D_3(s) = - x_1 x_2 s + [ 1 - (x_1 + x_2) x_3 ] m^2 \end{equation*}$$ Here, $g$ is the coupling constant, $s$ is a Mandelstam variable, $m$ is the mass of the particles, $\alpha = g^2/(4 \pi)^3$ and $x_1 + x_2 + x_3 = 1$. The branch cut of the logarithm is such that $\log(z) = \log(|z| ) - i \pi$ for $z< 0$.
I implement $V_3$ in Mathematica as

V3Naive[s_,g_]:= g ( 1 - 1 /2 \[Alpha][g] 2 NIntegrate[mylog[D3[x1,x2,1 - x1-x2,s]/m^2],{x1,0.,1.},{x2,0.,1.-x1}]) 

with

mylog[z_] /; z \[Element] Reals := If[z > 0, Log[z], If[z < 0, Log[Abs[z]] - I Pi]]

When calculating the integral naively, Mathematica warns about numerical issues. I illustrate this in the following graph, which is a picture of $V_3$ as function of $s$, calculated with V3Naive. If there are numerical problems the points are plotted in red. If there are no numerical problems, the points are plotted in green. All Mathematica code can be found at the bottom of this post.

Fig 1. $V_3$ as function of $s$, calculated with V3Naive
$g=1$ and $m=1$

For $s> 4 m^2$, Mathematica warns about numerical problems. The problem originates from the fact that there is a singularity in the integrand when $D_3$ becomes zero. Mathematica indeed notices that more points are needed around the singularity. Here is a graph of the integration points.

Fig 2. Integration points chosen in V3Naive
$s = 16$, $m=1$

The documentation of NIntegrate [2] explains that splitting the integration domain along the singular curve helps the algorithm. Better code is therefore

V3Better[s_,g_]:=  g ( 1 - 1 /2 \[Alpha][g] 2  (NIntegrate[mylog[D3[x1,x2,1 - x1-x2,s]/m^2]*Boole[- x1 x2 s + ( 1 - ( x1 +x2) (1 - (x1+x2)))m^2<0],{x1,0.,1.},{x2,0.,1.-x1}] +NIntegrate[mylog[D3[x1,x2,1 - x1-x2,s]/m^2]*Boole[- x1 x2 s + ( 1 - ( x1 +x2) (1 - (x1+x2)))m^2>0],{x1,0.,1.},{x2,0.,1.-x1}] ))

In this case Mathematica allocates many, but not extremely many, points around the singular curve, see next picture

Fig 3. Integration points chosen in V3Better
The red (blue) points are where $D_3$ is negative (positive).

I now produce a similar graph as figure 1, but using implementation V3Better instead of V3Naive

Fig 4. $V_3$ as function of $s$, calculated with V3Better

There are thus fewer warnings than before, but there are still problems for some values of $s$. I think the reason might be that the singular curve has an almost horizontal and vertical asymptotic for those values of $s$. I could add those points as information in the numerical integration, however, I do not want to spend more time on this. I also read in an article [3] that ''The calculation of the integrals over the Feynman parameters is straightforward in principle but very tedious in practice'' and ''during the process of calculating a diagram the most time consuming task is the numerical evaluation of the integrals over the Feynman parameters''. I did not know this before I started to calculate $V_3$. There are even books written about the calculation of integrals over Feynman diagrams [4].

References
[1] Mark Srednicki, Quantum Field Theory, pages 124-125
[2] NIntegrate Integration Strategies
[3] J. Calmet, A. Visconti, Computing Methods in Quantum Electrodynamics, Chapter in Field Theory, Quantization and Statistical Physics Volume 6 of the series Mathematical Physics and Applied Mathematics pp 33-57
[4] Vladimir Smirnov, Feynman Integral Calculus

The code that produces figure 1
(* inelegant code *)
tt =Table[s,{s,-50,100,1}];
For [i=1,i<= Length[tt],i++,
s = tt[[i]];
result= Check[tmp = V3Naive[s,1];"nowarning","warning"];
tt[[i]] = If[result ===  "nowarning", {s,Re[tmp],Green},{s,Re[tmp],Red}]
] 

Graphics[{#3, PointSize[0.01], Point[{#1, #2}]} & @@@ tt, Axes -> True, AspectRatio -> Full, PlotRange -> All]

The last bit of code was inspired by a question on Mathematica Stack Exchange.

The code that produces figure 2
I found the following code at the documentation of NIntegrate [2] ss = (2 2)^2
gr={Red,Point/@N@Reap[NIntegrate[mylog[D3[x1,x2,1 - x1-x2,ss]/m^2],{x1,0.,1.},{x2,0.,1.-x1},EvaluationMonitor:>Sow[{x1,x2}]]][[2,1]]}; 

Graphics[{PointSize[0.006],gr},Axes->True,AxesOrigin->{0,0}] 

The code that produces figure 3
I found this code again at [2].
gr1={Red,Point/@N@Reap[NIntegrate[mylog[D3[x1,x2,1 - x1-x2,ss]/m^2]*Boole[- x1 x2 ss + ( 1 - ( x1 +x2) (1 - (x1+x2)))m^2<0],{x1,0.,1.},{x2,0.,1.-x1},EvaluationMonitor:>Sow[{x1,x2}]]][[2,1]]};

gr2={Blue,Point/@N@Reap[NIntegrate[mylog[D3[x1,x2,1 - x1-x2,ss]/m^2]*Boole[- x1 x2 ss + ( 1 - ( x1 +x2) (1 - (x1+x2)))m^2>0],{x1,0.,1.},{x2,0.,1.-x1},EvaluationMonitor:>Sow[{x1,x2}]]][[2,1]]};

Graphics[{PointSize[0.006],gr1,gr2},Axes->True,AxesOrigin->{0,0}]