  Both collision geometry and event-by-event fluctuations are encoded in the
experimentally observed flow harmonic distribution $p(v_n)$ and $2k$-particle
cumulants $c_n\{2k\}$. In the present study, we systematically connect these
observables to each other by employing Gram-Charlier A series. We quantify the
deviation of $p(v_n)$ from Bessel-Gaussianity in terms of flow harmonic
fine-splitting. Subsequently, we show that the corrected Bessel-Gaussian
distribution can fit the simulated data better than the Bessel-Gaussian
distribution in the more peripheral collisions. Inspired by Gram-Charlier A
series, we introduce a new set of cumulants $q_n\{2k\}$ that are more natural
to study distributions near Bessel-Gaussian. These new cumulants are obtained
from $c_n\{2k\}$ where the collision geometry effect is extracted from it. By
exploiting $q_2\{2k\}$, we introduce a new set of estimators for averaged
ellipticity $\bar{v}_2$ which are more accurate compared to $v_2\{2k\}$ for
$k>1$. As another application of $q_2\{2k\}$, we show we are able to restrict
the phase space of $v_2\{4\}$, $v_2\{6\}$ and $v_2\{8\}$ by demanding the
consistency of $\bar{v}_2$ and $v_2\{2k\}$ with $q_2\{2k\}$ equation. The
allowed phase space is a region such that $v_2\{4\}-v_2\{6\}\gtrsim 0$ and $12
v_2\{6\}-11v_2\{8\}-v_2\{4\}\gtrsim 0$, which is compatible with the
experimental observations.
