In Billiards in Rectangles with Barriers, Eskin-Masur-Schmoll count the number of primitive square-tiled surfaces with two cone points with angle $4\pi$ for one cone point with angle $6\pi$: (See also):

$$ \begin{array}{crlr} N_d^P(1,1) &=& \frac{1}{3} d^3(d-1) \prod_{p|d} \left( 1 - \frac{1}{p^2}\right) \\ N_d^P(2) &=& \frac{3}{8} d^2(d-2) \prod_{p|d} \left( 1 - \frac{1}{p^2}\right) \end{array}$$

Where $\sum_{r|d} \frac{\mu(r)}{r^2} = \prod_{p|d} \left( 1 - \frac{1}{p^2}\right) $

It is mentioned these exact formulas can be derived from The Character of the Infinite Wedge by Bloch and Okounkov.

Their main result is stated very generally thru an n-point function, related to the theta function

$$ \Theta(t):= \eta(q)^{-3}\sum_{n \in \mathbb{Z}}(-1)^n q^{\frac{(n + \frac{1}{2})^2}{2} } x^{n + \frac{1}{2} }= (q)_\infty^{-2}\left(\sqrt{x} - \frac{1}{\sqrt{x}}\right)(qx)_\infty (q/x)_\infty$$

The original derivation in EMS is rather oblique using the Siegel-Veech formula

\[ \frac{1}{\tilde{\mu}(\mathcal{M}(S))}\int_{\mathcal{M}(S)}\hat{f}d\tilde{\mu} = k(S) \int_{\mathbb{R}^2} f \]

Are there derivations of this fact, possibly without resorting to moduli spaces or vertex operators that make the product structure of these counts particularly clear?

I am guessing the product over primes has to do with the primitive condition. In fact, I've asked about these functions before (and also here)