Suppose we have the $L$-loop amplitude of the form
$$\mathcal{I}_L=\int \prod_{i=1}^L \frac{d^D q_i}{(2 \pi)^D} \frac{1}{q_i^2} \frac{1}{(p-\sum_{i=1}^L q_i)^2}.$$
Introducing Feynman parameters to merge the denominators as usual, we may write the amplitude in the form
$$\mathcal{I}_L = (N-1)! \int_0^1 \prod_{j=1}^{L+1} dx_j \delta \left( \sum_{i=1}^{L+1} x_i-1\right) \int \prod_{i=1}^L \frac{d^D q_i}{(2 \pi)^D} \left[ q_i q_j M_{ij}-2 q_j K_j+J \right]^{-(L+1)}$$
where $M$ is a symmetric matrix. We can evaluate the $q$ integrals using the formula
$$\int \frac{d^D l}{(2 \pi)^D}\frac{1}{(l^2-\Delta)^N} = \frac{(-1)^Ni}{(4 \pi)^{D/2}}\frac{\Gamma(N-D/2)}{\Gamma(N)}\left( \frac{1}{ \Delta}\right)^{N-D/2},$$
It can be shown that we can write $\mathcal{I}_L$ as
$$\mathcal{I}_L=\frac{(-1)^L \Gamma(L+1-D)}{(4 \pi)^D} \int_0^1 \prod_{i=1}^{L+1} dx_i \delta \left( \sum_{i=1}^{L+1} x_i-1\right) \frac{\mathcal{U}^{L+1-3D/2}}{\mathcal{F}^{L+1-D}},$$
where we have defined $\mathcal{U}=\det M$ and $ \mathcal{F}=\det M \left( K_i M^{-1}_{ij} K_j -J \right).$
By setting $D=4-2 \epsilon$, my goal is to calculate the divergent part
$$\mathcal{I}_L = \frac{c_L}{\epsilon} + \mathcal{O}(\epsilon^0)$$
for any $L$. With that aim in mind, let's first look at $L=2$. Doing the above computations and calculating $\mathcal{U}$ and $\mathcal{F}$ explicitly, amounts to
$$c_L \propto \int dx_1 dx_2 dx_3 \delta \left( x_1+x_2+x_3-1\right) \frac{x_1 x_2 x_3}{(x_1 x_2 +x_2 x_3 + x_1 x_3)^3}.$$
The form for larger $L$ follows a similar pattern (product of $x_i$ in the numerator to some power, and sum of products of $x_i$ with one missing in each term, just like above), so evaluating the integral for $L=2$ would probably lead the way for a more general result. However, how does one treat integrals such as this? Any ideas about how one might try to evaluate it completely?