Thanks, GH! Let me have another go. I think the following is the right way to go about things, at least if one wants something self-contained and with good, explicit constants. (The latter more or less implies the former, given that almost all of the literature is non-explicit.)
We want to estimate
$$\frac{1}{2\pi i} \int_{\sigma-i\infty}^{\sigma+i \infty} \left|\frac{1}{s} -G(s)\right|^2 |\zeta(s)|^2 ds,$$
where $G(s)$ is the Mellin transform of a well-chosen function $g:\lbrack 0,\infty)\to \mathbb{R}$. It is easy to see that $G(s) \zeta(s)$ is the Mellin transform of $x\mapsto \sum_n g(n x)$. We will choose $g$ so that (a) $G(\sigma+it)$ is small for $|t|\geq T$, (b) the "physical-space" estimation we are about to do is easy.
By Plancherel,
$$\frac{1}{2\pi i} \int_{\sigma-i\infty}^{\sigma+i \infty} |1-G(s)|^2 \frac{|\zeta(s)|^2}{|s|^2} ds = \int_0^\infty |h(x)|^2 x^{2\sigma-1} dx,$$
where $$h(x) = \lfloor 1/x\rfloor - \sum_n g(n x).$$ (We will make sure that $G(1)=1$, so that there is no pole at $s=1$; in this way, the equation above will hold for $\Re(s)>0$, and not just for $\Re(s)>1$.)
First, let us show that $h(x)$ is bounded for all $x$. (This part is the same as what I had before.) By second-order Euler-Maclaurin,
$$\sum_n g(n x) = \frac{1}{x}
\int_0^\infty g(t) dt - \frac{g(0)}{2} - \frac{g'(0)}{12} + \textrm{err},$$
where $|\textrm{err}| \leq \frac{x}{24} \int_0^\infty |g''(t)| dt$. We will work with $g$ such that $g(0)=1$, $g'(0)=0$ and $\int_0^\infty g(t) dt = 1$. Then $$|h(x)| = \left|\lfloor 1/x\rfloor - \left(\frac{1}{x} - \frac{1}{2} + \text{err}\right)\right|\leq \frac{1}{2} + \textrm{err}$$ for all $x$.
We will now choose $g$ so that we can give a much better estimate for $x$ not too small. For starters, we will have $g(x)=1$ for $x\leq 1-\delta$ and $g(x)=0$ for $x\geq 1+\delta$. Then $$h(x) = 0$$ unless $n x \in \lbrack 1-\delta,1+\delta\rbrack$. Moreover, for each $x>2\delta$, there is at most one $n$ such that $n x$ is in that interval, since $\frac{1+\delta}{x} - \frac{1-\delta}{x} = \frac{2 \delta}{x} < 1$. Hence $$\begin{aligned}\int_0^\infty |h(x)|^2 x^{2\sigma-1} dx &\leq \int_0^{2\delta} c^2 x^{2\sigma-1} dx + \sum_{n\leq \frac{1}{2 \delta}} \int_{\frac{1-\delta}{n}}^{\frac{1}{n}} |1-g(n x)|^2 x^{2\sigma-1} dx \\ &+ \sum_{n\leq \frac{1+\delta}{2\delta}} \int_{\frac{1}{n}}^{\frac{1+\delta}{n}} |g(n x)|^2 x^{2\sigma-1} dx,\end{aligned}$$ where $c =1/2 + \frac{2\delta}{24} \int_0^\infty |g''(t)| dt$. Obviously $$\int_0^{2\delta} c^2 x^{2\sigma-1} dx = c^2 \frac{(2\delta)^{2\sigma}}{2\sigma}.$$ To estimate the two other integrals, we have to choose a convenient $g$ obeying our conditions. I will simply take $$g(x) = \begin{cases} 1 &\text{if $x\leq 1 - \delta$,}\\ \frac{1+\delta-x}{2\delta} &\text{if $1-\delta <x<1+\delta$,}\\0 &\text{if $x\geq 1 + \delta$.}\end{cases}$$ Hence $$\begin{aligned} \int_{\frac{1-\delta}{n}}^{\frac{1}{n}} |1-g(n x)|^2 x^{2\sigma-1} dx &+ \int_{\frac{1}{n}}^{\frac{1+\delta}{n}} |g(n x)|^2 x^{2\sigma-1} dx\\ &\leq \frac{2 \eta}{n^{2\sigma}} \int_1^{1+\delta} \left(\frac{1+\delta-t}{2\delta}\right)^2 dt = \frac{\delta \eta}{6 n^{2\sigma}},\end{aligned}$$ where $\eta=1$ if $1/2\leq \sigma\leq 1$ and $\eta = \frac{(1-\delta)^{2\sigma-1} + (1+\delta)^{2\sigma-1}}{2}$ if $0<\sigma\leq 1/2$. Now, for any $y\geq 0$, $$\sum_{n\leq y} \frac{1}{n^{2\sigma}} \leq \begin{cases} \zeta(2 \sigma)& \text{if $\sigma>1/2$,}\\ \log(2 y + 1) &\text{if $\sigma=1/2$,}\\ \frac{y^{1-2\sigma}}{1-2\sigma} & \text{if $\sigma<1/2$.}\end{cases}$$
Taking totals, we conclude that $$\int_0^\infty |h(x)|^2 x^{2\sigma-1} dx \leq c^2 \frac{(2\delta)^{2\sigma}}{2\sigma} + \begin{cases} \frac{\zeta(2\sigma)}{6} \delta& \text{if $\sigma>1/2$,}\\ \frac{\delta}{6} \cdot \log\left(\frac{1}{\delta} + 2\right) &\text{if $\sigma=1/2$,}\\ \frac{\eta (1+\delta)^{1-2\sigma}}{6\cdot 2^{1-2\sigma} (1-2\sigma)} \cdot \delta^{2\sigma} & \text{if $\sigma<1/2$.}\end{cases}$$
Now it just remains to estimate how much of the tail we captured. A quick calculation shows that we then have $$G(s) = \frac{(1+\delta)^{s+1} - (1-\delta)^{s+1}}{2\delta (s+1) s}.$$ Not unexpectedly, this is close to $1/s$ for $\Im(s)$ small. Even more to the point, for any $s$, $$|G(s)| \leq \frac{(1+\delta)+(1-\delta)}{2\delta |s+1| |s|} = \frac{1}{\delta |s+1| |s|},$$ and so $$|1-G(s) s|^2 \geq \left(1-\frac{1}{\delta |s+1|}\right)^2 \geq \left(1 - \frac{1}{\delta T}\right)^2$$ for every $s=\sigma+it$ with $|t|\geq T$. Hence $$\int_{\sigma-i\infty}^{\sigma-iT} + \int_{\sigma+iT}^{\sigma+i\infty} \frac{|\zeta(s)|^2}{|s|^2} ds \leq \frac{1}{\left(1 - \frac{1}{\delta T}\right)^2} \int_0^\infty |h(x)|^2 x^{2\sigma-1} dx .$$
Minimizing what will be the main term, we choose $\delta = 3/T$ if $\sigma\geq 1/2$, and $\delta = \frac{1+1/\sigma}{T}$ if $0<\sigma<1/2$. Then $\frac{1}{\left(1 - \frac{1}{\delta T}\right)^2}$ equals $9/4$ for $\sigma\geq 1/2$ and $(1+\sigma)^2$ for $0<\sigma<1/2$.
We conclude that $$\frac{1}{2\pi i} \int_{\sigma-i\infty}^{\sigma-iT} + \int_{\sigma+iT}^{\sigma+i\infty} \frac{|\zeta(s)|^2}{|s|^2} ds \leq \begin{cases} \frac{9 \zeta(2\sigma)}{8} \cdot \frac{1}{T} + \frac{9 c^2 6^{2\sigma}}{8\sigma} \cdot \frac{1}{T^{2\sigma}} & \text{if $\sigma>1/2$,}\\ \frac{9}{8} \frac{\log\left(\frac{T}{3}+2\right)+12 c^2}{T} &\text{if $\sigma=1/2$,}\\ \left(\frac{ \frac{1}{2} \left(\frac{1+\delta}{1-\delta}\right)^{1-2\sigma} + \frac{1}{2}}{6\cdot 2^{1-2\sigma} (1-2\sigma)} + c^2 \frac{2^{2\sigma}}{2\sigma}\right)\cdot \frac{(1+\sigma)^{2\sigma+2}}{\sigma^{2\sigma}}\cdot \frac{1}{T^{2\sigma}} & \text{if $\sigma<1/2$.}\end{cases}$$
Oh, by the way, $\int_0^\infty |g''(x)| = 1/\delta$, so $c = 7/12$.
No doubt the constants can still be improved, and any suggestions on how to simplify things further are very welcome. I think I can be happy now, though.