Without any more information about M,f or g, I'm pretty sure the best you can expect is that (f-g) is in the kernel of Hess. You could get an estimate on the size of (f-g) with a Poincare inequality or something (?). If you have boundary conditions then you could uniquely solve Hess(f-g)=0 for (f-g).
I'm not sure this is what you're looking for though (?)
If the manifold is Rn, and the functions are smooth, then Hess f = Hess g means that Hess(f-g) = 0 so grad(f-g) = constant. That is, f & g differ by a linear function. I'm not quite sure what this translates to for a general smooth Riemannian manifold.
To answer @David Steward's question, on a Riemannian manifold
Hess(f - g) = \nabla (d (f-g) ) = 0
means \theta = d(f-g) is a parallel 1-form, of which there typically are none other than the 0 form. Hence in the generic case d(f-g) = 0, and f and g differ by a locally constant function (i.e. a constant for every connected component).
However since \theta is also exact it is even more restricted than just being parallel.
The Hessian is the covariant derivative of the gradient. Then the gradient, $gr(f-g)$, of $(f-g)$ is a parallel vector field. If it is zero, then f and h differ by a constant. If not, by de Rhan decomposition theorem $M$ splits (locally or globally if $M$ is complete and simply connected) a line, i.e. $M= N\times \mathbb R$ (locally or globally). In this way $grad (f-g)$ is a constant field in the $\mathbb R$-factor of $M$. So $f$ differs from $g$ in a function wich is a multiple of the projection to the line factor $\mathbb R$.
Thus, if $M= \mathbb R ^k\times M'$, where $\mathbb R ^k$ is the Euclidean de Rham factor of $M$, $f-g = T\circ \pi_0$, where $T:\mathbb R ^k \to \mathbb R$ is affine and
$\pi _0$ is the projection from $M$ to the Euclidean de Rham factor $\mathbb R ^k$.