We suggest a framework based on the rainbow approximation with effective parameters adjusted to lattice data. The analytic structure of the gluon and ghost propagators of QCD in Landau gauge is analyzed by means of numerical solutions of the coupled system of truncated Dyson-Schwinger equations. We find that the gluon and ghost dressing functions are singular in complex Euclidean space with singularities as isolated pairwise conjugated poles. These poles hamper solving numerically the Bethe-Salpeter equation for glueballs as bound states of two interacting dressed gluons. Nevertheless, we argue that, by knowing the position of the poles and their residues, a reliable algorithm for numerical solving the Bethe-Salpeter equation can be established.