A model suited for calculating correlation functions in QCD from the ultraviolet to the infrared is reviewed. The model consist in standard Faddeev-Popov Lagrangian for Landau gauge with an extra mass term for gluons. It is shown that once this mass term is included, two and three point correlation functions can be calculated with good precision at one-loop order even at very low momenta in the quenched approximation. After that, the inclusion of quarks is analyzed. It is shown that in that case, the perturbative calculation only gives the qualitative behavior in some correlation functions. In particular, it is shown that the analysis of spontaneous chiral symmetry breaking requires to go beyond perturbation theory. A non-perturbative scheme but controlled by two small parameters is discussed and the corresponding results are shown to agree with high precision to Monte Carlo numerical simulations also in the quark sector