  Multiple Mellin-Barnes integrals are often used for perturbative calculations
in particle physics. In this context, the evaluation of such objects may be
performed through residues calculations which lead to their expression as
multiple series in powers and logarithms of the parameters involved in the
problem under consideration. However, in most of the cases, several series
representations exist for a given integral. They converge in different regions
of values of the parameters, and it is not obvious to obtain them. For twofold
integrals we present a method which allows to derive straightforwardly and
systematically: (a) different sets of poles which correspond to different
convergent double series representations of a given integral, (b) the regions
of convergence of all these series (without an a priori full knowledge of their
general term), and (c) the general term of each series (this may be performed,
if necessary, once the relevant domain of convergence has been found). This
systematic procedure is illustrated with some integrals which appear, among
others, in the calculation of the two-loop hexagon Wilson loop in N = 4 SYM
theory. Mellin-Barnes integrals of higher dimension are also considered.
