Here is a method that will allow one to find the exact upper and lower bounds on $g(z)$ over $z>0$ with any degree of accuracy.
Take any real $z>0$. Since
\begin{equation*}
\frac1y=\int_0^\infty dt\,e^{-y t}
\end{equation*}
for any real $y>0$, we have
\begin{align*}
\frac{g(z)}z
&=\int_z^{e z} dy\, \frac{\sin y}y \\
&=\int_0^\infty dt\,\int_z^{e z} dy\,e^{-y t}\sin y \\
&=\int_0^\infty dt\,
\Big( \frac{e^{-t z} (\cos z+t \sin z)}{t^2+1}
-\frac{e^{-e t z} (\cos ez+t \sin ez)}{t^2+1}\Big) \\
&=I_1(z) \cos z+I_2(z)\sin z -I_1(ez) \cos ez-I_2(ez)\sin ez, \tag{1}
\end{align*}
where
\begin{align*}
I_1(z)&:=\int_0^\infty dt\,\frac{e^{-t z}}{t^2+1}, \\
I_2(z)&:=\int_0^\infty dt\,\frac{e^{-t z}t}{t^2+1}.
\end{align*}
Next, letting $c_1$ and $c_2$ denote functions with values in $(0,1)$, we have
\begin{align*}
I_1(z)&=\frac1z\,\int_0^\infty du\,\frac{e^{-u}}{1+u^2/z^2} \\
&=\frac1z\,\int_0^\infty du\,e^{-u}
-\frac1z\,\int_0^\infty du\,\frac{u^2e^{-u}}{z^2+u^2} \\
&=\frac1z-\frac{2c_1(z)}{z^3};
\end{align*}
at the last step here, we used the inequality $z^2+u^2>z^2$ for $u>0$;
similarly,
\begin{align*}
I_2(z)&=\frac1{z^2}-\frac{3c_2(z)}{z^4}.
\end{align*}
So, by (1),
\begin{equation*}
g=h+r,
\end{equation*}
where
\begin{equation*}
h(z):=\cos z-\tfrac1e\,\cos ez
\end{equation*}
and
\begin{equation*}
r(z):=-\frac{2c_1(z)}{z^2}\, \cos z-\frac{3c_2(z)}{z^3}\,\sin z
+\frac{2c_1(ez)}{e^3z^2}\, \cos ez+\frac{3c_2(2z)}{e^4z^3}\,\sin ez
\end{equation*}
is the "remainder",
so that
\begin{equation*}
|r(z)|<\frac{2.1}{z^2}+\frac{3.1}{z^3},
\end{equation*}
which can be made however small if $z$ is large enough.
On the other hand, since $e$ is irrational, we will have \begin{equation*} \sup_{z>0}h(z)=-\inf_{z>0}h(z)=1+1/e=1.367\dots \end{equation*} (which is somewhat close to your value $1.4$).
So, to compute $\sup_{z>0}g(z)$ and $\inf_{z>0}g(z)$ with any degree of accuracy, it suffices to be able to compute $\sup_{z\in(0,a]}g(z)$ and $\inf_{z\in(0,a]}g(z)$ with any degree of accuracy for any given real $a>0$, which can be done by (say) the interval arithmetic method, using the formula $g(z)=z(\text{Si}(e z)-\text{Si}(z))$ and the monotonicity of the function $\text{Si}$ on each of the intervals of the form $[k\pi,(k+1)\pi]$ for $k=0,1,\dots$.