  We present an improved leading-order global DGLAP analysis of nuclear parton
distribution functions (nPDFs), supplementing the traditionally used data from
deep inelastic lepton-nucleus scattering and Drell-Yan dilepton production in
proton-nucleus collisions, with inclusive high-$p_T$ hadron production data
measured at RHIC in d+Au collisions. With the help of an extended definition of
the $\chi^2$ function, we now can more efficiently exploit the constraints the
different data sets offer, for gluon shadowing in particular, and account for
the overall data normalization uncertainties during the automated $\chi^2$
minimization. The very good simultaneous fit to the nuclear hard process data
used demonstrates the feasibility of a universal set of nPDFs, but also
limitations become visible. The high-$p_T$ forward-rapidity hadron data of
BRAHMS add a new crucial constraint into the analysis by offering a direct
probe for the nuclear gluon distributions -- a sector in the nPDFs which has
traditionally been very badly constrained. We obtain a strikingly stronger
gluon shadowing than what has been estimated in previous global analyses. The
obtained nPDFs are released as a parametrization called EPS08.
