In our case we have as the Monte Carlo sampling function the probability for finding the system in a state \( s \) given by $$ \begin{equation*} P_s=\frac{e^{-(\beta E_s)}}{Z}, \end{equation*} $$ with energy \( E_s \), \( \beta=1/kT \) and \( Z \) is a normalization constant which defines the partition function in the canonical ensemble. As discussed above $$ \begin{equation*} Z(\beta)=\sum_se^{-(\beta E_s)} \end{equation*} $$ is difficult to compute since we need all states.