In Euclidean $d=2k$ dimensions, we can take $\gamma^1, \gamma^2,\ldots, \gamma^{2k}$ to be hermitian. We can also arrange the usual $a_n, a^\dagger_n$ construction so that the $\gamma^i$ are symmetric for odd $i$ and antisymmetric for even $i$.
Assuming this to be the case, we define
$$
C_1 = \prod_{i{\rm\, odd}} \gamma^i, \quad C_2= \prod_{i{\rm \,even}} \gamma^i,
$$
and use them to construct ${\mathcal C}$ and ${\mathcal T}$ such that
$$
{\mathcal C}\gamma^i {\mathcal C}^{-1} = - (\gamma^i)^T,
$$
$$
{\mathcal T}\gamma^i {\mathcal T}^{-1} = + (\gamma^i)^T.\nonumber
$$
We find that for
{k=0, mod 4}:
${\mathcal C}=C_1$ symmetric, ${\mathcal T}=C_2$ symmetric. Both commute with $\gamma^{2k+1}$.
{k=1, mod 4}:
${\mathcal C}=C_2$ antisymmetric, ${\mathcal T}=C_1$ symmetric. Both anticommute with $\gamma^{2k+1}$.
{k=2, mod 4}:
${\mathcal C}=C_1$ antisymmetric, ${\mathcal T}=C_2$ antisymmetric. Both commute with $\gamma^{2k+1}$.
{k=3, mod 4}:
${\mathcal C}=C_2$ symmetric, ${\mathcal T}=C_1$ antisymmetric. Both anticommute with $\gamma^{2k+1}$.
Under a change of basis $\gamma^\mu \to A\gamma^\mu A^{-1}$ for the gamma matrices, ${\mathcal C}$ and ${\mathcal T}$ will no longer given by the explicit product expressions $C_1$ and $C_2$, but instead they transform as
$$
{\mathcal C} \to A^T {\mathcal C} A,\quad {\mathcal T} \to A^T {\mathcal T} A.
$$
The symmetry or antisymmetry of ${\mathcal C}$, ${\mathcal T}$ is unchanged, and is is thus a basis-independent property. Another way to think of this is by making use of the transpose of the the defining transformations. We then see that ${\mathcal C}^{-1}{\mathcal C}^T$ and ${\mathcal T}^{-1}{\mathcal T}^T$ commute with all $\gamma^\mu$ and so, by Schur's lemma, both ${\mathcal C}, {\mathcal T}$ are proportional to their transpose. Then, transposing again,
$$
C= \lambda C^T \Rightarrow C= \lambda^2 C
$$
so $\lambda=\pm 1$ in any representation. Since the Euclidean $\gamma^\mu$ are Hermitian, a similar argument shows that ${\mathcal C}^\dagger {\mathcal C}$ is proportional to the identity. As ${\mathcal C}^\dagger {\mathcal C}$ is a positive operator, the factor of proportionality is positive --- so ${\mathcal C}, {\mathcal T}$ can be scaled by a real factor so as to be unitary in all representations.
In odd dimensions we have to ask ourselves whether these objects commute or anticommute with the extra gamma matrix. We can construct (but not easy to display in MathJax) a table that displays the outcome of this examination. It shows whether the matrices exist, and whether they are symmetric (S), or antisymmetric (A). The table repeats mod 8.
When ${\mathcal T}$ is symmetric we can find a basis in which ${\mathcal T}= {\mathbb I}$ and all the Euclidean gamma matrices are real. When ${\mathcal C}$ is symmetric we can find a basis in which all the gamma matrices are pure imaginary.
Charge conjugation: We define the charge conjugate fields by
$$
\psi^c= {\mathcal C}^{-1}{\bar \psi}^T= {\mathcal C}^{-1}(\gamma^0)^T \psi^*,
$$
$$
\bar \psi^c = -\psi^T {\mathcal C}.\nonumber
$$
The transformation of $\bar \psi$ (which in Euclidean space is independent of $\psi$) is chosen to be consistent with the Minkowski version. To obtain this Minkowski signature version we first note that we can use the same ${\mathcal T}$ and ${\mathcal C}$ matrices as in Euclidean space --- the inclusion of factors of $i$ in the transposing $\gamma^\mu$ does not change the formula. Consider, for example, the mostly minus Minkowski metric $(+,-,-,\ldots)$ in which $\gamma^0$ is Hermitian and obeys $(\gamma^0)^2=1$. Then
$$
\psi^c= {\mathcal C}^{-1} (\gamma^0)^T \psi^*\, \Rightarrow (\psi^c)^\dagger = \psi ^T(\gamma^0)^T {\mathcal C},
$$
as ${\mathcal C}$ remains unitary in Minkowski space. Then we have
$$
\overline{(\psi^c)}\equiv (\bar \psi)^c = (\psi^c)^\dagger \gamma^0
$$
$$
= \psi^T (\gamma^0)^T {\mathcal C} \gamma^0 $$
$$
=- \psi^T \,{\mathcal C} \gamma^0 {\mathcal C}^{-1}{\mathcal C} \gamma^0$$
$$
= -\psi^T \,{\mathcal C}.\nonumber
$$
From these results, and with anticommuting $\psi$'s, we find that
$$
\bar \psi^c \gamma^\mu\psi^c =[ - \psi^T {\mathcal C}]\gamma^\mu [{\mathcal C}^{-1} (\gamma^0)^T\psi^*]= \psi^T (\gamma^\mu)^T(\gamma^0)^T \psi^* = - \psi^*\gamma^0 \gamma^\mu \psi =- \bar \psi \gamma^\mu\psi.
$$
Similarly, we have
$$
\bar \psi^c \psi^c =- \psi^T \bar \psi^T = \bar\psi\psi.
$$
The spin-current density transforms as
$$
\bar \psi^c \gamma^0 [\gamma^i,\gamma^j]\psi^c= -\bar\psi \gamma^0 [\gamma^j, \gamma^i]\psi = \bar \psi \gamma^0 [\gamma^i,\gamma^j]\psi.
$$
So, as a result of multiplying six (!) minus signs, the spin is left unchanged.
Using the anticommuting property of the Fermi fields, we now find that, in both Euclidean and Minkowski signatures,
$$
S= \int d^n x\, \bar \psi [\gamma^\mu (\partial_\mu+A)+m]\psi = \int d^n x\, {\bar\psi}^c [\gamma^\mu (\partial_\mu-A^T)+m]\psi^c.
$$
The $-A^T$ are the Lie algebra representation-valued fields in the the conjugate representation to that of $A$, and so $\psi_c$ has the opposite ``charge'' to $\psi$.
Minkowski-signature Marjorana Fermions: We have defined
$$
\psi^c= {\mathcal C}^{-1} (\gamma^0)^T \psi^*
$$
so, with ${\mathcal C}^T=\lambda {\mathcal C}$ we find
$$
(\psi^c)^c= {\mathcal C}^{-1} (\gamma^0)^T ({\mathcal C}^{-1} (\gamma^0)^T\psi^*)^*$$
$$
= {\mathcal C}^{-1}(\gamma^0)^T {\mathcal C}^T\gamma^0 \psi$$
$$
= \lambda \,{\mathcal C}^{-1}(\gamma^0)^T {\mathcal C}\gamma^0 \psi$$
$$
= -\lambda \gamma^0 \gamma^0\psi$$
$$
= -\lambda \psi.$$
We can therefore consistently impose the Majorana condition that $\psi^c=\psi$ only if ${\mathcal C}$ is antisymmetric: i.e. in 2, 3, 4 (mod 8) dimensions. An alternative way to view this is to regard the map
$\psi\to \psi^c$ as an antilinear map $J: {\rm Pin} \to {\rm Pin} $ where ${\rm Pin}$ is the gamma-matrix representation space. If $J^2=id$, this is real structure on the complex ${\rm Pin}$ space.
We can also define another conjugate
$$
\psi^t= {\mathcal T}^{-1}{\bar \psi}^T
$$
$$
\bar \psi^t = \psi^T {\mathcal T}.\nonumber
$$
This reverses the current, again leaves the spin unchanged, but flips the sign of $\bar\psi \psi$. If we demand that $(\psi^t)^t=\psi$, almost identical algebra shows that this defines a real structure only when ${\mathcal T}$ is symmetric hence in d= 8, 9 (mod 8). (The case d=10 (mod 8) is Majorana Weyl.) Some people regard gamma-matrix representations with this real structure also as being Marjorana. I think that these are what Jose} Figueroa-O'Farril calls pseudo-Majorana. Because this``conjugation'' flips $\bar\psi\psi$, these pseudo-Majorana fermions are necessarily massless. It is unfortunate that in the mostly-minus West-cast metric the gamma matrices of a pseudo-Majorana representation can be chosen to be real, while in a Majorana they can only be chosen pure imaginary. It is the other way around in the mostly-plus East-coast metric.
Euclidean Marjorana Fermions: Suppose that we have an eigenfunction of the (skew-adjoint) Dirac operator $D=\gamma^\mu \partial_\mu$ such that
$$
{D}u_n=i\lambda_n u_n,
$$
then
$$
{D}^* u^*=-i\lambda u^*_n.
$$
Assuming that any gauge fields are in real representations, this can be written as
$$
{\mathcal C} {D} {\mathcal C}^{-1} u^*_n = i\lambda_n u^*_n,
$$
or
$$
{D} {\mathcal C}^{-1} u^*_n = i\lambda_n {\mathcal C}^{-1} u^*_n.
$$
Thus $u_n$ and ${\mathcal C}^{-1} u_n^*$ are both eigenfunctions of ${D}$ with the same eigenvalue. They will be linearly independent when ${\mathcal C}$ is antisymmetric --- something that happens in $2,3,4$ (mod 8) Euclidean dimensions. These are precisely the dimensions in which {\it Minkowski space/} Majorana spinors can occur.
This suggests that we take the
Euclidean Majorana-Dirac action to be
$$
S[\psi]= \frac 12 \int d^nx \,\psi^T {\mathcal C}( {D}+m)\psi
$$
and setting
$$
\psi (x) = \sum_n[ \xi_n u_{n}(x)+ \eta_n( {\mathcal C}^{-1}u_{n}^*(x)],
$$
$$
\psi^T(x) = \sum_n [\xi_n u^T_n(x) - \eta_n (u_n^\dagger(x) {\mathcal C}^{-1})],\nonumber
$$
so that
$$
S= \frac 12 \int d^nx \sum_{i,j} ( \xi_i u^T_i - \eta_i u_i^\dagger {\mathcal C}^{-1}) {\mathcal C}(i\lambda_j+m)( \xi_j u_j+ \eta_j {\mathcal C}^{-1}u_j^*)
$$
$$
= \frac 12 \int d^nx \sum_{i, j}\left\{ \xi_i\eta_j (u^T_i u_j^*) + \xi_j \eta_i (u_i^\dagger u_j)\right\}(i\lambda_j+m)
$$
$$
= \sum_i \xi_i\eta_i (i\lambda_i+m).\nonumber
$$
Grassmann integration uses only one copy of the doubly degenerate eigenvalue, and so gives the square-root of the full Dirac determinant.
Best Answer
The transformation rule for $\hat{C}$ or that of $\hat{C}\gamma_{0}^{T}$ is different from the usual one.
The whole charge conjugation operation is given by $VK$, where $V\equiv\hat{C}\gamma_{0}^{T}$ is unitary and $K$ is complex conjugation. (i.e., $VK$ is antiunitary.) Then, under a basis transformation, $$ VK \ \rightarrow \ U^{\dagger}VK U = U^{\dagger}VU^{\ast} K. $$ Hence $V$ transforms as $V \ \rightarrow \ U^{\dagger}VU^{\ast}$. My guess is that this is the piece that was missing in your derivation.