  We discuss analytic continuation as a tool to extract the cumulants of the
quark number fluctuations in the strongly interacting medium from lattice QCD
simulations at imaginary chemical potentials. The method is applied to $N_f =
2+1$ QCD, discretized with stout improved staggered fermions, physical quark
masses and the tree level Symanzik gauge action, exploring temperatures ranging
from 135 up to 350 MeV and adopting mostly lattices with $N_t = 8$ sites in the
temporal direction. The method is based on a global fit of various cumulants as
a function of the imaginary chemical potentials. We show that it is
particularly convenient to consider cumulants up to order two, and that below
$T_c$ the method can be advantageous, with respect to a direct Montecarlo
sampling at $\mu = 0$, for the determination of generalized susceptibilities of
order four or higher, and especially for mixed susceptibilities, for which the
gain is well above one order of magnitude. We provide cumulants up to order
eight, which are then used to discuss the radius of convergence of the Taylor
expansion and the possible location of the second order critical point at real
$\mu$: no evidence for such a point is found in the explored range of $T$ and
for chemical potentials within present determinations of the pseudo-critical
line.
