One can take the continuum limit of your proof as $X \to \infty$, again using the prime number theorem, to obtain a proof that does not involve primes at all:
$$ \int f(t) \log \frac{1}{|t|}\ dt = \gamma - \sum_{\sigma = \pm 1} \int_0^\infty f(\sigma t) (\log t + \gamma)\ dt $$ $$ = \gamma - \lim_{\varepsilon \to 0} \sum_{\sigma = \pm 1} \int_0^\infty f(\sigma t) (\log \frac{t}{\varepsilon} + \gamma)\ dt + \log \frac{1}{\varepsilon}$$ $$ = \gamma - \lim_{\varepsilon \to 0} \sum_{\sigma = \pm 1} \int_0^\infty f(\sigma t) (\sum_{0 < s < t/\varepsilon} \frac{1}{s})\ dt + \log \frac{1}{\varepsilon}$$ $$ = \gamma - \lim_{\varepsilon \to 0} \sum_{\sigma = \pm 1} \sum_{s>0} \frac{1}{s} \int_{\varepsilon s}^\infty f(\sigma t)\ dt + \log \frac{1}{\varepsilon}$$ $$ = \gamma - \lim_{\varepsilon \to 0} \sum_{\sigma = \pm 1} \sum_{s>0} \int_{\varepsilon}^\infty f(\sigma s t)\ dt + \log \frac{1}{\varepsilon}$$ $$ = \gamma - \lim_{\varepsilon \to 0} \int_{\varepsilon}^\infty \sum_{s \in \mathbb{Z} \backslash \{0\}} f(s t)\ dt + \log \frac{1}{\varepsilon}$$ $$ = \gamma - A - \lim_{\varepsilon \to 0} \int_{\varepsilon}^1 \sum_{s \in \mathbb{Z} \backslash \{0\}} f(s t)\ dt + \log \frac{1}{\varepsilon}$$ $$ = \gamma - A - \lim_{\varepsilon \to 0} \int_{\varepsilon}^1 (\sum_{s \in \mathbb{Z} \backslash \{0\}} f(s t) - \frac{1}{t})\ dt $$ $$ = \gamma - A - \lim_{\varepsilon \to 0} \int_{\varepsilon}^1 (\frac{1}{t} \sum_{s \in \mathbb{Z} \backslash \{0\}} \hat f(s/t) - 1)\ dt $$ $$ = \gamma + 1 - A - B$$ where $$ A := \int_1^\infty \sum_{s \in \mathbb{Z} \backslash \{0\}} f(st)\ dt$$ $$ = \int_{|t| \geq 1} f(t) (\sum_{1 \leq s \leq |t|} \frac{1}{s})\ dt$$ and $$ B := \int_0^1 \sum_{s \in \mathbb{Z} \backslash \{0\}} \hat f(s/t)\ \frac{dt}{t}$$ $$ = \int_1^\infty \sum_{s \in \mathbb{Z} \backslash \{0\}} \hat f(st)\ \frac{dt}{t}$$ $$ = \int_{|t| \geq 1} \hat f(t) \frac{\lfloor |t| \rfloor}{t}\ dt.$$ [In the language of distributions, what this identity is saying I think is that the distributional Fourier transform of $\lfloor |t| \rfloor/t - 1$ is $\gamma - \log \frac{1}{|t|} 1_{|t| \leq 1} - \sum_{1 \leq s \leq |t|} \frac{1}{s}$.]
Since $A,B$ are clearly non-negative, this gives your inequality. This also shows that one is within $o(1)$ of equality if and only if one simultaneously has
$$ \int_{|t| \geq 1} f(t) (1 + \log |t|)\ dt = o(1)$$
and
$$ \int_{|t| \geq 1} \hat f(t)\ dt = o(1).$$
By the Hahn-Banach theorem, these estimates are incompatible with the hypotheses $f(0)=\hat f(0)=1$, $f(t), \hat f(t) \geq 0$ for $|t| \geq 1$ if and only if there exist non-negative measurable functions $a(t), b(t)$ supported on $|t| \geq 1$ with $\sup_t \frac{a(t)}{1+\log |t|}, \sup_t b(t) < \infty$ and numbers $\alpha,\beta$ not summing to zero, such that
$$ \alpha f(0) + \beta \hat f(0) = \int_{\mathbb R} f(t) a(t)\ dt + \int_{\mathbb R} \hat f(t) b(t)\ dt $$
for all Schwartz $f$, or equivalently that
$$ \alpha \delta + \beta = a + \check b$$
in the sense of tempered distributions, where $\delta$ is the Dirac delta. But the right-hand side is continuous at the origin, so $\alpha$ must vanish; the Fourier transform of the right-hand side has a continuous antiderivative at the origin, so $\beta$ must vanish, contradiction. This shows that one can make $A$ and $B$ simultaneously $o(1)$, so $1+\gamma$ is in fact optimal. (But the invocation of the Hahn-Banach theorem makes it difficult to explicitly construct $f$ that come close to equality!)
One can solve this equation as follows. By Lemma 3 of
Amrein, W.O.; Berthier, A.M., On support properties of Lsup(p)-functions and their Fourier transforms, J. Funct. Anal. 24, 258-267 (1977). ZBL0355.42015,
one can find, for any $R>0$, a non-zero function $f \in L^2({\bf R})$ such that $f$ and $\hat f$ both vanish on $[-R,R]$ (this is basically because the compact operator $1_{[-R,R]} {\mathcal F} 1_{[-R,R]}$ is a strict contraction on $L^2$, which in turn follows from the uncertainty principle that a function and its Fourier transform cannot be simultaneously compactly supported), in fact Proposition 6 gives an infinite-dimensional space of such functions. By convolving $f$ by a suitable approximation to the identity, and then multiplying by the Fourier transform of a suitable approximation to the identity (and shrinking $R$ slightly), one can make $f$ Schwartz.
When one takes a second antiderivative of $f$, one obtains a new Schwartz function $f_1$ which is equal to a linear function $a+\beta x$ on $[-R,R]$, while the Fourier transform still vanishes on $[-R,R]$. If $a=0$ (which can be achieved due to the infinite dimensional nature of the space of $f$), one can divide by $x$ and obtain a further Schwartz function $f_2$ that is equal to a constant $\beta$ on $[-R,R]$, while the Fourier transform is equal to a constant $\alpha$ on $[-R,R]$. This gives the identity $$ \alpha \delta + \beta = (\beta - f) + (\alpha - \hat f)^{\vee}$$ I think one can work a little harder to ensure that $\alpha,\beta$ can be arbitrary real numbers while simultaneously keeping $a=0$, and in particular can have non-zero sum (otherwise by Hahn-Banach there would be a way to express some nontrivial combination of polynomials on a halfline and Fourier transforms of polynomials on a halfline as functions supported on $[-R,R]$ plus a function with Fourier transform supported on $[-R,R]$, which should be easy to rule out by the argument in strikethrough). This gives a constraint of the desired form (taking $R=1$). So some improvement to $1+\gamma$ is in fact possible.