We present an improved method for determining the mass of neutron stars in eclipsing X-ray pulsar binaries and apply the method to six systems, namely, Vela X-1, 4U 1538-52, SMC X-1, LMC X-4, Cen X-3, and Her X-1. In previous studies to determine neutron star mass, the X-ray eclipse duration has been approximated analytically by assuming that the companion star is spherical with an effective Roche lobe radius. We use a numerical code based on Roche geometry with various optimizers to analyze the published data for these systems, which we supplement with new spectroscopic and photometric data for 4U 1538-52. This allows us to model the eclipse duration more accurately and thus calculate an improved value for the neutron star mass. The derived neutron star mass also depends on the assumed Roche lobe filling factor Delta *b of the companion star, where Delta *b = 1 indicates a completely filled Roche lobe. In previous work a range of Delta *b between 0.9 and 1.0 was usually adopted. We use optical ellipsoidal light-curve data to constrain Delta *b. We find neutron star masses of 1.77 ? 0.08 M for Vela X-,0.87 ? 0.07 M for 4U 1538-52 (eccentric orbit), 1.00 ? 0.10 M for 4U 1538-52 (circular orbit), 1.04 ? 0.09 M for SMC X-1, 1.29 ? 0.05 M for LMC X-4, 1.49 ? 0.08 M for Cen X-3, and 1.07 ? 0.36 M for Her X-1. We discuss the limits of the approximations that were used to derive the earlier mass determinations, and we comment on the implications our new masses have for observationally refining the upper and lower bounds of the neutron star mass distribution.