For $\alpha, \beta \geq 0$ and $\gamma > 0$ let
$$ f_{\alpha,\beta,\gamma} (z) = \frac{\exp \left(\mathrm{i} z \frac{z^2 - \alpha^2}{z^2 - \beta^2}\right)}{z^2 + \gamma^2} \, , \, z \in \mathbb{C} \setminus \{\pm\beta, \pm \mathrm{i} \gamma\} \, . $$
By symmetry we have
$$ I (\alpha,\beta,\gamma) = \frac{1}{2} \int \limits_{-\infty}^\infty f_{\alpha,\beta,\gamma} (x) \, \mathrm{d} x \, .$$
Clearly, the integral exists for every combination of parameters and is bounded by $\frac{\pi}{2 \gamma}$ .
We can compute its value using the residue theorem. First note that we have
$$ \operatorname{Res} (f_{\alpha,\beta,\gamma}, \mathrm{i} \gamma) = \frac{1}{2 \mathrm{i} \gamma} \exp \left(- \gamma \frac{\alpha^2 + \gamma^2}{\beta^2 + \gamma^2}\right) $$
and (by Jordan's lemma)
$$ \lim_{R \to \infty} \int \limits_{\Gamma_R} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z = 0 \, ,$$
where $\Gamma_R$ is a semi-circle of radius $R > 0$ in the upper half-plane. A naive application of the residue theorem therefore yields
$$ I (\alpha,\beta,\gamma) = \frac{1}{2} \left[ 2 \pi \mathrm{i} \operatorname{Res} (f_{\alpha,\beta,\gamma}, \mathrm{i} \gamma) - \lim_{R \to \infty} \int \limits_{\Gamma_R} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z \right]= \frac{\pi}{2 \gamma} \exp \left(- \gamma \frac{\alpha^2 + \gamma^2}{\beta^2 + \gamma^2}\right) \, .$$
While the result is indeed correct, this calculation is only valid for $\alpha = \beta$ . In this case we obtain (as already mentioned in the question) $I(\alpha,\alpha,\gamma) = \frac{\pi}{2 \gamma}\mathrm{e}^{- \gamma}$ .
If $\alpha \neq \beta$ , $f_{\alpha,\beta,\gamma}$ has essential singularities on the real axis, namely at $\pm \beta$ . Thus we need to deform our contour using small semi-circles $\gamma_\varepsilon (\pm \beta)$ of radius $\varepsilon > 0$ centred at these points and show that their contribution to the integral vanishes as $\varepsilon \to 0$ .
For simplicity, I will only discuss the case $\beta = 0$ in detail. We want to find
$$ \int \limits_{\gamma_\varepsilon (0)} f_{\alpha,0,\gamma} (z) \, \mathrm{d} z = \varepsilon \int \limits_0^\pi \frac{\mathrm{e}^{\mathrm{i} \left[\phi + \varepsilon \mathrm{e}^{\mathrm{i}\phi} - \alpha^2 \varepsilon^{-1} \mathrm{e}^{-\mathrm{i} \phi}\right]}}{\gamma^2+\varepsilon^2 \mathrm{e}^{2 \mathrm{i} \phi}} \, \mathrm{d} \phi = \frac{\varepsilon}{\gamma^2} \int \limits_0^\pi \mathrm{e}^{\mathrm{i} \left[\phi - \alpha^2 \varepsilon^{-1} \mathrm{e}^{-\mathrm{i} \phi}\right]} \left[1 + \mathcal{O} (\varepsilon) \right] \, \mathrm{d} \phi \, .$$
The leading-order term can actually be calculated analytically: for $z \in \mathbb{C}$ we have
$$ \int \limits_0^\pi \mathrm{e}^{\mathrm{i} \left[\phi - z \mathrm{e}^{-\mathrm{i} \phi}\right]} \, \mathrm{d} \phi = \mathrm{i} \left[(2 \operatorname{Si}(z) - \pi) z + 2 \cos(z)\right] \, . $$
The sine integral $\operatorname{Si}$ satisfies $\lim_{x \to \infty} \operatorname{Si}(x) = \frac{\pi}{2}$, so we obtain
$$ \int \limits_{\gamma_\varepsilon (0)} f_{\alpha,0,\gamma} (z) \, \mathrm{d} z \sim \frac{\mathrm{i} \varepsilon}{\gamma^2} \left[(2 \operatorname{Si}(\alpha^2 \varepsilon^{-1}) - \pi) \alpha^2 \varepsilon^{-1} + 2 \cos(\alpha^2 \varepsilon^{-1})\right] \stackrel{\varepsilon \to 0}{\longrightarrow} 0 $$
as desired. Now we are allowed to apply the residue theorem , which yields
\begin{align}
I(\alpha,0,\gamma) &= \frac{1}{2} \left[ 2 \pi \mathrm{i} \operatorname{Res} (f_{\alpha,0,\gamma}, \mathrm{i} \gamma) - \lim_{R \to \infty} \int \limits_{\Gamma_R} f_{\alpha,0,\gamma} (z) \, \mathrm{d} z - \lim_{\varepsilon \to 0} \int \limits_{\gamma_\varepsilon (0)} f_{\alpha,0,\gamma} (z) \, \mathrm{d} z \right] \\
&= \frac{\pi}{2 \gamma} \exp \left(- \frac{\alpha^2 + \gamma^2}{\gamma}\right) \, .
\end{align}
Almost the same calculation yields
\begin{align}
I(\alpha,\beta,\gamma) &= \frac{1}{2} \left[ 2 \pi \mathrm{i} \operatorname{Res} (f_{\alpha,\beta,\gamma}, \mathrm{i} \gamma) - \lim_{R \to \infty} \int \limits_{\Gamma_R} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z - \lim_{\varepsilon \to 0} \int \limits_{\gamma_\varepsilon (\beta)} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z - \lim_{\varepsilon \to 0} \int \limits_{\gamma_\varepsilon (-\beta)} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z \right] \\
&= \frac{\pi}{2 \gamma} \exp \left(- \gamma \frac{\alpha^2 + \gamma^2}{\beta^2 + \gamma^2}\right)
\end{align}
for $\beta > 0$ .
Combining these results we conclude that the integral is given by
$$ I(\alpha,\beta,\gamma) = \frac{\pi}{2 \gamma} \exp \left(- \gamma \frac{\alpha^2 + \gamma^2}{\beta^2 + \gamma^2}\right) $$
for every $\alpha,\beta \geq 0$ and $\gamma > 0$ .
