In response to the questioner's comment: one standard viewpoint is that distributions (or tempered distributions, similarly) are "weak-dual-topology" limits of test functions (or of Schwartz functions), where the weak-dual topology (also called "weak-star topology) is given by semi-norms $\mu_\varphi(u)=|u(\varphi)|$ for test functions (or Schwartz) $\varphi$, at first for other test functions (or Schwartz $\varphi$) and then on the (quasi-) completion, by continuity. Thus, in this viewpoint, the integral representation of $\delta$ is not a literal integral, but is an extension-by-continuity of literal integrals giving approximations to $\delta$. I think this viewpoint, of having integrals not being literal, but extensions-by-continuity, may be most consonant with a common physics-y viewpoint... and does (to my mind) convey a more natural sense of "generalized functions/distributions", not as "functionals", but as an extended notion of function.
Of course, the second viewpoint, perhaps more standard in mathematics per se, is that distributions are continuous linear functionals, and that the integrals written about don't make sense... again, at least not as literal integrals, but perhaps more... haughtily said? That is, the same expressions are obtained, but if/because one is insisting on literal integrals, implicitly disallowing writing integral signs for extensions-by-continuity, some new symbols must be made up. (As you can tell, I am not deeply moved by this viewpoint.)
Either way, the underlying mathematics and mechanisms are the same. That is, e.g., for each Schwartz function $\varphi$, the integrals against the truncated Fourier inversion expression for $\delta$ go to $\varphi(0)$ as the cut-off goes to infinity.
The non-uniformity issue/feature is that there is (provably) no estimate for this limit-approach uniform in $\varphi$. It depends, as it turns out, on $\varphi$.
Is this to the point?
EDIT: I forgot to make explicit the point that, whether or not the integrals are pointwise/numerically valued, literal-or-not, etc., the operator/functional indicated by the integral (as extended by continuity, say) has certain expected properties, so that the usual manipulations which would be correct for numerical/literal integrals are also correct for the non-literal (extended-by-continuity) maps written as integrals... because they extend those integrals, and have the same properties... apart, perhaps, from non-critical numerical/pointwise evaluability.
This is why many physics-y heuristics succeed, after all: what is written is not what is literally meant (from the most exclusionary viewpoint), and it is "still ok" because things extend by continuity, with the same correct transformation properties of whatever-it-is operator denoted by "integral" (and extended from a literal integral).
Yes, on good days, we have our cake and eat it, too. This is (to my mind) the marvelous-ness of mathematics, both internally and in applications.
True, on a bad day bad things happen which will be inexplicable to the (naive?) optimist. Yet, crazily enough, even those bad things can often be repurposed to interesting ends. That's the second miracle in this business, I think.
Edit-edit: following up on the edit to the original question, yes, in that example, for $f$ represented by its Taylor-Maclaurin series in a fixed interval $[-\varepsilon_o,\varepsilon_o]$, as $\varepsilon\to 0$ we can estimate the rate of approach to $f(0)$ uniformly.
But if we look at the aggregate of the classes of such functions over all $\varepsilon_o>0$, this uniform-ness is lost.
Also, although I'd not advocate making this one's first worry, there is an old theorem of E. Borel that says that an arbitrary sequence of real (or complex) numbers can be the Taylor-Maclaurin coefficients at $0$ of a smooth function (obviously, then, not necessarily represented by the infinite series, but the finite-series-with-error-term still works...). This shows, strikingly, that there's no way to hope for true uniformity.
Also, as a sort of paraphrase of the Borel resuls: in the classification of distributions supported at $0$ (essentially proven by Taylor-Maclaurin expansions), we find that they are only the finite linear combinations of $\delta$ and its derivatives. Thus, there is no sequence $c_n$ such that $\sum_{n\ge 1} c_n\cdot \delta^{(n)}$ converges as a distribution.
If we use 3D rotational symmetry and define
$$x^{\mu}~=~(T,0,0,\underbrace{R}_{\geq 0})\quad\text{and}\quad x^{\prime \mu}~=~(t,r\sin\theta\cos\phi,r \sin\theta \sin\phi,r\cos\theta),$$
and
$$\rho~=~\sqrt{r^2+\alpha^2}\quad\text{and}\quad \Delta~=~R^2-T^2, $$
then OP's integral (1) becomes
$$ \begin{align}
I(R,T,\alpha)~=~& \int_{\mathbb{R}} \!\mathrm{d}t \int_{\mathbb{R}_+} \!r^2 \mathrm{d}r \int_0^{\pi} \!\sin\theta \mathrm{d}\theta\int_0^{2\pi} \!\mathrm{d}\phi~\cr
&\delta(\rho^2+R^2-2Rr\cos\theta-(T-t)^2)~\Theta(T-t) ~\delta(\rho^2-t^2)~\Theta(t)\cr
~\stackrel{R\neq 0}{=}& 2\pi \Theta(T) \int_0^T \!\mathrm{d}t \int_{\mathbb{R}_+} \!r^2 \mathrm{d}r \int_{-1}^1 \! \mathrm{d}\cos\theta~\cr
&\frac{1}{2Rr}\delta\left(\frac{\rho^2+R^2-(T-t)^2}{2Rr}-\cos\theta\right) ~\frac{\delta(\rho-t)}{\rho+t}\cr
~=~& 2\pi \Theta(T) \int_{\mathbb{R}_+} \!\frac{r^2 \mathrm{d}r}{2Rr~2\rho} ~\Theta\left(\left|\frac{\rho^2+R^2-(T-\rho)^2}{2Rr}\right|<1\right) ~\Theta(0<\rho<T)\cr
~=~& 2\pi \Theta(T) \int_{\alpha}^{\infty} \!\frac{ \mathrm{d}\rho}{4R} ~\Theta(|\Delta +2T\rho|< 2Rr) ~\Theta(\rho<T)\cr
~=~& \frac{\pi}{2R} \Theta(T-\alpha) \int_{\alpha}^T \! \mathrm{d}\rho ~\Theta(4R^2(\rho^2-\alpha^2)-(\Delta+2T\rho)^2)\cr
~=~& \frac{\pi}{2R} \Theta(T-\alpha) \int_{\alpha}^T \! \mathrm{d}\rho ~\Theta(A\rho^2+B\rho+C), \end{align}$$
where
$$A~=~4\Delta, \quad B~=~-4T\Delta \quad\text{and}\quad C~=~-\Delta^2 -4R^2\alpha^2.$$
The discriminant is
$$ D~=~B^2-4AC~=~16R^2\Delta(\Delta +\alpha^2).$$
The roots are
$$ \rho_{\pm}~=~\frac{-B\pm\sqrt{D}}{2A}~=~\frac{T\pm R\sqrt{1+\alpha^2/\Delta}}{2}. $$
Case $-\alpha^2\leq \Delta\leq 0$: Then $D\leq 0$ and the integral $I(R,T,\alpha)=0$ vanishes.
Case $\Delta> 0$: Then $D>0$ and the roots $\rho_{\pm}$ are outside the integration region $[\alpha,T]$, so the integral $I(R,T,\alpha)=0$ vanishes.
Case $\Delta \leq -\alpha^2$: Then $R\leq T$, and therefore $R\sqrt{1+\alpha^2/\Delta}\leq T$, so the root $\rho_+< T$.
$$I(R,T,\alpha)~=~\frac{\pi}{2R} \Theta(T-\alpha)~\Theta(-\Delta-\alpha^2)~\underbrace{m([\rho_-,\rho_+]\cap [\alpha,T])}_{= \text{ length of the interval}}
,$$
where $m$ denotes the Lebesgue measure.
Let us for the remainder of this answer assume that $\alpha=0$. Then the discriminant $D\geq 0$. The roots are
$$ \rho_{\pm}~=~\frac{-B\pm\sqrt{D}}{2A}~=~\frac{T\pm R}{2}. $$
It follow that the integral is supported in the future light-cone:
$$I(R,T,\alpha=0)~=~\frac{\pi}{2R} \Theta(T)~R\Theta(-\Delta)~=~\frac{\pi}{2}\Theta(T)\Theta(-\Delta).$$
(The case $R=0$ follows e.g. from continuity of the limit $R\to 0^+$.)
Best Answer
Let us for concreteness consider a 3D position space $\mathbb{R}^3$. (The generalization to other dimensions is straightforward.) Consider position operators $\hat{\bf r}_1$ and $\hat{\bf r}_2$ for 2 non-identical particles.
In a rigged Hilbert space we define eigenkets $$|{\bf r}_1,{\bf r}_2\rangle, \tag{1} $$ such that $$ \hat{\bf r}_1 |{\bf r}_1,{\bf r}_2\rangle~=~{\bf r}_1|{\bf r}_1,{\bf r}_2\rangle, \qquad \hat{\bf r}_2 |{\bf r}_1,{\bf r}_2\rangle~=~{\bf r}_2|{\bf r}_1,{\bf r}_2\rangle, \tag{2}$$ and $$ \langle {\bf r}_1,{\bf r}_2 | {\bf r}^{\prime}_1,{\bf r}^{\prime}_2\rangle ~=~\delta^3({\bf r}_1\!-\!{\bf r}^{\prime}_1)\delta^3({\bf r}_2\!-\!{\bf r}^{\prime}_2) .\tag{3}$$
Define operators $$ \hat{\bf r}_{\rm cm}~:=~\frac{\hat{\bf r}_1+\hat{\bf r}_2}{2},\qquad \hat{\bf r}_{\rm rel}~:=~\hat{\bf r}_1-\hat{\bf r}_2, \tag{4} $$ and corresponding kets $$|{\bf r}_{\rm cm},{\bf r}_{\rm rel}\rangle\!\rangle~:=~|{\bf r}_1,{\bf r}_2\rangle, \qquad {\bf r}_{\rm cm}~\equiv~\frac{{\bf r}_1+{\bf r}_2}{2},\qquad {\bf r}_{\rm rel}~\equiv~{\bf r}_1-{\bf r}_2. \tag{5} $$
One may use eqs. (2), (3), (4) & (5) to show that the kets (5) are eigenkets $$ \hat{\bf r}_{\rm cm} |{\bf r}_{\rm cm},{\bf r}_{\rm rel}\rangle\!\rangle~=~{\bf r}_{\rm cm}|{\bf r}_{\rm cm},{\bf r}_{\rm rel}\rangle\!\rangle, \qquad \hat{\bf r}_{\rm rel} |{\bf r}_{\rm cm},{\bf r}_{\rm rel}\rangle\!\rangle~=~{\bf r}_{\rm rel}|{\bf r}_{\rm cm},{\bf r}_{\rm rel}\rangle\!\rangle, \tag{6}$$ and normalized $$ \langle\!\langle {\bf r}_{\rm cm},{\bf r}_{\rm rel} | {\bf r}^{\prime}_{\rm cm},{\bf r}^{\prime}_{\rm rel}\rangle\!\rangle ~=~\delta^3({\bf r}_{\rm cm}\!-\!{\bf r}^{\prime}_{\rm cm})\delta^3({\bf r}_{\rm rel}\!-\!{\bf r}^{\prime}_{\rm rel}). \tag{7}$$ In eq. (7) we used how higher-dimensional Dirac delta distributions transform under coordinate transformations, cf. e.g. my Math.SE answer here. It is important that the Jacobian of the transformation (4) is one.
Similarly, one may show that the overlap is $$ \langle {\bf r}_1,{\bf r}_2 | {\bf r}^{\prime}_{\rm cm},{\bf r}^{\prime}_{\rm rel}\rangle\!\rangle ~=~\delta^3(\frac{{\bf r}_1+{\bf r}_2}{2}\!-\!{\bf r}^{\prime}_{\rm cm})\delta^3({\bf r}_1-{\bf r}_2\!-\!{\bf r}^{\prime}_{\rm rel}), \tag{8}$$ etc.