  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 lightcurve data to constrain beta. We find neutron star masses of
1.77 +/- 0.08 M_{sun} for Vela X-1, 0.87 +/- 0.07 M_{sun} for 4U 1538-52
(eccentric orbit), 1.00 +/- 0.10 M_{sun} for 4U 1538-52 (circular orbit), 1.04
+/- 0.09 M_{sun} for SMC X-1, 1.29 +/- 0.05 M_{sun} for LMC X-4, 1.49 +/- 0.08
M_{sun} for Cen X-3, and 1.07 +/- 0.36 M_{sun} 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.
