Numerically calculating inverse Laplace via the inverse Laplace transformation formula

I'm trying to simulate a control system whose transfer function is $H(s)$. I'm comparing different numerical methods for this. I have already used these two methods:
- Converting the transfer function to the state space form, then iterating then system.
- Converting the transfer function into the discrete form by bilinear transformation, then solving the difference equation.

Next, I want to do the simulation by directly solving the inverse Laplace transformation formula:

$\mathcal{L}^{-1} \{H(s)\} = f(t) = \frac{1}{2\pi i}\lim_{T\to\infty}\int_{\sigma-iT}^{\sigma+iT}e^{st}H(s)\,ds, \quad \mbox{(where } \quad \sigma = Re\{s\} \mbox{)}$

The problem is that I don't know about complex integration. A typical transfer function I'm dealing with is $H(s) = \frac{1}{s+4}$. For this transfer function, what values of $\sigma$ and $T$ should I choose? Assuming that the system will always be real, should the integration return a real number (absolute value of the result), or a complex number?

I'm going to use C++ language for the implementation, and use rectangular rule to solve the integral. What must be the maximum step size of the integration for not seeing any error by eye when I look at the plotting of the output signal?

-

$\sigma$ has to be chosen such that it is greater than the real part of all singularities of $H$, in your case greater than $-4$. $T$ should just be large enough (your $H$ has an obvious decay for $\Im(s) \to \infty$, you could try to work out estimates on the error for the truncation of the integral).