In this talk, we introduce our recently completed next-to-leading order (NLO) global analysis of the nuclear parton distribution functions (nPDFs) called EPS09 - a higher order successor to the well-known leading-order (LO) analysis EKS98 and also to our previous LO work EPS08. As an extension to similar global analyses carried out by other groups, we complement the data from deep inelastic $l+A$ scattering and Drell-Yan dilepton measurements in p+$A$ collisions by inclusive midrapidity pion production data from d+Au collisions at RHIC, which results in better constrained gluon distributions than before. The most important new ingredient, however, is the detailed error analysis, which employs the Hessian method and which allows us to map out the parameter-space vicinity of the best-fit to a collection of nPDF error sets. These error sets provide the end-user a way to compute how the PDF-uncertainties will propagate into the cross sections of his/her interest. The EPS09 package to be released soon, will contain both the NLO and LO results for the best fits and the uncertainty sets.