$\require{cancel}$I) OP is considering Dirac fermions in a curved spacetime. OP's action has various shortcomings. The correct action reads$^1$
$$ S~=~\int\!d^nx~ {\cal L}, \qquad {\cal L} ~=~e L, \qquad L~=~T-V,\qquad e~:=~\det(e^a{}_{\mu})~=~\sqrt{|g|}, $$
$$ T~=~\frac{i}{2} \bar{\psi}
\stackrel{\leftrightarrow}{\cancel{\nabla}} \psi, \qquad
V~=~ \alpha j^a \eta_{ab} j^b,
\qquad j^a~:=~ \bar{\psi} \gamma^a\psi,\qquad \bar{\psi}~:=~\psi^{\dagger}\gamma^0,$$
$$ \bar{\psi}\stackrel{\leftrightarrow}{\cancel{\nabla}}\psi
~:=~ \bar{\psi}\stackrel{\leftrightarrow}{\cancel{\partial}}\psi
+\frac{1}{2} \omega_{c, ab}~\gamma^{cab}\psi
~=~\bar{\psi}\left[\gamma^c\stackrel{\rightarrow}{\nabla_c}
-\stackrel{\leftarrow}{\nabla_c}\gamma^c\right]\psi, $$
$$\stackrel{\rightarrow}{\nabla_c}\psi
~:=~ \stackrel{\rightarrow}{\partial_c}\psi +\frac{1}{4} \omega_{c, ab}~\gamma^{ab}\psi, \qquad
\bar{\psi}\stackrel{\leftarrow}{\nabla_c}
~:=~ \bar{\psi}\stackrel{\leftarrow}{\partial_c} -\frac{1}{4} \bar{\psi}~\gamma^{ab}\omega_{c, ab},$$
$$\stackrel{\leftrightarrow}{\cancel{\partial}} ~:=~ \gamma^c\stackrel{\rightarrow}{\partial_c} - \stackrel{\leftarrow}{\partial_c}\gamma^c, \qquad
\stackrel{\rightarrow}{\partial_c}~:=~E^{\mu}{}_c \stackrel{\rightarrow}{\partial_{\mu}}, \qquad
\stackrel{\leftarrow}{\partial_c}~:=~ \stackrel{\leftarrow}{\partial_{\mu}}E^{\mu}{}_c,$$
$$
\partial_{\mu}~:=~\frac{\partial}{\partial x^{\mu}}, \qquad \gamma^{ab}~:=~\frac{1}{2}[\gamma^a,\gamma^b], \qquad
\gamma^{abc}~:=~\frac{1}{2}\{\gamma^a,\gamma^{bc}\}_+.
\tag{1} $$
II) The main point is that in order to write down a covariant kinetic term for a Dirac fermion in curved spacetime, we should use a covariant derivative $\nabla_{\mu}\psi$ of a spinor $\psi$, and hence we need a spin connection $\omega_{\mu}{}^a{}_b$. In turn, we need a vielbein
$$g_{\mu\nu}~=~e^a{}_{\mu} ~\eta_{ab}~e^b{}_{\nu}, \qquad
e^a{}_{\mu}~ E^{\mu}{}_b~=~\delta^a_b, \qquad
E^{\mu}{}_a~e^a{}_{\nu}~=~\delta^{\mu}_{\nu}, \tag{2} $$
which (we for simplicity will assume) is covariantly conserved
$$0~=~(\nabla_{\mu}e)^{a}{}_{\nu}~=~\partial_{\mu}e^{a}{}_{\nu} +\omega_{\mu}{}^a{}_b ~e^b{}_{\nu}- e^a{}_{\lambda}~\Gamma_{\mu\nu}^{\lambda}.\tag{3} $$
Hence the spin connection is complete determined
$$ 2\omega_{\mu, ab}~=~2\left(-\partial_{\mu}e_{a\nu}
+e_{a\lambda}~\Gamma_{\mu\nu}^{\lambda}\right)E^{\nu}{}_b
~=~-\left(\partial_{\mu}e_{a\nu}
+\partial_a g_{\mu\nu}\right)E^{\nu}{}_b -(a\leftrightarrow b)$$
$$ ~=~-\partial_{\mu}e_{a\nu}~E^{\nu}{}_b-\partial_a e_{b\mu} + g_{\mu\nu}~\partial_a E^{\nu}{}_b -(a\leftrightarrow b),\tag{4} $$
and
$$ 2\omega_{c, ab}~:=~2E^{\mu}{}_c~\omega_{\mu, ab}
~=~-f_{cab}-f_{abc}-f_{acb}-(a\leftrightarrow b), \tag{5}$$
where we defined
$$f_{abc}~:=~\partial_a e_{b\nu}~E^{\nu}{}_c . \tag{6}$$
III) The kinetic term becomes
$$ T~=~\frac{i}{2}\bar{\psi}\stackrel{\leftrightarrow}{\cancel{\nabla}}\psi
~=~\frac{i}{2}\bar{\psi}\stackrel{\leftrightarrow}{\cancel{\partial}}\psi -\frac{i}{4}\bar{\psi}~f_{abc}~\gamma^{cab}~\psi $$
$$ ~=~\frac{i}{2}\bar{\psi}\left[\gamma^c~E^{\mu}{}_c\stackrel{\rightarrow}{\partial_{\mu}} -\stackrel{\leftarrow}{\partial_{\mu}}E^{\mu}{}_c~\gamma^c
\right]\psi
-\frac{i}{4}\bar{\psi}~E^{\mu}{}_a~\partial_{\mu} e_{b\nu}~E^{\nu}{}_c
~\gamma^{cab}~\psi. \tag{7}$$
IV) The natural generalization of the Hilbert SEM tensor
$$T^{\mu\nu}~=~- \frac{2}{\sqrt{|g|}}\frac{\delta S}{\delta g_{\mu\nu}}, \qquad
T_{\mu\nu}~=~\frac{2}{\sqrt{|g|}}\frac{\delta S}{\delta g^{\mu\nu}},\qquad(\leftarrow\text{Not applicable!})\qquad
\tag{8} $$
to fermions is given by the formula
$$T^{\mu\nu}~=~-\frac{E^{\mu c}}{2e}\frac{\delta S}{\delta e^c{}_{\nu}}+(\mu\leftrightarrow \nu), \qquad T_{\mu\nu}~=~\frac{e_{c\mu}}{2e}\frac{\delta S}{\delta E^{\nu}{}_c}+(\mu\leftrightarrow \nu).\tag{9}$$
Formula (9) reduces to the standard Hilbert SEM tensor (8) if the action only depends on the vielbein through the metric (2). However formula (9) is more general and is necessary in the case of fermions in curved spacetime.
V) The Hilbert SEM tensor with flat indices then becomes
$$T_{cd}~:=~ E^{\mu}{}_c~ T_{\mu\nu} ~E^{\nu}{}_d~\stackrel{(9)}{=}~-\frac{e_{c\nu}}{2e}\frac{\delta S}{\delta e^d{}_{\nu}} +(c\leftrightarrow d)
~=~\frac{E^{\nu}{}_c}{2e}\frac{\delta S}{\delta E^{\nu d}} +(c\leftrightarrow d)$$
$$~\stackrel{(7)}{=}~\frac{i}{4}\bar{\psi}\left[\gamma_c\stackrel{\rightarrow}{\partial_d} -\stackrel{\leftarrow}{\partial_c}\gamma_d
+\frac{1}{2} (f_{cba}-f_{abc}-f_{acb})~\gamma_d{}^{ab} \right]\psi
-\frac{1}{2}\eta_{cd}L+(c\leftrightarrow d)$$
$$~\stackrel{(5)}{=}~\frac{i}{4}\bar{\psi}\left[\gamma_c\stackrel{\rightarrow}{\partial_d} -\stackrel{\leftarrow}{\partial_c}\gamma_d
+\frac{1}{2} \omega_{c,ab}~\gamma_d{}^{ab} \right]\psi
-\frac{1}{2}\eta_{cd}L+(c\leftrightarrow d)$$
$$~\stackrel{(1)}{=}~\frac{i}{4}\bar{\psi}\left[\gamma_c\stackrel{\rightarrow}{\nabla_d} -\stackrel{\leftarrow}{\nabla_c}\gamma_d\right]\psi
-\frac{1}{2}\eta_{cd}L+(c\leftrightarrow d).
\tag{10} $$
Eq. (10) is the formula for the (generalized) Hilbert SEM tensor of a Dirac fermion in curved spacetime. This is the appropriate matter source term in the EFE, cf. OP's title question (v3). For further details, see also my Phys.SE answers here and here.
--
$^1$ One may show that the Lagrangian density (1) is real using
$$ (\gamma^a)^{\dagger}~=~ \gamma^0\gamma^a\gamma^0,\qquad (\gamma^0)^2~=~{\bf 1}.\tag{11} $$
Conventions: In this answer, we will use $(+,-,-,-)$ Minkowski sign convention, and Clifford algebra
$$\{\gamma^a,\gamma^b\}_+~=~2\eta^{ab}{\bf 1}.\tag{12}$$
Greek indices $\mu,\nu,\lambda, \ldots,$ are so-called curved indices, while Roman indices $a,b,c, \ldots,$ are so-called flat indices.