Let $\sigma$ be represented by a PD matrix $A$ and $\rho$ by $A+B$. Note that $|B|\ge B$ (in the sense of PSD matrices) and has the same $1$-norm. Also $\Pi |B|\Pi\ge \Pi B \Pi$, so to dominate $\Pi B\Pi $ is always easier than to dominate $\Pi|B|\Pi$. Thus, we can replace $B$ by $|B|$ and assume that $B$ is PD as well, so its trace norm is just its trace.

Now let's take $a<1$ and consider the spectral decomposition of $H$ in which we unite all the eigenvalues of $A$ between $a^{k+1}$ and $a^k$, i.e., we write $H=\oplus H_k$ such that $H_k$ is an invariant subspace of $A$ and $a^{k+1}I\le A\le a^k I$ on $H_k$. Let $d_k$ be the dimension of $H_k$ and let $P_k$ be the projector to $H_k$.

Now we will construct our projector $\Pi$ as the sum of $\Pi_k P_k$ where $\Pi_k$ is a projection within $H_k$, i.e., if $H\ni x=\sum_k x_k$ is the orthogonal decomposition of $x$ with $x_k\in H_k$, we'll have $\Pi x=\sum_k \Pi_k x_k$. Notice that then $\Pi A\Pi\le a^{-1}A$ regardless of the choice of $\Pi_k$ and, if you think a bit, you will realize that this is essentially almost all you can do to ensure that inequality.

Now take a single $H_k$ and consider the operator $P_kBP_k$. Let its trace be $\mu_k$ so that $\sum_k\mu_k=\varepsilon$. Take $\lambda>0$ and remove from $H_k$ the eigenvectors of this operator with eigenvalues greater than $\lambda a^{k+1}$. The dimension of the removed space will be $\le \frac{\mu_k}{\lambda a^{k+1}}$.

On the remaining subspace $V_k\subset H_k$, we have $\langle Bx_k,x_k\rangle\le\lambda \langle Ax_k,x_k\rangle$. Let now $x=\sum_k x_k$, $x_k\in V_k$. Then, choosing some integer $K$, we can write
$$
\langle Bx,x\rangle\le \sum_{k,m} |\langle Bx_k,x_m\rangle|=
\\
\sum_k \langle Bx_k,x_k\rangle
+2\sum_{k,m:k-K\le m\le k-1} |\langle Bx_k,x_m\rangle|
\\
+2\sum_{k,m:m< k-K} |\langle Bx_k,x_m\rangle|
\\
=\Sigma_1+\Sigma_2+\Sigma_3\,.
$$
We have $\Sigma_1\le\lambda\langle Ax,x\rangle$. Also
by Cauchy-Schwarz and the positive definiteness of $B$, we have $\Sigma_2\le 2K\Sigma_1$. Thus, the sum of the first two terms is at most $(2K+1)\lambda\langle Ax,x\rangle$. The question is, of course, what to do with $\Sigma_3$. And the answer is that we will just kill it entirely by further reducing $x_k$ to the intersection of the kernels of the corresponding $P_mB$ on $H_k$. Since $P_mB$ is an operator of rank $d_m$ at most, that will reduce the dimension of $V_k$ by at most $\sum_{m:m<k-K}d_m$.

Then we shall have $\Pi B\Pi\le (2K+1)\lambda \Pi A\Pi$ and
$$
\Pi(A+B)\Pi\le [1+(2K+1)\lambda]\Pi A\Pi\le [1+(2K+1)\lambda]a^{-1}A\,.
$$
On the other hand, the codimension of the final $V_k$ in $H_k$ is at most $D_k=\frac{\mu_k}{\lambda a^{k+1}}+\sum_{m:m<k-K}d_m$, so the trace lost is at most
$$
\sum_k a^k D_k=\sum_k \frac{\mu_k}{\lambda a}+\sum_{k,m:m<k-K}a^k d_m
\\
=\frac{\varepsilon}{\lambda a}+\sum_m d_m a^{m+1}\frac{a^K}{1-a}
\\
\le \frac \varepsilon{\lambda a}+\frac{a^K}{1-a}
$$
because $\sum_m d_m a^{m+1}\le \operatorname{tr} A=1$.

Now one can just play with the parameters to balance.
Suppose we want the domination with the constant $1+C\delta$ and the trace loss at most $C\Delta$ with $C$ being some absolute constant (say, $5$). Then we are forced to take $a=1-\delta$ and $(2K+1)\lambda\le 3\delta$. We also need $\frac{a^K}{1-a}<\Delta$, which calls for $K=\delta^{-1}\log\frac{1}{\delta\Delta}$ and $\lambda=\delta^2 \left(\log\frac{1}{\delta\Delta}\right)^{-1}$, so our $\varepsilon$ should be less that $\Delta\delta^2\left(\log\frac{1}{\delta\Delta}\right)^{-1}$ to make the conditions compatible. If $\delta=\Delta$, then we get them both around $\sqrt[3]{\varepsilon\log\frac 1\varepsilon}$. This is, probably, not optimal, but you requested just *some* speed of tending to $0$ with $\varepsilon$ and the power bound is not terribly bad, so I'll stop here.