Let S be a closed, orientable 2d manifold and G a finite group. Since a principal G-bundle over S is specified by maps $\phi : \pi_1(S) \rightarrow G$ modulo the adjoint action by G, the way to count how many such bundles exist is 
\begin{align}
\dfrac{ | Hom(\pi_1(S),G) |}{|G|} = \sum_{\text{irreps V of }G} \bigg( \dfrac{d(V)}{|G|}  \bigg)^{2-2g} \; ,
\end{align}
where $d(V)$ is the dimension of the irrep V, and g is the genus of S. 
My question is how the formula above generalizes when we want to count principal G-bundles over S, when S has n boundaries with prescribed holonomies for the bundle specified by conjugacy classes $\{ k_i : i=1,...,n \}$ of G. Morally we should have something like
\begin{align}
\dfrac{ | Hom(\pi_1(S_{g,k_i}),G) |}{|G|} = \sum_{\text{irreps V of }G} \bigg( \dfrac{d(V)}{|G|}  \bigg)^{2-2g-n} \prod_{i=1}^n \chi_q(k_i)  \; ,
\end{align}
where $\chi_q(k_i)$ is the character of the representation q evaluated at any representative of the conjugacy class $k_i$. 
So, is the above formula completely correct? Are there factors missing? I would really appreciate any references that prove what the generalization is.