Using the recent world average alpha(s)(M(Z)(2)) = 0.118 +/- 0.006, we give the first direct extraction from the Psi and Upsilon data of the values of the running heavy quark masses within QCD spectral sum rules to two-loops in the <(MS)over bar>-scheme: (m) over bar(b)(M(b)(PT2)) = (4.23(-0.04)(+0.03) +/- 0.02) GeV and (m) over bar(c)(M(c)(PT2)) = (1.23(-0.04)(+0.03) +/- 0.03) GeV, (the errors are respectively due to alpha(s) and to the gluon condensate),and the corresponding value of the short-distance perturbative pole masses to two-loops: M(b)(PT2) = (4.62 +/- 0.02) GeV, M(c)(PT2) = (1.42 +/- 0.03) GeV, which we compare with the updated values of the non-relativistic pole masses re-extracted directly from the two-loop non-relativistic sum rules: M(b)(NR) = (4.69(-0.01)(+0.02) +/- 0.02) GeV and M(c)(NR) = (1.45(-0.03)(+0.04) +/- 0.03) GeV. It is also informative to compare the three-loop values of the short-distance pole masses: M(b)(PT3) = (4.87 +/- 0.05 +/- 0.02) GeV and M(c)(PT3) = (1.64(-0.07)(+0.10) +/- 0.03) GeV, with the dressed mass M(b)(nr) = (4.94 +/- 0.10 +/- 0.03) GeV, entering into the non-relativistic Balmer formula including higher order alpha(s) corrections. The small mass-differences M(b)(NR) - M(b)(PT2) similar or equal to M(b)(nr) - M(b)(PT3) similar or equal to 70 MeV and M(c)(NR) - M(c)(PT2) similar or equal to (30 +/- 30) MeV can measure the size of the non-perturbative effect induced by renormalon type-singularities. Finally, the b and c quark-pole mass difference is found to delta M(bc) = M(b)-M(c) = (3.22 +/- 0.03) GeV. An analogous analysis is pursued for the heavy-light mesons, where a simultaneous re-fit of the B and B* masses from relativistic sum rules leads to: M(b)(PT2) = (4.63 +/- 0.08) GeV, while the average of the results from full-QCD and HQET sum rules in the large mass limit gives the meson-quark mass difference to two-loops: delta M(b)(infinity) = Lambda = (M(B) - M(b)(NR))(infinity) similar or equal to (0.58 +/- 0.05) GeV. A comparison of these new and accurate results with the existing ones in the literature is done. As a consequence, the updated values of the pseudoscalar decay constants to two-loops are: f(D) = (1.35 +/- 0.04 +/- 0.06) f(pi) and f(B) = (1.49 +/- 0.06 +/- 0.05) f(pi), which lead to f(B) root B-B = (149 +/- 0.14)f(pi).