I don't know how to avoid the use of series, but this would be something with them:
Split the integral into two integrals, one over ]-\infty,\beta], and one over [\beta,\infty[. Then substitute the geometric series
<br />
\frac{1}{e^{x-\beta} + e^{-(x-\beta)}} = \frac{1}{e^{-(x-\beta)}}\sum_{n=0}^{\infty} (-1)^n e^{2n(x-\beta)},\quad\quad\quad x<\beta<br />
and
<br />
\frac{1}{e^{x-\beta} + e^{-(x-\beta)}} = \frac{1}{e^{x-\beta}}\sum_{n=0}^{\infty} (-1)^n e^{-2n(x-\beta)},\quad\quad\quad x>\beta.<br />
If I looked this right, now you should get such series for the integrand, that you know how to integrate each term in the series. Of course there's lot of work to be carried out, and in the end the result is in a form of series, so this is not the most desirable way to get the result... I'll be waiting eagerly to see if somebody has better ideas.
edit: oh no. I made one mistake. It is simple to integrate symmetric Gaussian peak over a domain ]-\infty, 0] or [0,\infty[, but actually splitting the integral of the original problem at \beta results is something more difficult. I don't think my idea is working. But it could work if \beta=0.
edit edit: On the other hand, if \beta=0, then it is clear that the integral is zero, so actually I didn't help in anything
