After developing a Resonance Chiral Lagrangian (R?L) model to describe hadronic ? lepton decays, the model was confronted with experimental data. This was accomplished by using a fitting framework that was developed to take into account the complexity of the model and to ensure numerical stability for the algorithms used in the fitting. Since the model used in the fit contained 15 parameters and there were only three one-dimensional distributions available, we could expect multiple local minima or even whole regions of equal potential to appear. Our methods had to thoroughly explore the whole parameter space and ensure (as well as possible) that the result is a global minimum. This paper is focused on the technical aspects of the fitting strategy used. The first approach was based on a re-weighting algorithm published in article Shekhovtsova et al. and produced results in about two weeks. A later approach, with an improved theoretical model and a simple parallelization algorithm based on Inter-Process Communication (IPC) methods of UNIX system, reduced computation time down to 2–3 days. Additional approximations were introduced to the model, decreasing the necessary time to obtain the preliminary results down to 8 hours. This allowed us to better validate the results, leading to a more robust analysis published in article Nugent et al.