Lax and Phillips (J. Funct. Anal. vol 46 , 280--350 (1982)) 1982), 280--350) showed, for any crystallographic group $\Gamma$ in the Euclidean plane, that $$ N(r;x,x_0)= \frac{\pi r^2}{|F|} + O(s^{2/3O(r^{2/3} (\log s)^{1/2})r)^{1/2}), $$ as $r\to+\infty$, where $|F|$ denotes the volume of the fundamental domain of $\Gamma$, $x,x_0\in\mathbb R^2$ and $N(r;x,x_0)$ is the number of elements $\gamma\in\Gamma$ such that $$ |x-\gamma (x_0)|\leq r. $$
Later, Levitan (Russian Math. Surveys vol 42:3 (1987), 13--42) improved the error term to $O(r^{2/3})$.
In your particular case, $\Gamma$ is the subgroup of the isometries of the plane generated by the translations for $(1,0)$ and $(1/2,\sqrt{3}/2)$, thus $|F|=\sqrt{3}/2$.
Both papers works in higher dimensions and in hyperbolic spaces.

