By De Moivre's formula $\sin(x)=\frac{e^{ix}-e^{-ix}}{2i}$ we have the following Fourier sine series:
$$\frac{\sin(3x)\sin(4x)\sin(5x)\cos(6x)}{\sin^2(x)}\\= -\frac{1}{2} \sin(2x)-\frac{1}{2}\sin(4x)+\sin(8x)+\frac{3}{2}\sin(10x)+\frac{3}{2}\sin(12x)+\sin(14x)+\frac{1}{2}\sin(16 x)$$
and:
$$I(n)=\int_{0}^{+\infty}\frac{\sin(2nx)}{x\cosh(x)}\,dx = 2\arctan\left(\tanh\frac{\pi n}{2}\right) $$
follows by differentiation under the integral sign. The original integral can so be expressed in terms of the Gudermannian function:
$$ I = \frac{1}{2} \big(-\text{gd}(\pi)- \text{gd}(2\pi) +
2 \text{gd}(4\pi) + 3 \text{gd}(5\pi) +
3 \text{gd}(6\pi) + 2 \text{gd}(7\pi) +
\text{gd}(8\pi)\big) \approx 7.11363 $$
If $\;a<0,\;$ then $\;\sin(e^{ax})<e^{ax},\;$ and
$$I<\int\limits_{17}^\infty e^{ax}\,\text dx,$$
i.e. converges.
If $\;a=0,\;$ then
$$I \ge \dfrac{17\sin1}{17+e} \int\limits_{17}^\infty \dfrac{\ln\ln x}xdx
= \dfrac{17\sin1}{17+e} \int\limits_{17}^\infty \ln\ln x\,\text d\ln x,$$
i.e. diverges.
If $\;a>0,\; x\ge 17 > e^e,\;$ then
$$\;\dfrac{d}{dx}\dfrac{\ln \ln x}{x+e} = \dfrac{e-x(\ln x\cdot\ln\ln x-1)}{x(x+e)^2\ln x} < 0.\tag1$$
Applying substitution $\;x=\dfrac1a\ln y,\; k_0=\left\lceil\dfrac{e^{17a}}{2\pi}\right\rceil,\;$ easily to get
$$I = \dfrac1a\int\limits_{\large e^{17a}}^\infty\; \dfrac{\ln\ln\left(\dfrac1a\ln y\right)}{\dfrac1a\ln y+e}\,\dfrac{\sin y}y\,\text dy,$$
$$aI = \int\limits_{e^{17a}}^{2k_0\pi}\; \dfrac{\ln\ln\left(\dfrac1a\ln y\right)}{\dfrac1a\ln y+e}\,\dfrac{\sin y}y\,\text dy
+ \sum\limits_{k=k_0}^\infty\;\int\limits_{2\pi k}^{2\pi k+2\pi}\; \dfrac{\ln\ln\left(\dfrac1a\ln y\right)}{\dfrac1a\ln y+e}\,\dfrac{\sin y}y\,\text dy.\tag2 $$
From $(1)$ should, that the sinusoids in the integrals under the sum have decreasing weight, and for $\;y\in(2k\pi,2k\pi+\pi),\;k=k_0\dots\infty\;$
$$\dfrac{\ln\ln\left(\dfrac1a\ln y\right)}{\dfrac1a\ln y+e}\,\dfrac{\sin y}y
\ge\left|\dfrac{\ln\ln\left(\dfrac1a\ln(y+\pi)\right)}{\dfrac1a\ln(y+\pi)+e}\,\dfrac{\sin (y+\pi)}{y+\pi}\right|.$$
Then
$$\int\limits_{2\pi k}^{2\pi k+2\pi}\; \dfrac{\ln\ln\left(\dfrac1a\ln y\right)}{\dfrac1a\ln y+e}\,\dfrac{\sin y}y\,\text dy>0.$$
Therefore, should exist the constants $\;C_1,\,C_2,\;$ such as
$$I < C_1 + C_2\int\limits_{2\pi k_0}^\infty\dfrac{\sin y}y\,\text dy,\tag3$$
i.e. converges.
Finally, the given integral converges if $\;a\in\mathbb R\setminus \{0\}.$
Best Answer
First method. Using Poisson summation formula
If a continuous integrable function $\varphi$ and its Fourier transform are rapidly going to zero at infinity (check its Wikipedia page for more details) then $$ \sum_{n=-\infty}^\infty \varphi(n) = \sum_{n=-\infty}^\infty \hat \varphi(n) $$ Since the Fourier transform tends to transform rapidly decreasing (to zero) functions to slowly decreasing (to zero) functions, and vice versa, Poisson's formula is well suited to the calculations of certain series.
You can check that for $\varphi(x) = \pi e^{-2\pi |x|}$ we have $\hat \varphi(s) = \frac{1}{1+s^2}$. Poisson tells us that $$ \sum_{n=-\infty}^\infty \frac{1}{n^2 + 1} = \sum_{n=-\infty}^\infty \pi e^{-2\pi |n|} $$ Use parity and geometric series to complete the computation.
Second method. Using Parseval's identity for Fourier series
Let $e_n(t) = e^{2\pi i nt}$ be the Fourier orthonormal system of $L^2[0,1]$ and recall Parseval's identity : $$ \forall f,g \in L^2[0,1] ~~,~~ \langle f,g\rangle = \int_0^1 f(t)\overline{g(t)} \,\mathrm dt= \sum_{n=-\infty}^\infty \langle f,e_n\rangle \overline{\langle g,e_n\rangle} $$ Take $s \in \mathbf C\setminus \mathbf Z$ and let $z=\overline{s}$. Compute the Fourier coefficients of $f(t) = e^{2 \pi i st}$ and $g(t)=e^{2\pi i zt}$ $$ \langle f,e_n\rangle = \frac{1}{2\pi i} \frac{e^{2\pi i s}-1}{s-n} ~~\text{ and } ~~ \overline{\langle g,e_n\rangle} = -\frac{1}{2\pi i} \frac{e^{-2\pi i s}-1}{s-n} $$ Using Parseval's identity, it comes, after some simplifications : $$ \frac{\pi^2}{\sin^2(\pi s)} = \sum_{n=-\infty}^\infty \frac{1}{(s-n)^2} \tag{1} $$ Break the sum into three pieces, and write $$\frac{\pi^2}{\sin^2(\pi s)} = \frac{1}{s^2} + \sum_{n=1}^\infty \frac{1}{(s-n)^2} + \frac{1}{(s+n)^2} \tag{2} $$ Integrate $(2)$ both sides, then multiply by $-1$ to get the famous formula : $$ \pi \cot(\pi s) = \frac{1}{s} + \sum_{n=1}^\infty \frac{2s}{s^2-n^2} $$ Replacing $s$ by $is$ leads to : $$ \pi \coth (\pi s) = \frac{1}{s} + \sum_{n=1}^\infty \frac{2s}{s^2 +n^2} $$ If you really want to avoid the use of any complex variable, replace $s$ by $is$ in $(1)$ or in $(2)$ then integrate in the « real » sense. Again, it is easy to recover the result given by WolframAlpha.
Also, using Fourier series and Dirichlet's theorem (for pointwise convergence) you may check other topics on MSE, for example. :)
ADDENDUM. We can reach our goal a bit faster. Apply Parseval's identity to $f(t)= e^{2\pi ist}$ and $g(t) = \overline{f(t)}$ : $$ \int_0^1 f(t)^2 \mathrm dt = \sum_{n=-\infty}^\infty \langle f,e_n\rangle \langle f,e_{-n}\rangle $$ Substitute things... $$ \frac{1}{2\pi i} \frac{e^{4 \pi i s}-1}{2s} = \sum_{n=-\infty}^\infty \bigg(\frac{e^{2\pi is} - 1}{2\pi i} \bigg)^2 \frac{1}{s^2-n^2} $$ Simplify some terms... $$ 2\pi \cot(\pi s) = \sum_{n=-\infty}^\infty \frac{2s}{s^2-n^2} $$ Finally, by symmetry in $n\neq 0$ $$ \pi \cot(\pi s) = \frac{1}{s} + \sum_{n=1}^\infty \frac{2s}{s^2-n^2} $$