Evaluating 4-dim'l Euclidean integral using spherical coordinates Scenario: I have an integral originally expressed in Cartesian coordinates which should--in principle--converge, but I have difficulty evaluating it explicitly (due to apparent divergences) in spherical coordinates using the residue theorem.
Setup: I want to evaluate the following 4-dimensional integral in Euclidean space, 
\begin{equation}
\mathcal{I}_L(x_1,x_2)=\int_{|y-x_2|\le L} d^4y \frac{1}{(y-x_1)^2 (y-x_2)^2}, \tag{1}\label{1}
\end{equation}
where $y,x_1,x_2$ are 4-vectors; "$||$" denotes the Euclidean distance function;  and $(y-x_1)^2\equiv(y-x_1)\cdot(y-x_1)$ with "$\cdot$" the dot product.  Also assume $L>0$ but finite.
First note that although the integrand diverges at $y=x_1,x_2$, the integral should be expected to converge on dimensional grounds.
Also note that shifting the integration variable, $y\to y+x_2$, in \eqref{1} gives,
\begin{align}
\mathcal{I}_L(x_1,x_2)=\int_{|y|\le L} d^4y \frac{1}{(y-(x_1-x_2))^2 y^2} \tag{2}\label{2}\\
\end{align}
Now this 4-dim'l integral can be simplified to a 2-dim'l integral by switching to 4-dim'l "spherical coordinates", 
\begin{align}
y^1&=|y|\cos(\phi_1)\\
y^2&=|y|\sin(\phi_1)\cos(\phi_2)\\
y^3&=|y|\sin(\phi_1)\sin(\phi_2)\cos(\phi_3)\\
y^4&=|y|\sin(\phi_1)\sin(\phi_2)\sin(\phi_3),\\
\end{align}
where $\phi_1,\phi_2$  range over $[0,\pi]$ while $\phi_3$ ranges over $[0,2\pi]$. [Notation: superscripts denote components of Cartesian 4-vector--i.e. $y=(y^1,y^2,y^3,y^4)$] 
Choosing to orient our axes such that $\phi_3$ coincides with the angle between the 4-vectors $y$ and $x_1-x_2$ allows us to write \eqref{2} as 
\begin{align}
\mathcal{I}_L(\sigma)&= \int_0^L d|y|\,|y|^3 \int_{S^3} d\Omega_3 \frac{1}{(|y|^2-2|y|\sigma \cos{\phi_3}+\sigma^2)|y|^2}\\
&=\pi \int_0^L d|y| \int_0^{2\pi} d\phi_3 \frac{|y|}{|y|^2-2|y|\sigma \cos{\phi_3}+\sigma^2}, \label{3}\tag{3}\\
\end{align}
where I define the shorthand $\sigma\equiv |x_1-x_2|$.
Now, the angular part of the integral in \eqref{3} looks like a good candidate for a calculus of residues approach, so I take $z=e^{i\phi}$ and find
\begin{align}
\mathcal{I}_L(\sigma)&=\pi i \int_0^L d|y| \oint_{|z|=1}dz \frac{|y| }{(|y|z-\sigma)(\sigma z-|y|)}\\
&=\frac{\pi i}{\sigma}\int_0^L d|y|\oint_{|z|=1}dz\frac{1}{(z-\sigma/|y|)(z-|y|/\sigma)}.\label{4}\tag{4}\\
\end{align}
Using the residue theorem, I find
\begin{align}
\oint_{|z|=1}dz\frac{dz}{(z-\sigma/|y|)(z-|y|/\sigma)}=2\pi i \frac{\sigma |y|}{|y|^2-\sigma^2}\left\{ \theta (\sigma-|y|)-\theta (|y|-\sigma)\right\}, \tag{5}\label{5}
\end{align}
where $\theta(x>0)=1$ and $\theta(x)=0$ otherwise.
Plugging \eqref{5} into \eqref{4}, yields
\begin{align}
\mathcal{I}_L(\sigma)=-2\pi^2\int_0^L d|y|\frac{|y|}{|y|^2-\sigma^2}\left\{ \theta (\sigma-|y|)-\theta (|y|-\sigma)\right\} \label{*} \tag{*}
\end{align}
Confusion/Contradiction: For $L<\sigma$, \eqref{*} does seem to give a convergent result:
\begin{align}
\mathcal{I}_{L<\sigma}(\sigma)=-\pi^2 \ln\left[1-L^2/\sigma^2\right].\\
\end{align}
[Recall: I have defined the shorthand $\sigma\equiv |x_1-x_2|$] 
However, for $L\ge\sigma$, $\eqref{*}$ seems to diverge which contradicts the obvious convergence of $\eqref{1}$.
More explicitly, for $L\ge\sigma$, the integral \eqref{*} naturally splits into 2 parts:
\begin{align}
\mathcal{I}_L(\sigma)&=-2\pi^2\left[ \int_0^\sigma-\int_\sigma^L\right] d|y|\frac{|y|}{|y|^2-\sigma^2}\\
&=-\pi^2 \left[ \ln{(|y|^2-\sigma^2)}\right]^\sigma_0+\pi^2\left[\ln{(|y|^2-\sigma^2)}\right]^L_\sigma\\
&\sim -\infty\\
\end{align}
Edit: (elaboration on H.H.Rugh's answer)
As correctly pointed out by H.H.Rugh below, $\phi_3$ cannot be chosen to coincide with the angle between $y$ and $x_1-x_2$.  The nature of my mistake can be easily visualized in 3-dimensional Euclidean space, where it is conventional to choose one angle $\phi:[0,2\pi]$ to denote rotations in the "$x,y$"-plane around the "$z$-axis" and one angle $\theta:[0,\pi]$ to denote rotations away from the $z$-axis as follows:
\begin{align}
x&=r \sin\theta \cos\phi\\
y&=r \sin\theta \sin\phi\\
z&=r \cos\theta,\\
\end{align}
where $r^2=x^2+y^2+z^2$.
Clearly the $\theta$-component of an arbitrary 3-vector in spherical coordinates, $(r,\phi,\theta)$ coincides with the angle between that 3-vector and one lying along the $z$-axis.  However, the $\phi$-component of that same arbitrary 3-vector does not necessarily coincide with the angle between that 3-vector and one lying along the $x$ or $y$-axis.  This is analogous to my mistake of choosing $\phi_3$ to coincide with the angle between the 4-vectors $y$ and $x_1-x_2$.
For completeness, I also explicitly finish the computation started by H.H.Rugh, beginning with 
\begin{align}
\mathcal{I}_L(\sigma)&= 2\pi \int_0^L dr \int_0^{2\pi} d\phi  \frac{r \sin^2 \phi}{r^2 - 2 r \sigma \cos \phi_1 + \sigma^2 } \\
&=\frac{i\pi}{2\sigma}\int_0^L dr \oint_{z=1} dz \frac{(z^2-1)^2}{z^2(z-\sigma/r)(z-r/\sigma)}\\
\end{align}
There is now an additional relevant residue at $z=0$.  Applying the residue theorem, I find
$$
\mathcal{I}_{L}(\sigma)=
\begin{cases}
\pi^2 L^2/\sigma^2 & \text{if } L<\sigma  \\
\pi^2 \left(1+\ln(L^2/\sigma)\right) & \text{if } L\ge \sigma\\
\end{cases}
$$
 A: I think that the divergence comes from a wrong setup for the spherical coordinates relative to $\Delta x=x_1-x_2$. I don't see how you can make  $\phi_3$ the angle between $y$ and $\Delta x$? Instead it seems natural to let 
$\Delta x$ be in the (positive) $y^1$ direction so the angle is $\phi_1$.
The volume element in 4-dim spherical coordinates is $$r^3 \sin^2(\phi_1) \sin(\phi_2) dr \; d\phi_1\; d\phi_2\; d\phi_3$$ and the square distance between $y$ and $\Delta x$ is (with $\sigma=|\Delta x|$ and the described choice of angle)
$$ (r \cos \phi_1 - \sigma )^2 + (r \sin \phi_1)^2 = r^2 - 2 r \sigma \cos \phi_1 + \sigma^2 $$ 
For $r=\sigma$ this in the denominator indeed gives rise to a singularity  as $\phi_1$ goes to zero of order $1/\phi_1^2$ but
when you integrate you have now (due to dim 4) the important factor $(\sin \phi_1)^2$ which will remove the singularity. The integral after carrying out the trivial part should look something like (modulo errors in my calculations)
$$ 4\pi \int_0^\pi dr \int_0^\pi d\phi  \frac{r \sin^2 \phi}{(r-\sigma)^2+ 4r\sigma (1-\cos\phi)} $$
You may replace by half the integral over $[0,2\pi]$ and try to evaluate using residues but in any case the integrals are  now integrals of bounded functions. 
