UPDATE:
A generalization of this integral appears on page 429 of the textbook A Treatise on the Theory of Bessel Functions, namely
$$\small \int_{0}^{\infty} \frac{x}{x^{2}+\alpha^{2}} \, J_{\mu} (\beta x) \left(\cos \left(\frac{\pi(\mu-\nu)}{2} \right)J_{\nu}(rx) + \sin \left(\frac{\pi(\mu-\nu)}{2} \right)Y_{\nu}(rx) \right) \, \mathrm dx = I_{\mu}(\beta \alpha) K_{\nu}(r \alpha) $$ where $r \ge \beta >0$, and $\mu$ and $\nu$ are nonnegative real parameters such that $\mu > \nu -2$.
As in my previous answer, we can exploit properties of the Hankel function of the first kind.
Let $H_{0}^{(1)}(z)$ be the Hankel function of first kind of order zero defined as $$H_{0}^{(1)}(z) = J_{0}(z) + i Y_{0}(z), $$ where $Y_{0}(z)$ is the Bessel function of the second kind of order zero.
The principal branch of $H_{0}^{(1)}(z)$ has a branch cut on the negative real axis.
Since $$Y_{0}(xe^{i \pi})= Y_{0}(x) + 2i J_{0}(x), \quad x >0,$$ it follows that
$$H_{0}^{(1)}(xe^{i \pi}) = - J_{0}(x) + iY_{0}(x), \quad x >0. $$
And since $\frac{x}{x^{2}+\alpha^{2}}$ is an odd function, we have $$\int_{0}^{\infty} \frac{x}{x^{2}+\alpha^{2}} \, J_{0}(x) J_{0}(rx) \, \mathrm dx = \frac{1}{2} \, \Re \int_{-\infty}^{\infty} \frac{x}{x^{2}+\alpha^{2}} \, H_{0}^{(1)} (x) J_{0}(rx) \, \mathrm dx,$$ where the integration form $-\infty$ to $0$ is done on the upper side of the branch cut.
Now let's integrate the function $$f(z) = \frac{z}{z^{2}+\alpha^{2}} \, H_{0}^{(1)}(z) J_{0}(rz), \quad 0 < r \le 1 \, , $$ around a contour consisting of the upper side of the branch cut from $-R$ to $-\epsilon$, a small clockwise- oriented semicircle about the origin, the real axis from $\epsilon$ to $R$, and the upper half of the circle $|z|=R$.
Since $\lim_{z \to 0} f(z) = 0$, the contribution from the small semicircle about the origin vanishes as $\epsilon \to 0$.
And as $|z| \to \infty$ in the upper half plane, $|f(z)|$ is asymptotic to $$\frac{1}{\pi \sqrt{r}} \frac{e^{(r-1)\Im(z)}}{|z|^{2}}. $$
(See here.)
Since $0 < r \le 1$, the integral vanishes on the upper half of the circle $|z|=R$ as $R \to \infty$ (by the estimation lemma), and we have $$ \begin{align} \int_{-\infty}^{\infty} \frac{x}{x^{2}+\alpha^{2}} \, H_{0}^{(1)} (x) J_{0}(rx) \, \mathrm dx &= 2 \pi i \operatorname*{Res}_{z=i \alpha} f(z) \\ &= 2 \pi i \lim_{z \to i \alpha} \frac{z}{z+i \alpha} \, H_{0}^{(1)}(z) J_{0}(rz) \\ &= i \pi \, H_{0}^{(1)}(i \alpha)J_{0}(i r \alpha) \\ &\overset{(\spadesuit)}{=} i \pi \left(\frac{2K_{0}(\alpha)}{i \pi} \right) I_{0}(r \alpha) \\&= 2 K_{0}(\alpha) I_{0}(r \alpha). \end{align}$$
Equating the real parts on both sides of the equation, we have $$\int_{0}^{\infty} \frac{x}{x^{2}+ \alpha^{2}} \, J_{0}(x) J_{0}(rx) \, \mathrm dx = K_{0}(\alpha) I_{0}(r \alpha), \quad 0 < r \le 1. $$
This result also holds for $r=0$.
$\spadesuit$ https://dlmf.nist.gov/10.27#E8
To show that $$\int_{0}^{\infty} \frac{x}{x^{2}+ \alpha^{2}} \, J_{0}(x) J_{0}(rx) \, \mathrm dx = I_{0}(\alpha) K_{0}(r \alpha), \quad r \ge 1, $$ integrate $$g(z) = \frac{z}{z^{2}+\alpha^{2}} \, J_{0}(z) H_{0}^{(1)}(rz) $$ around the same contour.
Best Answer
Since given function is even, we can just evaluate the integral from 0 to infinity and double it.
Let $\ t = \frac{1}{1+x^2}$.
Then, $\ dx = \frac{-dt}{2t^{\frac{3}{2}}(1-t)^{\frac{1}{2}}}$
Substituting $\ t$ throughout, we get:
$$\ I = \int_{\ 1}^{0} t^3 \cdot \frac{-dt}{t^{\frac{3}{2}}(1-t)^{\frac{1}{2}}}$$
$$ = \int_{\ 0}^{\ 1} t^{\frac{3}{2}}(1-t)^{\frac{-1}{2}}dt$$ $$ = \ B(\frac{5}{2},\frac{1}{2})$$
Where $\ B(x,y)$ is the Beta function.
Evaluating, we get
$$\ I = \frac{\Gamma(\frac{5}{2})\Gamma(\frac{1}{2})}{\Gamma(3)}$$
$$\ I = \frac{\frac{3\sqrt{\pi}}{4} \cdot \sqrt{\pi}}{2}$$
$$\ I = \frac{3\pi}{8}$$