I will show the infimum is $0$ and give a bound for the supremum. Since $f$ is clearly positive, its supremum is positive, so $f$ is not constant. I don't know if there is any way to estimate the exact value of the supremum other than numerics.
Note that $f(g)$ is clearly right invariant by $SL_2(\mathbb Z)$ and left invariant by $SO_2(\mathbb R)$, i.e. it defines a function on the modular curve $X(1)$. My upper bound will show that $f(g)$ decreases to $0$ as $g$ approaches the cusp.
By the left $SO(2)$-invariance, we may assume $g = \begin{pmatrix} \lambda^{-1} & \alpha \\ 0 & \lambda \end{pmatrix}$ for $\lambda$ a positive real number. Then for $\gamma = \begin{pmatrix} a & b\\ c& d\end{pmatrix} \in SL_2(\mathbb Z)$ we have $$g\gamma = \begin{pmatrix} \lambda^{-1} a + \alpha c & \lambda^{-1} b + \alpha d \\ \lambda c & \lambda d \end{pmatrix}$$ so $$\|g \gamma\|^4= \left( \lambda^2 c^2 +\lambda^2 d^2 + (\lambda^{-1} a+ \alpha c)^2 + (\lambda^{-1} b + \alpha d)^2 \right)^2$$
Let's consider what happens to the most complicated term $(\lambda^{-1} a+ \alpha c)^2 + (\lambda^{-1} b + \alpha d)^2 $ if we replace $\begin{pmatrix} a & b\\ c& d\end{pmatrix}$ by $\begin{pmatrix} a+r c & b + rd\\ c& d\end{pmatrix}$ for $r \in \mathbb Z$. This gives
$$(\lambda^{-1} a + \lambda^{-1} rc + \alpha c)^2 + (\lambda^{-1} b + \lambda^{-1} rd + \alpha d)^2 $$
If we let $v_1 = \begin{pmatrix} \lambda^{-1} c \\ \lambda^{-1} d \end{pmatrix}$ and $v_2 = \begin{pmatrix} \lambda^{-1} a+ \alpha c \\ \lambda^{-1} b + \alpha d \end{pmatrix} $ then this is
$$ \| rv_1 + v_2\|^2 = r^2 \|v_1\|^2 + 2 r (v_1 \cdot v_2) + \|v_2\|^2 = \left(r + \frac{ v_1 \cdot v_2}{ \|v_1\|^2} \right)^2 \| v_1\|^2 + \| v_2\|^2 - \frac{ (v_1 \cdot v_2)^2}{ \|v_1\|^2} \geq \left(r + \frac{ v_1 \cdot v_2}{ \|v_1\|^2} \right)^2 \| v_1\|^2 =\left(r + \frac{ v_1 \cdot v_2}{ \|v_1\|^2} \right)^2 \lambda^{-2} (c^2+d^2) $$
Thus the sum over $r$ of $\frac{1}{ \| g\gamma\|^4}$ is bounded by
$$ \sum_{r\in \mathbb Z} \frac{1}{ \left( \lambda^2 c^2 +\lambda^2 d^2 + \left(r + \frac{ v_1 \cdot v_2}{ \|v_1\|^2} \right)^2 \lambda^{-2} (c^2+d^2) \right)^2 } $$
We can reparameterize the sum over $r$ so that $k=0$ is the closest value of $r$ to $-\frac{ v_1 \cdot v_2}{ \|v_1\|^2}$, $k=1$ is the second-closest, and so on. Then it's easy to see that $ \left| r + \frac{ v_1 \cdot v_2}{ \|v_1\|^2} \right| \geq \frac{k}{2}$ which means the sum is bounded by
$$ \sum_{k=0}^{\infty} \frac{1}{ \left( \lambda^2 c^2 +\lambda^2 d^2 + \frac{k^2}{4} \lambda^{-2} (c^2+d^2) \right)^2 } $$ $$= \sum_{k=0}^{\lceil 2\lambda^2 \rceil } \frac{1}{ \left( \lambda^2 c^2 +\lambda^2 d^2 + \frac{k^2}{4} \lambda^{-2} (c^2+d^2) \right)^2 } + \sum_{k= \lceil 2\lambda^2 \rceil+1 } ^{\infty} \frac{1}{ \left( \lambda^2 c^2 +\lambda^2 d^2 + \frac{k^2}{4} \lambda^{-2} (c^2+d^2) \right)^2 } $$$$\leq \sum_{k=0}^{\lceil 2\lambda^2 \rceil } \frac{1}{ \left( \lambda^2 c^2 +\lambda^2 d^2 \right)^2 } + \sum_{k= \lceil 2\lambda^2 \rceil+1 } ^{\infty} \frac{1}{ \left( \frac{k^2}{4} \lambda^{-2} (c^2+d^2) \right)^2 } $$
$$\leq \frac{ 2 \lambda^2+2}{ \left( \lambda^2 c^2 +\lambda^2 d^2 \right)^2 } + \int_{ \lceil 2\lambda^2 \rceil}^{\infty} \frac{1}{ \left( \frac{k^2}{4} \lambda^{-2} (c^2+d^2) \right)^2 } dk $$ $$ \leq \frac{ 2 \lambda^2+2}{ \left( \lambda^2 c^2 +\lambda^2 d^2 \right)^2 } + \int_{ 2\lambda^2 }^{\infty} \frac{1}{ \left( \frac{k^2}{4} \lambda^{-2} (c^2+d^2) \right)^2 } dk $$ $$ = \frac{ 2 \lambda^2+2}{ \left( \lambda^2 c^2 +\lambda^2 d^2 \right)^2 } + \frac{1}{ 3 (2\lambda^2)^3 \left( \frac{1}{4} \lambda^{-2} (c^2+d^2) \right)^2 } $$ $$ = \frac{ \frac{8}{3} \lambda^{-2}+2 \lambda^{-4} }{ \left( c^2 + d^2 \right)^2 } $$
Summing over $c,d$ this gives an upper bound of $$ \left( \frac{8}{3} \lambda^{-2}+2 \lambda^{-4} \right) \cdot \sum_{ \substack{ c,d \in \mathbb Z^2 \\ \gcd(c,d) =1 }} \frac{1}{ (c^2+d^2)^2 } $$
which clearly goes to $0$ as $\lambda \to \infty$, showing the infimum is zero, and showing the supremum is bounded by $$ \frac{64}{9} \cdot \sum_{ \substack{ c,d \in \mathbb Z^2 \\ \gcd(c,d) =1 }} \frac{1}{ (c^2+d^2)^2 } $$
as we can always take $\lambda \geq \frac{\sqrt{3}}{2}$ using the standard fundamental domain.
