There is a good reason to use $z$ instead of $e^{sT}$. Before starting with any analysis, let me remind you that in analysis of signals and systems we are interested in analyzing the frequency spectrum of the signal, i.e. the Laplace transform on the imaginary line $s = jw$. And since your signal $x[n]$ is discrete, then its frequency spectrum its periodic, so its more general define $s=jwT$.
Now, let $X(z) = \mathcal{Z}\{x[n]\}$ of a causal or non-causal discrete signal $x[n]$, i.e.
$$ X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}. $$
Since $z\in\mathbb{C}$ we have $z = |z| e^{j\arg z}$. Without loss of generality we rewrite $|z| = r$ and $\arg z = wT$, i.e. $z=r e^{jwT}$ (note that not necessarily $r=1$). Then
$$ \begin{aligned}
X(z) &= \sum_{n=-\infty}^{\infty} x[n] z^{-n}\\
&= \sum_{n=-\infty}^{\infty} x[n] (r e^{jwT})^{-n}\\
% &= \sum_{n=-\infty}^{\infty} (x[n] r^{-n}) e^{-njwT}\\
&= \sum_{n=-\infty}^{\infty} (x[n] r^{-n}) (e^{jwT})^{-n},
\end{aligned} $$
which implies $X(z) = \left. \mathcal{L}\{x[n]r^{-n}\} \right\rvert_{s=jwT} = \mathcal{F}\{x[n]r^{-n}\}$. As a consequence, $X(z)$ is a Fourier transform more generic than the Fourier transform $X(e^{jwT}) = \mathcal{F}\{x[n]\}$ of our signal of interest.
So, if the convergence radius of $X(z)$ is less than unity then $X(e^{jwT})$ does not exist and therefore its Fourier transform does not either, which represents a problem because there are many signals with this problem of convergence, e.g. non-causal signals such as a digital image filter. Therefore, it is convenient (and even necessary in non-causal signals) to use the Z-transform.
Or informally, use $z$ instead of $\left. e^{sT} \right\rvert_{s=jwT} = e^{jwT}$ whenever you can.
We also recommend to see this link about radius convergence.
At this point, it is clear that the Z-transform has the same objective as the Laplace transform: ensure the convergence of the transform in some region of $\mathbb{C}$, where the Z-transform does it for discrete signals and Laplace transform for continuous signals.
Related techniques (I), (II). Using the fact about the Laplace transform $L$ that
$$ L(f*g)=L(f)L(g)=F(s)G(s)\implies (f*g)(t)=L^{-1}(F(s)G(s)) .$$
In our case, given $ H(s)=\frac{1}{(s+1)}\frac{s}{(s^2+4)}$
$$F(s)=\frac{1}{s+1}\implies f(t)=e^{-t},\quad G(s)=\frac{s}{s^2+4}\implies g(t)=\cos(2t).$$
Now, you use the convolution as
$$ h(t) = \int_{0}^{t} e^{-(t-\tau)}\cos(2\tau) d\tau . $$
Best Answer
The LT of a convolution of two functions is the product of the individual LTs. If you recall, we found that
$$\hat{f_T}(p) = \frac{2 - e^{-p T}}{p}$$
Thus, the LT of the above convolution is simply
$$\hat{h_{S*T}}(p) = \left ( \frac{2 - e^{-p S}}{p} \right )\left ( \frac{2 - e^{-p T}}{p} \right )$$