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 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 \beta of the companion star , where \beta = 1 indicates a completely filled Roche lobe . In previous work a range of \beta between 0.9 and 1.0 was usually adopted . We use optical ellipsoidal light-curve data to constrain \beta . We find neutron star masses of 1.77 \pm 0.08 ~ { } M _ { \odot } for Vela X-1 , 0.87 \pm 0.07 ~ { } M _ { \odot } for 4U 1538-52 ( eccentric orbit ) , 1.00 \pm 0.10 ~ { } M _ { \odot } for 4U 1538-52 ( circular orbit ) , 1.04 \pm 0.09 ~ { } M _ { \odot } for SMC X-1 , 1.29 \pm 0.05 ~ { } M _ { \odot } for LMC X-4 , 1.49 \pm 0.08 ~ { } M _ { \odot } for Cen X-3 , and 1.07 \pm 0.36 ~ { } M _ { \odot } 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 .