A precise knowledge of nuclear parton distribution functions (nPDFs) is -- among other things -- important for the unambiguous interpretation of hard process data taken in pA and AA collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). The available fixed target data for deep inelastic scattering (DIS) and Drell-Yan (DY) lepton pair production mainly constrain the light quark distributions. It is hence crucial to include more and more collider data in global analyses of nPDFs in order to better pin down the different parton flavors, in particular the gluon distribution at small x. To help constrain the nuclear gluon PDF, we extend the nCTEQ15 analysis by including single inclusive hadron (SIH) production data from RHIC (PHENIX and STAR) and LHC (ALICE). In addition to the DIS, DY and SIH data sets, we will also include LHC W/Z production data. As the SIH calculation is dependent on hadronic fragmentation functions (FFs), we use a variety of FFs available in the literature to properly estimate this source of uncertainty. We study the impact of these data on the PDFs, and compare with both the nCTEQ15 and nCTEQ15WZ sets. The calculations are performed using a new implementation of the nCTEQ code (nCTEQ++) including a modified version of INCNLO which allows faster calculations using pre-computed grids. The extension of the nCTEQ15 analysis to include the SIH data represents an important step toward the next generation of PDFs.