your Lagrangian is almost correct. But you also need to use a mass term that is conserved, which won't be the case of the mass term with your Lagrangian.
If you use: $$ \mathcal{L}_m = \sum_p m_p \gamma_p^{-1}(-g)^{-1/2}\delta^{(3)}(x^j-x^j_p(\tau_p)),$$ where $\gamma_p=dx^0/cd\tau_p$ the Lorentz factor, $u^\mu_p=dx^\mu_p/cd\tau$ the 4-velocity of the particle (such that $u^\mu u_\mu=-1$) and $x^j_p(\tau_p)$ the trajectory of the p-th particle, then you have all what you need. Indeed, you can directly compute that it reduces to $$S=-\sum_p m_p \int d \tau_p,$$ where $m_p$ is conserved.
Note that this is intimately related to the so-called conserved density $\rho^*$ that is often used in practical applications of general relativity (i.e.g. in cellestial mechanics). $\rho^*=\sqrt{-g} \gamma_p \rho$, where $\rho$ is the density appearing in the stress-energy tensor. Then, one has the "Newtonian" conservation of the so-called conserved density (i.e. $\partial_0 \rho^*+\partial_i(\rho^* v^i)=0$, with $v^i = u^i/\gamma_p$). (You can derive it from the conservation equation $\nabla_\sigma(\rho u^\sigma)=0$).