Computing $\int_{-\infty}^{\infty} \frac{\cos x}{x^{2} + a^{2}}dx$ using residue calculus
Let $\gamma$ be the path along the real axis then circling back counter-clockwise through the upper half-plane, letting the circle get infinitely big.
$$
\begin{align}
\int_{-\infty}^\infty\frac{\cos(x)}{x^2+a^2}\mathrm{d}x\tag{1}
&=\Re\left(\int_{-\infty}^\infty\frac{\exp(ix)}{x^2+a^2}\mathrm{d}x\right)\\\tag{2}
&=\Re\left(\int_{\gamma}\frac{\exp(ix)}{x^2+a^2}\mathrm{d}x\right)\\\tag{3}
&=\Re\left(2\pi i\,\mathrm{Res}\left(\frac{\exp(ix)}{x^2+a^2},ia\right)\right)\\\tag{4}
&=\Re\left(2\pi i\,\lim_{z\to ia}\frac{\exp(ix)}{x+ia}\right)\\\tag{5}
&=\Re\left(2\pi i\,\frac{\exp(-a)}{2ia}\right)\\[3pt]
&=\frac{\pi \exp(-a)}{a}\tag6
\end{align}
$$
$(1)$: $\Re(\exp(ix))=\cos(x)$
$(2)$: The integral along the circle back through the upper half-plane vanishes as the circle gets bigger.
$(3)$: There is only one singularity of $\dfrac{\exp(ix)}{x^2+a^2}$ in the upper half-plane, at $x=ia$. The integral along $\gamma$ is the residue of $\dfrac{\exp(ix)}{x^2+a^2}$ at $x=ia$.
$(4)$: The singularity of $\dfrac{\exp(ix)}{x^2+a^2}$ at $x=ia$ is a simple pole. We can compute the residue as $\displaystyle\lim_{x\to ia}(x-ia)\frac{\exp(ix)}{x^2+a^2}=\lim_{x\to ia}\frac{\exp(ix)}{x+ia}$.
$(5)$: plug in $x=ia$.
$(6)$: evaluate
Following the same strategy, $$ \begin{align} \int_{-\infty}^\infty\frac{x\sin(x)}{x^2+a^2}\mathrm{d}x &=\Im\left(\int_{-\infty}^\infty\frac{x\exp(ix)}{x^2+a^2}\mathrm{d}x\right)\\ &=\Im\left(\int_{\gamma}\frac{x\exp(ix)}{x^2+a^2}\mathrm{d}x\right)\\ &=\Im\left(2\pi i\,\mathrm{Res}\left(\frac{x\exp(ix)}{x^2+a^2},ia\right)\right)\\ &=\Im\left(2\pi i\,\lim_{z\to ia}\frac{x\exp(ix)}{x+ia}\right)\\ &=\Im\left(2\pi i\,\frac{ia\exp(-a)}{2ia}\right)\\[6pt] &=\pi \exp(-a) \end{align} $$
Although the OP is searching for a way forward using contour integration and the residue theorem, I thought it might be instructive to present an approach that uses real analysis only. To that end, we proceed.
Let $g(a)$ be given by the convergent improper integral
$$g(a)=\int_{-\infty}^\infty \frac{\cos(x)}{x^2+a^2}\,dx \tag1$$
Exploiting the even symmetry of the integrand and enforcing the substitution $x\to ax$ reveals
$$g(a)=\frac2a\int_0^\infty \frac{\cos(ax)}{x^2+1}\,dx$$
Now let $f(a)=\frac a2 g(a)=\int_0^\infty \frac{\cos(ax)}{x^2+1}\,dx$
Since the integral $\int_0^\infty \frac{x\sin(ax)}{x^2+1}\,dx$ is uniformly convergent for $|a|\ge \delta>0$, we may differentiate under the integral in $(3)$ for $|a|>\delta>0$ to obtain
$$\begin{align} f'(a)&=-\int_0^\infty \frac{x\sin(ax)}{x^2+1}\,dx\\\\ &=-\int_0^\infty \frac{(x^2+1-1)\sin(ax)}{x(x^2+1)}\,dx\\\\ &=-\int_0^\infty \frac{\sin(ax)}{x}\,dx+\int_0^\infty \frac{\sin(ax)}{x(x^2+1)}\,dx\\\\ &=-\frac{\pi}{2}+\int_0^\infty \frac{\sin(ax)}{x(x^2+1)}\,dx\tag4 \end{align}$$
Again, since the integral $\int_0^\infty \frac{\cos(ax)}{x^2+1}\,dx$ converges uniformly for all $a$, we may differentiate under the integral in $(4)$ to obtain
$$f''(a)=\int_0^\infty \frac{\cos(ax)}{x^2+1}\,dx=f(a)\tag 5$$
Solving the second-order ODE in $(5)$ reveals
$$f(a)=C_1 e^{a}+C_2 e^{-a}$$
Using $f(0)=\pi/2$ and $f'(0)=-\pi/2$, we find that $C_1=0$ and $C_2=\frac{\pi}{2}$ and hence $f(a)=\frac{\pi e^{-a}}{2}$. Multiplying by $2/a$ yields the coveted result
$$\bbox[5px,border:2px solid #C0A000]{\int_{-\infty}^\infty \frac{\cos(x)}{x^2+a^2}\,dx=\frac{\pi}{ae^a}}$$
as expected!
Let us integrate along the contour $\Gamma$, which is the semicircle of radius $R$ described above. Let $C_R$ be the arc of the contour with radius $R$. We may say that
$$\oint_{\Gamma}\frac{\cos z}{z^2+a^2}\, dz=\int_{-R}^{R}\frac{\cos x}{x^2+a^2}\, dx+\int_{C_R}\frac{\cos z}{z^2+a^2}\, dz$$
Note that $\cos x = \operatorname{Re}\,(e^{ix})$. Thus
$$ f(z)=\frac{\operatorname{Re}\,(e^{iz})}{z^2+a^2} = \operatorname{Re}\,\left(\frac{e^{iz}}{(z+ia)(z-ia)}\right) $$
Because $f(z)$ can be written ias $g(z)e^{iaz}$ and suffices all of the conditions of Jordan's lemma, we see that the integral along the arc tends to zero when $R\to \infty$. Thus
$$ \lim_{R\to\infty}\oint_{\Gamma}\frac{\cos z}{z^2+a^2}\, dz=\lim_{R\to\infty} \int_{-R}^{R}\frac{\cos x}{x^2+a^2}\, dx+0 $$
To solve the LHS, we find the residues. The residue, similar to the one you found is $$ b=\frac{e^{i^2a}}{ia+ia}=\frac{e^{-a}}{2ia} $$
Thus
$$ \displaystyle\int_{-\infty}^{\infty} \frac{\cos x}{x^{2} + a^{2}}\, dx= \lim_{R\to\infty}\oint_{\Gamma}\frac{\cos z}{z^2+a^2}\, dz=\operatorname{Re}\,\left(\displaystyle\int_{-\infty}^{\infty} \frac{e^{iz}}{z^{2} + a^{2}}\, dz\right) =\operatorname{Re}\,(2\pi ib)=\operatorname{Re}\,(2\pi i\frac{e^{-a}}{2ia}) = \frac{\pi e^{-a}}{a} $$
Similarly, with $\displaystyle\int_{-\infty}^{\infty} \frac{x \sin x}{x^{2} + a^{2}}\ dx$ we change $\sin x$ to $\operatorname{Im}\,(e^{ix})$ and make the function complex-valued.
$$g(z)=\frac{z\sin z}{z^2+a^2}=\operatorname{Im}\,\left(\frac{ze^{iz}}{(z+ia)(z-ia)}\right)$$
Jordan's lemma again can be used to show that the integral around the arc tends to zero as $R\to\infty$. We then proceed to find the residues. The residue of pole $z_1=ia$ is
$$b=\frac{iae^{-a}}{2ia}=\frac{e^{-a}}{2}$$
$z_1$ is the unique pole in the contour, so upon multiplying its residue, $b$, with $2\pi i$ we find: $$ \displaystyle\int_{-\infty}^{\infty} \frac{x \sin x}{x^{2} + a^{2}}\, dx= \lim_{R\to\infty}\oint_{\Gamma} \frac{x \sin x}{x^{2} + a^{2}}\, dz= \operatorname{Im}\,\left(\displaystyle\int_{-\infty}^{\infty} \frac{ze^{iz}}{(z+ia)(z-ia)}\, dz\right)= \operatorname{Im}\,(2\pi ib)= \operatorname{Im}\,(2\pi i \frac{e^{-a}}{2})= \operatorname{Im}\,(\pi e^{-a}i)= \pi e^{-a} $$