Hi dp,
Using the approach I mentioned above, that is splitting the $\csc(\pi z)$ term into its pole part $\frac{1}{\pi z}$ and the rest, I get an exact answer for the pole part and a series for the rest. The pole term yields a principal value integral that can be computed analytically, plus the contribution from integrating around an infinitesimal semi-circle, giving
$$ I_1 = i\pi + i\pi\ \mbox{erf} \left( \frac{a}{2 b} \right), $$
where $\mbox{erf}()$ is the error function.
The remaining contribution can be evaluated as an asymptotic
series for large $b$. We write
$$ \csc(\pi u ) - \frac{1}{\pi u} = \sum_{m=1} c_{2 m-1} ( u)^{2m -1} ,$$
which has a radius of convergence of $\vert u\vert \le 1.$ Therefore, if
we substitute this in the integrand, the resulting series will be an
asymptotic series for large $b$,
since the series has a finite radius of convergence. Since this integral has no
singularity on the real axis (please note, I have changed variables using
$z = i u$, so the integration is along the real axis),
we may take $\epsilon=0$ and write
$$
I_2 \approx i\pi \int_{-\infty}^{\infty} \exp(i a u)
\exp(-b^2 u^2) \sum_{m=1}
c_{2m-1} i^{2m-1} u^{2m-1} .
$$
For clarity write $b^2 \rightarrow g.$ We use
$u^{2m} = \left( -\frac{\partial^m}{\partial g^m}\right)
\exp(-g u^2).$ Reversing the order of the sum and integral and performing
the integral gives
$$
I_2 \approx i \pi \sum_{m=1} c_{2m-1}
\frac{\partial^m}{\partial g^m}
\pi\ \mbox{erf}\left( \frac{a}{2 \sqrt{g}} \right) ,
$$
where $g = b^2$.
May I ask, from what problem does this integral arise? Is it from an inverse
Mellin transform used to perform some other integral? Do you expect the parameter $b$ to be significantly larger than one?
I add this - numerically integrating your original integral
in Mathematica for $a=2, b=2$ gives $I \approx 4.56 i$ while
the main term above gives $I_1 \approx 4.78 i $. For $a=2, b=4$ the
results are $3.97 i, 4.01 i.$
I will include additional details here later, with some results for the second term $I_2.$
EDIT - For small $b$, one can evaluate the integral by expanding the term $e^{b^2 z^2}$ in a series.
The term coming from the principal value integral is kept the same as shown above, while the other integral
can be written as
$$ I_2 = \pi \int_{-i\infty}^{i\infty} dz e^{a z + b^2 z^2} \left( \csc(\pi z) -\frac{1}{\pi z} \right), $$
where we have set $\epsilon=0$ since the singularity is removed.
Putting $z = i u$ leads to
$$ I_2 = 2 i \pi \int_0^{\infty} du\ e^{-b^2 u^2} \sin(a u) \left(\mbox{csch}(\pi u) - \frac{1}{\pi u}\right) .$$
Now expand the exponential, and use the relationship
$$ \frac{\partial^{2k}}{\partial a^{2k}} \sin(a u) = (-1)^k u^{2k} sin(a u) $$
to obtain
$$ I_2 = 2 i \pi \sum_{k=0}^{\infty} \frac{b^{2k}}{k!} \frac{\partial^{2k}}{\partial a^{2k}} \int_0^{\infty} du\
\sin(a u) \left( \mbox{csch} (\pi u) - \frac{1}{\pi u} \right). $$
The integral can be evaluated to give $-1/(e^a + 1)$. Replacing it in the series leads to a useful series for small $b$,
so that
$$ I = i\pi + i\pi\ \mbox{erf}\left( \frac{a}{2 b} \right) - 2 i \pi \sum_{k=0}^{\infty} \frac{b^{2 k}}{k!} \frac{\partial^{2k}}{\partial a^{2k}} \frac{1}{e^a+1} .$$
I (or you) need to analyze the radius of convergence of the series; numerical tests I have done show it
to converge quickly for small $b/a$ much less than 1, while for larger values it does not seem to converge. Will get back to you on that.
May I ask, what are you going to use this for?
EDIT 2: The best I have been able to come up with is the following (for small $b$): your integral is equal to, in the limit of zero $\epsilon$,
$$ I = i\pi + 2 i \pi \int_0^{\infty} du\ \frac{\sin(a u)}{\sinh (\pi u)} e^{-b^2 u^2} .$$
The integral can be evaluated for $b=0$:
$$ \int_0^{\infty} du\ \frac{\sin(a u)}{\sinh (\pi u)} = \frac{1}{2} \tanh\left( \frac{a}{2} \right) $$.
Now, if $b$ is small, we can expand the exponential and obtain the required integrals by differentiating wrt to $a$; these leads to a series that converges rapidly if $a$ is large enough. That is, the series is useful for small $b/a$. The result is
$$ I \approx i\pi + i\pi \sum_{k=0} \frac{b^{2k}}{k!} \frac{\partial^{2k}}{\partial a^{2k}} \tanh\left( \frac{a}{2} \right) $$
For example, for $a=2$, if we compute the difference between the numerically integrated integral and the sum, taken to 11 terms, for the values $b=1/m, m=1,2,\cdots, 10$, we find
14.568696715793862860,
2.8724585353959*10^-6,
2.410511554*10^-10,
2.718473*10^-13,
1.3485*10^-15,
2.55*10^-17,
2.53*10^-17,
4.93*10^-17,
7.79*10^-17,
1.077*10^-16
So, for $a=2, b=1$ the series is only asymptotic; you get a much better answer using just a couple of terms, but for small $b/a$ it works quite well. I suggest doing some numerical experiments yourself to see if this answer meets your needs.
Hope this helps,
Tom