  We present a next-to-leading order (NLO) global DGLAP analysis of nuclear
parton distribution functions (nPDFs) and their uncertainties. Carrying out an
NLO nPDF analysis for the first time with three different types of experimental
input -- deep inelastic $\ell$+A scattering, Drell-Yan dilepton production in
p+$A$ collisions, and inclusive pion production in d+Au and p+p collisions at
RHIC -- we find that these data can well be described in a conventional
collinear factorization framework. Although the pion production has not been
traditionally included in the global analyses, we find that the shape of the
nuclear modification factor $R_{\rm dAu}$ of the pion $p_T$-spectrum at
midrapidity retains sensitivity to the gluon distributions, providing evidence
for shadowing and EMC-effect in the nuclear gluons. We use the Hessian method
to quantify the nPDF uncertainties which originate from the uncertainties in
the data. In this method the sensitivity of $\chi^2$ to the variations of the
fitting parameters is mapped out to orthogonal error sets which provide a
user-friendly way to calculate how the nPDF uncertainties propagate to any
factorizable nuclear cross-section. The obtained NLO and LO nPDFs and the
corresponding error sets are collected in our new release called {\ttfamily
EPS09}. These results should find applications in precision analyses of the
signatures and properties of QCD matter at the LHC and RHIC.
