  We analyse the probability distribution functions of the energy and
magnetisation of the two-dimensional Blume-Capel and Baxter-Wu models with spin
values $S \in \{1/2,1, 3/2\}$ in the presence of a crystal field $\Delta$. By
employing extensive single-spin flip Monte Carlo simulations and a recently
developed method of studying the zeros of the energy probability distribution
we are able to probe, with a good numerical accuracy, several critical
characteristics of the transitions. Additionally, the universal aspects of
these transitions are scrutinised by computing the corresponding probability
distribution functions. The energy distribution has been underutilised in the
literature when compared to that of the magnetisation. Somewhat surprisingly,
however, the former appears to be more robust in characterising the
universality class for both models upon varying the crystal field $\Delta$ than
the latter. Finally, our analysis suggests that in contrast to the Blume-Capel
ferromagnet, the Baxter-Wu model appears to suffer from strong finite-size
effects, especially upon increasing $\Delta$ and $S$, that obscure the
application of traditional finite-size scaling approaches.
