This integral is equivalent to the special case $\alpha=\gamma=p$, $\phantom.\beta=\delta=1/2$ of the known contour integral
$$
\int_{-i\infty}^{i\infty} \Gamma(\alpha+s) \phantom. \Gamma(\beta+s) \phantom. \Gamma(\gamma-s) \phantom. \Gamma(\delta-s) \phantom. ds
= 2\pi i \frac{\Gamma(\alpha+\gamma) \phantom. \Gamma(\alpha+\delta) \phantom. \Gamma(\beta+\gamma) \phantom. \Gamma(\beta+\delta)} {\Gamma(\alpha+\beta+\gamma+\delta)},
$$
valid for all $\alpha,\beta,\gamma,\delta$ of positive real part. This is formula 641.2 on page 655 of Gradshteyn and Ryzhik [GR], which in turn cites formula 302(32) in Erdélyi et al. [E]. [**EDIT** Thanks to Richard Borcherds for recognizing this integral as a result of Barnes 1908. See postscript below about references and proofs for this integral.]

Indeed, the factor $z \phantom./\sinh(2\pi z)$ of the integrand factors as
$$
\frac{1}{2} \frac{iz}{\sin(\pi i z)} \frac1{\cos(\pi i z)} =
\frac1{2\pi^2} \Gamma(1+iz) \phantom. \Gamma(1-iz) \phantom. \Gamma\bigl(\frac12+iz\bigr) \phantom.\Gamma\bigl(\frac12-iz\bigr),
$$
while the factor $\prod_{j=1}^{p-1} (j^2+z^2)$ is
$$
\prod_{j=1}^{p-1} (j+iz) \phantom. \prod_{j=1}^{p-1} (j-iz)
= \frac{\Gamma(p+iz)}{\Gamma(1+iz)} \frac{\Gamma(p-iz)}{\Gamma(1-iz)}.
$$
Hence the integral is
$$
\frac1{2\pi^2} \int_{-\infty}^{\infty} \Gamma(p+iz) \phantom. \Gamma(p-iz) \phantom. \Gamma\bigl(\frac12+iz\bigr) \phantom.\Gamma\bigl(\frac12-iz\bigr) \phantom. dz.
$$
Now take $z=is$ to obtain the known integral, multiplied by $1/2\pi^2 i$.

The formula for that integral, multiplied by the same factor $1/2\pi^2 i$, then yields
$$
\frac1\pi \frac{\Gamma(2p) \phantom. \Gamma(p+(1/2))^2 \phantom. \Gamma(1)}{\Gamma(2p+1)} ,
$$
which is equivalent to the desired answer because $\Gamma(1)=1$ and $\Gamma(2p+1) = 2p \phantom. \Gamma(2p)$.

*Postscript* about the contour integral formula: Gary couldn't locate the Erdélyi reference; I found it, but it turns out the book just displays the formula without proof or source. Fortunately Richard Borcherds (in a comment to this answer) recognized it as a result of Barnes [B]; see Lemma 15 on pages 154-155. There's even a Wikipedia page on Barnes' work in this area, including this result which is called the "first Barnes lemma". Barnes' proof uses a hypergeometric identity of Gauss, which is natural in the context of his paper. Here's an alternative approach using Fourier convolutions. For $\mu,\nu$ of positive real part, and real $t>1$, we have

$$
\int_{-i\infty}^{i\infty} \Gamma(\mu+s) \phantom. \Gamma(\nu-s) \phantom. t^s \phantom. ds
= 2\pi i \phantom. \Gamma(\mu+\nu) \phantom. t^\nu / (1+t)^{\mu+\nu}.
$$
[This is basically [GR, p.657, 6.422#3], and can be proved by shifting the contour integral to the left: the poles of $\Gamma(\mu+s)$ yield the terms in the Laurent expansion of the right-hand side about $t = \infty$.] Now take $(\mu,\nu) = (\alpha,\gamma)$ and $(\delta,\beta)$ and multiply, finding that
$$
\int_{-i\infty}^{i\infty} F(s) \phantom. t^s \phantom. ds = (2\pi i)^2 \Gamma(\alpha+\gamma)\phantom.\Gamma(\beta+\delta)\phantom.t^{\beta+\gamma} / (1+t)^{\alpha+\beta+\gamma+\delta}
$$
where $F(\cdot)$ is the convolution on the imaginary $s$-axis of $\Gamma(\alpha+s) \phantom. \Gamma(\gamma-s)$ with $\Gamma(\delta+s) \phantom. \Gamma(\beta-s)$.
Since this integral uniquely determines $F$, it follows that
$$
F(s) = 2\pi i \frac{\Gamma(\alpha+\gamma) \phantom. \Gamma(\beta+\delta)} {\Gamma(\alpha+\beta+\gamma+\delta)} \Gamma(\alpha+\delta+s) \phantom. \Gamma(\beta+\gamma-s).
$$
Taking $s=0$ recovers the Barnes formula.

**References**

B] Barnes, E.W.: A new development of the theory of the hypergeometric functions. *Proc. LMS* (1908) s2-6(1): 141–177.

[E] Erdélyi, A., et al.: *Table of Integral Transforms* II. New York: McGraw Hill, 1954.

[GR] Gradshteyn, I.S., and Ryzhik, I.M.: *Table of Integrals, Series, and Products* (4th ed.). New York: Academic Press, 1980.