We study the covariance matrix of the cluster mass function in cosmology . We adopt a two-line attack : firstly , we employ the counts-in-cells framework to derive an analytic expression for the covariance of the mass function . Secondly , we use a large ensemble of N -body simulations in the \Lambda CDM framework to test this . Our theoretical results show that the covariance can be written as the sum of two terms : a Poisson term , which dominates in the limit of rare clusters ; and a sample variance term , which dominates for more abundant clusters . Our expressions are analogous to those of for multiple cells and a single mass tracer . Calculating the covariance depends on : the mass function and bias of clusters , and the variance of mass fluctuations within the survey volume . The predictions show that there is a strong bin-to-bin covariance between measurements . In terms of the cross-correlation coefficient , we find r \gtrsim 0.5 for haloes with M \lesssim 3 \times 10 ^ { 14 } h ^ { -1 } M _ { \odot } at z = 0 . Comparison of these predictions with estimates from simulations shows excellent agreement . We use the Fisher matrix formalism to explore the cosmological information content of the counts . We compare the Poisson likelihood model , with the more realistic likelihood model of , and all terms entering the Fisher matrices are evaluated using the simulations . We find that the Poisson approximation should only be used for the rarest objects , M \gtrsim 3 \times 10 ^ { 14 } h ^ { -1 } M _ { \odot } , otherwise the information content of a survey of size V \sim 13.5 h ^ { -3 } { Gpc } ^ { 3 } would be overestimated , resulting in errors that are \sim 2 times smaller . As an auxiliary result , we show that the bias of clusters , obtained from the cluster-mass cross-variance , is linear on scales > 50 h ^ { -1 } { Mpc } , whereas that obtained from the auto-variance is nonlinear .