$\newcommand{\Ei}{\operatorname{Ei}}
\newcommand{\Ci}{\operatorname{Ci}}
\newcommand{\si}{\operatorname{si}}
\newcommand{\Li}{\operatorname{Li}_4}
$
We have $E(4),E(5)$ in closed forms.
$$\begin{aligned}
&E(2)=\frac{2}{\pi} ,\\
&E(3)=1,\\
&E(4)=\frac{28}{\pi^3}\zeta(3)+\frac{4}{\pi},\\
&E(5)=-\frac{1280}{3\pi^4}\Li\left(\frac12\right)+\frac{235}{27}-\frac{160}{9\pi^4}\ln(2)^4+\frac{160}{9\pi^2}\ln(2)^2,
\end{aligned}
$$
where $\operatorname{Li}_n(z)$ denotes polylogarithms.
Proof. The moment integral, as has shown by another answer, is equivalent to evaluate a sort of integrals with powered trigonometric integral functions.
$$
\begin{aligned}
E(n) & = \left ( \frac{2}{\pi} \right )^n
\int_{[0,\pi/2]^n}\frac{\tan x_1...\tan x_n}{
\tan x_1+...+\tan x_n}\text{d}x_1...\text{d}x_n,\\
&=\left ( \frac{2}{\pi} \right )^n
\int_{\left [ 0,\infty \right ]^n }
\frac{1}{(1+x_1^2)...(1+x_n^2)} \frac{x_1...x_n}{x_1+...+x_n}
\text{d}x_1...\text{d}x_n,\\
&=\left ( -1 \right )^n \left ( \frac{2}{\pi} \right )^n
\int_{0}^{\infty} A(x)^n\text{d}x( A(x) \text{ to be defined as some trigonometric integrals}).
\end{aligned}
$$The last line is by replacing $\frac1a=\int_0^\infty e^{-at}\text{d}t$, switching the order of integration and using
$$
A(\alpha)=-\int_{0}^{\infty} \frac{xe^{-\alpha x}}{1+x^2}\text{d}x
,\quad \alpha\in\mathbb{R}^{+}.$$
We use complex tricks to convert these trigonometric integrals into exponential integrals, thence to the substantial reduction in computation complexity. Denote two exponential analogues(because they behave like exponent functions),
$$
E_1(z)=e^{z}\operatorname{Ei}_1(-z),E_2(z)=e^{-z}\left(\operatorname{Ei}_2(z)-\pi i\right),
$$
where $\Ei_{1,2}(.)$ is the exponent integral function with branch cuts chosen. Use the definition,
$$
\Ei_1(-z)=
\gamma+\log(z)+\int_{0}^{1} \frac{e^{-zt}-1}{t} \text{d}t,\quad
\Ei_2(z)=
\gamma+\log(z)+\int_{0}^{1} \frac{e^{zt}-1}{t} \text{d}t,
$$
and assume $\log(z)=\ln(|z|)+i\arg(z),\arg(z)\in[0,2\pi)$ with branch cut $[0,\infty).$ Note that for a postive real $z$, one have $\Ei_1(-z)=\Ei(-z),\Ei_2(z)=\Ei(z)$.
Next we show that, defining
$$A(x)=\Ci(x)\cos(x)+\si(x)\sin(x),\quad B(x)=\si(x)\cos(x)-\Ci(x)\sin(x),
$$
the following relations are valid for $x\in\mathbb{R}^{+}$. $$
E_1(ix)=A(x)-i B(x),E_2(ix)=A(x)+iB(x).
$$
We start as follows.
$$
\begin{aligned}
\Ei_1(-ix)
&=\gamma+\log(ix)+\int_{0}^{1} \frac{e^{-ixt}-1}{t} \text{d}t,\\
&=\left(\gamma+\log(x)+\int_{0}^{1} \frac{\cos(xt)-1}{t} \text{d}t\right)
+i\left ( \frac\pi2-\int_{0}^{1} \frac{\sin(xt)}{t} \text{d}t\right ),\\
&=\Ci(x)-i \si(x).\\
\Rightarrow E_1(-ix)&=e^{ix}E_1(-ix),\\
&=\left ( \cos(x)+i \sin(x) \right ) \left(\Ci(x)-i \si(x)\right),\\
&=A(x)-iB(x).
\end{aligned}
$$The other can be done simultaneously.
And it's worth noticing that for a positive real $\alpha$
$$
A(\alpha)=-\int_{0}^{\infty} \frac{xe^{-\alpha x}}{1+x^2}\text{d}x
,B(\alpha)=-\int_{0}^{\infty} \frac{e^{-\alpha x}}{1+x^2}\text{d}x.
$$
A short proof is given by integrating respectively
$$
f(z)=\frac{ze^{-\alpha z}\left(\Ei_2(\alpha z)-\pi i\right)}{1+z^2},\frac{e^{-\alpha z}\left(\Ei_2(\alpha z)-\pi i\right)}{1+z^2}
$$
along the real-axis. Compare the imaginary part, giving
$$
\begin{aligned}
& \int_{0}^{\infty} \frac{xe^{-\alpha x}}{1+x^2}\text{d}x
=\left ( -\frac{1}{\pi} \right ) \Im\left ( 2\pi i\operatorname{Res} (f(z),i)\right )
=-\Re E_2(i\alpha)=-A(\alpha),\\
&\int_{0}^{\infty} \frac{e^{-\alpha x}}{1+x^2}\text{d}x
=\left ( -\frac{1}{\pi} \right ) \Im\left ( 2\pi i\operatorname{Res} (f(z),i)\right )
=-\Im E_2(i\alpha)=-B(\alpha).
\end{aligned}
$$
Then we apply some simple contour integration. Define the contour enclosed in the first quadrant with $R$ sufficiently large: $$C:[0,R)\cup\left \{ Re^{i\theta} \mid
\theta\in\left[0,\frac{\pi}{2}\right]\right \}
\cup\left \{ ix\mid x\in[0,R) \right \},$$ and integrate the function
$$
f(z)=E_1(z)^mE_2(z)^n,m,n\in\mathbb{N},m+n\ge2,
$$
giving
$$
\int_{0}^{\infty}e^{(m-n)x}\Ei(-x)^m\left ( \Ei(x)-\pi i \right )^n
\text{d} x
=i\int_{0}^{\infty}\left ( A(x)-i B(x) \right )^m
\left ( A(x)+iB(x) \right )^n\text{d}x.
$$
The LHS are, by definition, some common logarithmical integrals.(If interested, click here.)
To calculate $E(4)$, substituting $(m,n)=(4,0),(3,1),(2,2)$ to obtain
$$
\begin{aligned}
&\int_{0}^{\infty}e^{4x}\Ei(-x)^4\text{d} x
=i\int_{0}^{\infty}\left ( A(x)-i B(x) \right )^4\text{d}x,\\
&\int_{0}^{\infty}e^{2x}\Ei(-x)^3\left ( \Ei(x)-\pi i \right )
\text{d} x
=i\int_{0}^{\infty}\left ( A(x)-i B(x) \right )^2
\left ( A(x)^2+B(x)^2\right )\text{d}x,\\
&\int_{0}^{\infty}\Ei(-x)^2\left ( \Ei(x)-\pi i \right )^2
\text{d} x
=i\int_0^\infty\left ( A(x)^2+B(x)^2 \right )^2\text{d}x.
\end{aligned}
$$
Discard the real part of every equality, solve the linear system and make use of these known evaluations,
$$
\int_{0}^{\infty}e^{2x}\Ei(-x)^3\text{d}x = -\frac72 \zeta(3),\qquad\int_{0}^{\infty}\Ei(-x)^2\Ei(x)\text{d}x =- \frac{\pi^2}{3}.
$$
We obtain
$$
\color{blue}{\int_{0}^{\infty} A(x)^4\text{d}x
=\frac{7\pi}{4}\zeta(3)+\frac{\pi^3}4},\qquad
\int_{0}^{\infty} A(x)^2B(x)^2\text{d}x
=\frac{\pi^3}{12},\qquad \int_{0}^{\infty} B(x)^4\text{d}x
=-\frac{7\pi}{4}\zeta(3)+\frac{\pi^3}4.
$$
The first integral is exactly $(\frac\pi2)^4E(4)$, which completes the proof.
To calculate $E(5)$, substituting $(m,n)=(5,0),(4,1),(3,2)$ to obtain
$$
\begin{aligned}
&\int_{0}^{\infty}e^{5x}\Ei(-x)^5\text{d} x
=i\int_{0}^{\infty}\left ( A(x)-i B(x) \right )^5\text{d}x,\\
&\int_{0}^{\infty}e^{3x}\Ei(-x)^4\left ( \Ei(x)-\pi i \right )
\text{d} x
=i\int_{0}^{\infty}\left ( A(x)-i B(x) \right )^3
\left ( A(x)^2+B(x)^2\right )\text{d}x,\\
&\int_{0}^{\infty}e^{x}\Ei(-x)^3\left ( \Ei(x)-\pi i \right )^2
\text{d} x
=i\int_{0}^{\infty}\left ( A(x)-i B(x) \right )
\left ( A(x)^2+B(x)^2 \right )^2\text{d}x.
\end{aligned}
$$
Discard the real part of every equality, solve the linear system and make use of these known evaluations,
$$
\begin{aligned}
&\int_{0}^{\infty}e^{3x}\Ei(-x)^4\text{d}x = \frac{26\pi^4}{135} +\frac{4\pi^2}{9}\ln(2)^2
-\frac{4}{9}\ln(2)^4-\frac{32}{3}\Li\left ( \frac12 \right ) ,\\
&\int_{0}^{\infty}e^{x}\Ei(-x)^3\Ei(x)\text{d}x = \frac{61\pi^4}{360} +\frac{\pi^2}{3}\ln(2)^2
-\frac{1}{3}\ln(2)^4-8\Li\left ( \frac12 \right ).
\end{aligned}
$$
We obtain
$$
\begin{aligned}
&\color{purple}{\int_{0}^{\infty} A(x)^5\text{d}x
=\frac{40\pi}{3}\,\text{Li}_4\left(\frac{1}{2}\right)-\frac{235\pi^5}{864}+\frac{5\pi}{9}\ln(2)^4-\frac{5\pi^3}{9}\ln(2)^2},\\
&\int_{0}^{\infty} A(x)^3B(x)^2\text{d}x
=\frac{4\pi}{3}\,\text{Li}_4\left(\frac{1}{2}\right)-\frac{131\pi^5}{4320}+\frac{\pi}{18}\ln(2)^4-\frac{\pi^3}{18}\ln(2)^2,\\
&\int_{0}^{\infty} A(x)B(x)^4\text{d}x
=-\frac{\pi^5}{160} .
\end{aligned}
$$
The first integral is exactly $-(\frac\pi2)^5E(5)$, which completes the proof.
Remark. It would be exhaustive to evaluate each integral separately. Moreover, the result of the last integral shows great simplicity. In fact, using integral representations we have
$$
\int_{0}^{\infty} A(x)B(x)^4\text{d}x
=-\int_{[0,\infty]^5}
\frac{1}{\left ( 1+x_1^2\right )
\left ( 1+x_2^2\right )\left ( 1+x_3^2\right )
\left ( 1+x_4^2\right )\left ( 1+x_5^2\right ) }
\frac{x_1}{x_1+x_2+x_3+x_4+x_5}\text{d}x_1...\text{d}x_5=-\frac{1}{5} \left ( \frac{\pi}{2} \right )^5
$$
by symmetry.
