The paper presents a global analysis of nuclear parton distribution functions (nPDFs) at next-to-leading order (NLO) accuracy in perturbative QCD (pQCD). The analysis includes three types of experimental data: deep inelastic lepton-nucleus scattering (DIS), Drell-Yan dilepton production in p+A collisions, and inclusive pion production in d+Au and p+p collisions at RHIC. The authors find that these data can be well described within a conventional collinear factorization framework. The analysis also includes inclusive pion production, which has not been traditionally included in global analyses, and finds that the shape of the nuclear modification factor $R_{\text{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. The uncertainties in the nPDFs are quantified using the Hessian method, which maps out the sensitivity of $\chi^2$ to variations in the fitting parameters to orthogonal error sets. The obtained NLO and LO nPDFs, along with their corresponding error sets, are collected in the new release called EPS09. These results are expected to be useful in precision analyses of QCD matter at the LHC and RHIC. The paper discusses the framework and analysis method, including the definition of nPDFs, experimental input, and cross-sections. It also presents the results and comparisons with data, including DIS, Drell-Yan, and inclusive pion production data. The authors compare their results with earlier global analyses and demonstrate the importance of the pion d+Au data in constraining the large-$x$ gluons. Finally, the paper applies the obtained EPS09 parametrization to a cross-section that was not included in the fit, specifically inclusive negative hadron production at forward rapidities in p+p and d+Au collisions, to illustrate how the parametrization should be applied in practice. The results show good agreement between the theory and the data, with the exception of the d+Au cross-section at $\eta = 3.2$, which lies below the pQCD prediction.The paper presents a global analysis of nuclear parton distribution functions (nPDFs) at next-to-leading order (NLO) accuracy in perturbative QCD (pQCD). The analysis includes three types of experimental data: deep inelastic lepton-nucleus scattering (DIS), Drell-Yan dilepton production in p+A collisions, and inclusive pion production in d+Au and p+p collisions at RHIC. The authors find that these data can be well described within a conventional collinear factorization framework. The analysis also includes inclusive pion production, which has not been traditionally included in global analyses, and finds that the shape of the nuclear modification factor $R_{\text{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. The uncertainties in the nPDFs are quantified using the Hessian method, which maps out the sensitivity of $\chi^2$ to variations in the fitting parameters to orthogonal error sets. The obtained NLO and LO nPDFs, along with their corresponding error sets, are collected in the new release called EPS09. These results are expected to be useful in precision analyses of QCD matter at the LHC and RHIC. The paper discusses the framework and analysis method, including the definition of nPDFs, experimental input, and cross-sections. It also presents the results and comparisons with data, including DIS, Drell-Yan, and inclusive pion production data. The authors compare their results with earlier global analyses and demonstrate the importance of the pion d+Au data in constraining the large-$x$ gluons. Finally, the paper applies the obtained EPS09 parametrization to a cross-section that was not included in the fit, specifically inclusive negative hadron production at forward rapidities in p+p and d+Au collisions, to illustrate how the parametrization should be applied in practice. The results show good agreement between the theory and the data, with the exception of the d+Au cross-section at $\eta = 3.2$, which lies below the pQCD prediction.