$\newcommand{\expn}{\mathbb{E}}
\newcommand{\varn}{\operatorname{Var}}$
Following up on Brendan McKay's approach in the comments, let $v, w$ be two uniformly random vectors on the sphere. We are interested in the following probability:
\begin{align}
p = P(v > 0 | w > 0, v \cdot w = \alpha) = \frac{P(v, w > 0 | v \cdot w = \alpha)}{P(w > 0 | v \cdot w = \alpha)} = 2^d P(v, w > 0 | v \cdot w = \alpha).
\end{align}
Here we used that the events $\{w > 0\}$ and $\{v \cdot w = \alpha\}$ are pairwise independent. Bayes' law states that
\begin{align}
P(v, w > 0 | v \cdot w = \alpha) = \frac{P(v, w > 0)}{P(v \cdot w = \alpha)} \, P(v \cdot w = \alpha | v, w > 0).
\end{align}
Since clearly $P(v > 0) = P(w > 0) = 2^{-d}$, and $\{v > 0\}$, $\{w > 0\}$, $\{v \cdot w = \alpha\}$ are pairwise independent events, this is equal to
\begin{align}
p &= \frac{2^d P(v, w > 0)}{P(v \cdot w = \alpha)} \, P(v \cdot w = \alpha | v, w > 0) = \frac{P(v \cdot w = \alpha | v, w > 0)}{2^d P(v \cdot w = \alpha)}.
\end{align}
Note that for two vectors on the sphere, the density $P(v \cdot w = \alpha)$ is well-known to be proportional to $\sin(\arccos \alpha)^d = (1 - \alpha^2)^{d/2}$ by computing the volume of the corresponding spherical cap.
For computing $P(v \cdot w = \alpha | v, w > 0)$, we use a Gaussian approximation: a $d$-dimensional vector $v$ or $w$ on the unit sphere (in the positive orthant) roughly follows the same distribution as sampling a $d$-dimensional vector where each entry is drawn from a(n absolute) Gaussian distribution with mean $0$ and variance $\frac{1}{d}$. This is the approximation which makes the end result possibly incorrect.
So, for computing $P(v \cdot w = \alpha | v, w > 0)$, we basically want to compute the density of inner products between $v$ and $w$ at $\alpha$, given that both $v$ and $w$ are sampled from an absolute Gaussian distribution $|N(0,\frac{1}{d})|^d$ (also known as the half-normal distribution). The inner product is described by $S = \sum_{i=1}^d V_i W_i$ where the $2d$ random variables $V_i, W_i \sim |N(0, \frac{1}{d})|$ are all independent, and as each product $X_i = V_i W_i$ follows the same distribution (and is independent from other values $X_j$), by the CLT the sum $S$ of these $d$ i.i.d.\ random variables is approximately normally distributed with a certain mean and variance. As $\expn(V_i) = \expn(W_i) = \sqrt{2/(\pi d)}$ and $\expn(V_i^2) = \expn(W_i^2) = \frac{1}{d}$, for $S$ we obtain the following expectation
\begin{align}
\expn(S) = d \cdot \expn(V_i W_i) = d \cdot \expn(V_i) \cdot \expn(W_i) = \frac{2}{\pi}.
\end{align}
Similarly, for the variance we get
\begin{align}
\varn(S) = d \cdot \left[\expn(V_i^2 W_i^2) - \expn(V_i W_i)^2\right] &= \frac{1}{d} \left(1 - \frac{4}{\pi^2}\right).
\end{align}
As a result, the probability $P(v \cdot w = \alpha | v, w > 0)$ roughly corresponds to the density of a Gaussian random variable $Z$ with mean $2/\pi$ and variance $(1 - 4/\pi^2) / d$ at the point $\alpha$. This density is equal to
\begin{align}
P(v \cdot w = \alpha | v, w > 0) \approx P\left(N\left(\frac{2}{\pi}, \frac{1}{d} - \frac{4}{\pi^2 d}\right) = \alpha\right) \propto \exp\left(-\frac{(\alpha - \frac{2}{\pi})^2}{\frac{2}{d} (1 - \frac{4}{\pi^2})}\right) = \exp\left(-\frac{d (\alpha \pi - 2)^2}{2 \pi^2 - 8}\right).
\end{align}
To summarize, by using exact arguments for $P(v \cdot w = \alpha)$ (i.e.\ not using a Gaussian approximation), and using approximate arguments for $P(v \cdot w = \alpha | v, w > 0)$, we obtain the end result:
\begin{align}
\boxed{\ p \approx \exp\left[d \left(- \ln 2 - \frac{1}{2} \ln(1 - \alpha^2) - \frac{(\alpha \pi - 2)^2}{2 \pi^2 - 8} \right) + o(d)\right].}
\end{align}
Substituting $\alpha = \frac{1}{2}$, this results in
\begin{align}
g\left(\frac{\pi}{3}\right) \approx \frac{\sqrt{3}}{3} \, \exp\left(-\frac{(\pi -4)^2}{8 \pi^2 - 32}\right) \approx {\color{blue}{0.5684}}.
\end{align}
Alternative approach
By also using a Gaussian approximation for the density $P(v \cdot w = \alpha)$, instead of using exact arguments there we can follow a similar path as for the probability $P(v \cdot w = \alpha | v, w > 0)$ by approximating the spherical vectors with Gaussian vectors, and invoking the CLT to approximate the probabilities in the limit for large $d$. Similar computations show that in that case the inner product has mean $0$ and variance $\frac{1}{d}$, and so the probability $P(v \cdot w = \alpha)$ can be estimated as
\begin{align}
P(v \cdot w = \alpha) \approx P\left(N\left(0, \tfrac{1}{d}\right) = \alpha\right) \propto \exp\left(-\frac{\alpha^2 d}{2}\right).
\end{align}
In that case, using a Gaussian approximation for both densities, we obtain the approximation
\begin{align}
\boxed{\ p \approx \exp\left[d \left(- \ln 2 + \frac{\alpha^2}{2} - \frac{(\alpha \pi - 2)^2}{2 \pi^2 - 8} \right) + o(d)\right].}
\end{align}
Substituting $\alpha = \frac{1}{2}$ here, we obtain
\begin{align}
g\left(\frac{\pi}{3}\right) \approx \frac{1}{2} \, \exp\left(\frac{2 \pi - 5}{2 \pi^2 - 8}\right) \approx {\color{blue}{0.5578}}.
\end{align}
Concluding, together with Carlo's result, this gives three different answers (all Gaussian approximations) to the same question of computing $g(\pi/3)$. It is not clear which of them is correct, or whether they are perhaps all wrong...
Update: Now community wiki -- please feel free to update the answer with correct arguments, not using the CLT.
Using large deviations
Instead of using the CLT for $P(v \cdot w = \alpha | v, w > 0)$, let us use Cramer's theorem. Let $S = \frac{1}{d} \sum_{i=1}^d \tilde{V}_i \tilde{W}_i$
be the sum over the products of normalized half-normal random variables $\tilde{V}_i, \tilde{W}_i \sim |N(0,1)|$. Then as computed above, for
$X_i = \tilde{V}_i \tilde{W}_i$ we have $\expn(X_i) = 2/\pi$ and $\varn(X_i) = 1 - 4/\pi^2$. With some effort we can further compute the MGF of $X_1$ as
\begin{align}
\expn(e^{t X_1}) &= \frac{2}{\pi} \int_0^{\infty} \int_0^{\infty} \exp(\tfrac{-v^2 - w^2}{2} + t v w) \ dv \, dw = \frac{\pi + 2 \arcsin t}{\pi \sqrt{1 - t^2}}.
\end{align}
This last expression is conditioned on $t < 1$, as otherwise the integral does not converge. Using Cramer's theorem, we have that $P(S > \alpha)$ is
proportional to $\exp(-d I(\alpha))$ where $I(z)$ is the function
\begin{align}
I(z) = \sup_{t \in \mathbb{R}_{+}} \left[t z - \ln \expn(e^{t X_1})\right].
\end{align}
Substituting the above, this is equal to
\begin{align}
I(z) = \sup_{t \in (0,1)} \left[t z - \ln \left(\frac{\pi + 2 \arcsin t}{\pi \sqrt{1 - t^2}}\right)\right].
\end{align}
The function of $t$ on the right is decreasing on $(0,\frac{2}{\pi})$ so the supremum is attained at $t \to 0^+$, in which case the limit is also $0$. This leads
to $I(z) = 0$ and so
\begin{align}
P(v \cdot w > \alpha | v, w > 0) &\sim \exp(-d \cdot 0) = 1.
\end{align}
The problem in this derivation is that Cramer's theorem only applies when $\expn(e^{t X_1}) < \infty$ for all $t \in \mathbb{R}_+$ which is not the case here. So this first try does not work.
