  Using action-based distribution function for the dynamical model of the Milky
Way we have estimated its total mass and its density profile. Constraints are
coming from the globular cluster proper motions from Gaia EDR3, from the
rotation curve based on Gaia DR2 data, and from the vertical force data. We use
Bayesian Markov chain Monte Carlo method to explore the parameters, for which
the globular cluster distribution function and the Galactic potential are fully
constrained. Numerical simulations are used to study the uncertainties on the
potential constraint if considering a possible massive Large Magellanic Could
(LMC). We found that a massive LMC (1.5$\times10^{11}$ M$_{\odot}$) will affect
the MW mass measurement at large radius, which includes both the Milky Way and
the LMC. We also use the FIRE2 Latte cosmological hydrodynamic simulations to
make mock data set from a Milky-Way like galaxy that includes many unrelaxed
substructures. We test the effect of these unrelaxed substructures on the final
results, and found that the measured rotation curve fluctuated around input
value within 5 percent. By keeping a large freedom in choosing a priori mass
profile for both baryonic and dark matter leads a total mass of the MW that
ranges from $5.36_{-0.68}^{+0.81}\times10^{11}$ M$_{\odot}$ to
$7.84_{-1.97}^{+3.08} \times 10^{11}$ M$\odot$. This includes the contribution
of a putative massive LMC and significantly narrows the MW total mass range
published earlier. Such total mass leads to dark matter density at solar
position of $0.34_{-0.02}^{+0.02}$ GeV cm$^{-3}$.
