2D Integral of inverse of product of absolute values. I am interested in evaluating integrals of the following form:
$$ \int_\mathbb{C}\frac{1}{|z-z_1| \dots |z-z_6|}\mathrm|{d}z|^2$$
where $|\mathrm{d}z|^2$ is the standard area form on $\mathbb C$ and the $z_j$ are pairwise distinct complex numbers. As the integrand is $o(|z|^{-6})$ for $ |z| \to \infty$ and $o(|z-z_j|^{-1})$ for $z \to z_j$, we know that this integral is finite.
One may rewrite it as a double integral as
$$\int_\mathbb{R} \int_\mathbb{R} \frac{1}{\sqrt{((x-x_1)^2+(y-y_1)^2)\dots ((x-x_6)^2+(y-y_6)^2)}}\mathrm{d}x\mathrm{d}y,$$
where $z_j = x_j + y_j\mathrm{i}$.
Is there any literature on these types of integrals or any elementary approach to calculate these? I've tried doing it numerically (for $z_j = j$) but got no results. For context, these integrals arise when computing the special Kähler structure on the moduli space of Higgs-bundles on a five-punctured $\mathbb{CP}^1$. So the $z_1, \dots, z_5$ are the punctures, while the $|z-z_6|$-term stems from a point in the moduli space. In fact, using a Mobius transform, we may assume that $z_1,z_2, z_3$ are fixed, for example to $0,1,2$ or whatever makes the evaluation of the integral easiest.

Feynman parametrization allows us to write:
$$\frac{1}{\sqrt{A_{1}(x,y)\cdots A_{n}(x,y)}}=\frac{\Gamma(\frac{n}{2})}{\Gamma\left(\frac{1}{2}\right)^{n}}\int_{0}^{1}\prod(du_{i}u_{i}^{-\frac{1}{2}})\frac{\delta(1-\sum u_{i})}{(\sum u_{i}A_{i})^{\frac{n}{2}}}$$
with
$A_{i}(x,y)=(x-x_{i})^{2}+(y-y_{i})^{2}=x^{2}+y^{2}-2xx_{i}-2yy_{i}+x_{i}^{2}+y_{i}^{2}$
expanding the denominator and noting that $\left(\sum u_{i}\right)=1$ due to the constraint in the integration:
$$\sum u_{i}A_{i}=\left(\sum u_{i}\right)(x^{2}+y^{2})-2\left(\sum u_{i}x_{i}\right)x-2\left(\sum u_{i}y_{i}\right)y+\sum u_{i}(x_{i}^{2}+y_{i}^{2})$$ $$\equiv (x^{2}+y^{2})-X(u)x-Y(u)y+Z(u)$$
so the $x,y$ integration is over
$$1\over ((x^{2}+y^{2})-X(u)x-Y(u)y+Z(u))^{n\over2} .$$
Since the integration is over all space, we can shift $(x^2-Xx)=(x-{1\over 2} X)^2 - {X^2\over 4}$ etc to write
$$1\over ((x'^{2}+y'^{2})+{4Z(u)-X(u)^2-Y(u)^2 \over 4})^{n\over2} .$$
The $x',y'$ integral can be massaged into the form (up to prefactors):
$$\int \int{dx' dy' \over ((x'^{2}+y'^{2})+{4Z(u)-X(u)^2-Y(u)^2 \over 4})^{n\over2}}\rightarrow \int_0^{\infty} {dp \over (p+{4Z(u)-X(u)^2-Y(u)^2 \over 4})^{n\over2}}$$ and eventually
$$({4Z(u)-X(u)^2-Y(u)^2 \over 4})^{1-{n\over2}}\int_0^{\infty} {dp \over (p + 1)^{n\over2}}$$
The $p$ integral again is just a prefactor and we are left with
a power of $4Z(u)-X(u)^2-Y(u)^2$, however this final, multiple integral seems just as involved:
$$ \int_{0}^{1}\prod(du_{i}u_{i}^{-\frac{1}{2}})\delta(1-\sum u_{i})({4Z(u)-X(u)^2-Y(u)^2 \over 4})^{1-{n\over2}}$$
