The shapes of the inclusive photon spectra in the processes $J/\stackrel{\ensuremath{\rightarrow}}{\ensuremath{\psi}}\ensuremath{\gamma}X$ and $\stackrel{\ensuremath{\rightarrow}}{\ensuremath{\Upsilon}}\ensuremath{\gamma}X$ have been analyzed using all available experimental data. Relativistic, higher order QCD, and gluon mass corrections were taken into account in the fitted functions. Only on including the gluon mass corrections were consistent and acceptable fits obtained. Values of ${0.721}_{\ensuremath{-}0.068}^{+0.016} \mathrm{GeV}$ and ${1.18}_{\ensuremath{-}0.29}^{+0.09} \mathrm{GeV}$ were found for the effective gluon masses (corresponding to Born level diagrams) for $J/\ensuremath{\psi}$ and $\ensuremath{\Upsilon},$ respectively. The width ratios $\ensuremath{\Gamma}(\stackrel{\ensuremath{\rightarrow}}{V}\mathrm{hadrons})/\ensuremath{\Gamma}(\stackrel{\ensuremath{\rightarrow}}{V}\ensuremath{\gamma}+\mathrm{hadrons})V=J/\ensuremath{\psi},\ensuremath{\Upsilon},$ were used to determine ${\ensuremath{\alpha}}_{s}(1.5 \mathrm{GeV})$ and ${\ensuremath{\alpha}}_{s}(4.9 \mathrm{GeV}).$ Values consistent with the current world average ${\ensuremath{\alpha}}_{s}$ were obtained only when gluon mass correction factors, calculated using the fitted values of the effective gluon mass, were applied. A gluon mass $\ensuremath{\simeq}1 \mathrm{GeV},$ as suggested by these results, is consistent with previous analytical theoretical calculations and independent phenomenological estimates, as well as with a recent, more accurate, lattice calculation of the gluon propagator in the infrared region.