(* Version May 9, 2022
Introduced Jordan normal form to handle non-diagonalisable cases.
Canged normalisation of SW term operators.
Added If[GW,...] to treat matching for Wilson and GW quarks differently.

Version July 20, 2021
Fixed a mistake in massive mixing for the operator (Subscript[g, 0])²(ΨΨ)².
Added last two massive operators, which do not mix with any other operators at 1-loop but are required for a complete basis.
Normalised tr(M^n) by Nf to treat the quenched case properly.
*)

(* non diagonal basis: *)
basis={{"1/(Subscript[g, 0])² tr((Subscript[D, μ])(Subscript[F, μρ])(Subscript[D, ν])(Subscript[F, νρ]))","1/(Subscript[g, 0])²(∑_μ)tr((Subscript[D, μ])(Subscript[F, μρ])(Subscript[D, μ])(Subscript[F, μρ]))","(∑_μ)(Subscript[Ψγ, μ])(Subscript[D, μ])³Ψ","(Subscript[g, 0])²((Subscript[Ψγ, μ])Ψ)²","(Subscript[g, 0])²((Subscript[Ψγ, μ])(Subscript[γ, 5])Ψ)²","(Subscript[g, 0])²((Subscript[Ψγ, μ])(T^a)Ψ)²","(Subscript[g, 0])²((Subscript[Ψγ, μ])(Subscript[γ, 5])(T^a)Ψ)²","(Subscript[g, 0])²(ΨΨ)²","(Subscript[g, 0])²((Subscript[Ψγ, 5])Ψ)²","(Subscript[g, 0])²(Ψ(Subscript[σ, μν])Ψ)²","(Subscript[g, 0])²(((ΨT)^a)Ψ)²","(Subscript[g, 0])²((Subscript[Ψγ, 5])(T^a)Ψ)²","(Subscript[g, 0])²(Ψ(Subscript[σ, μν])(T^a)Ψ)²","i/4 (Subscript[ΨMσ, μν])(Subscript[F, μν])Ψ","i tr(M)/(4 Nf) Ψ(Subscript[σ, μν])(Subscript[F, μν])Ψ","1/(Subscript[g, 0])² tr(M²)/Nf tr((Subscript[F, μν])(Subscript[F, μν]))","1/(Subscript[g, 0])² tr(M)²/Nf² tr((Subscript[F, μν])(Subscript[F, μν]))","ΨM³Ψ","tr(M)/Nf ΨM²Ψ","tr(M²)/Nf ΨMΨ","tr(M)²/Nf² ΨMΨ","tr(M³)/Nf ΨΨ","tr(M²)tr(M)/Nf² ΨΨ","tr(M)³/Nf³ ΨΨ"}}//Transpose;

(* different subblocks of the mixing matrix (alp=alpha/(4πϵ), 2ϵ=4-D) *)

massless=(({{1+((7 alp Nc)/3)+((2 alp Nf)/3), 0, 0, 0, ((3 alp)/2)-((3 alp)/(2 (Nc^2))), ((2 alp)/(3 Nc))+((alp Nc)/6), ((6 alp)/Nc)-((3 alp Nc)/2), 0, 0, 0, 0, 0, 0}, {-((7 alp Nc)/15), 1+((21 alp Nc)/5)+((2 alp Nf)/3), ((11 alp)/(60 Nc))-((11 alp Nc)/60), 0, (alp/2)-(alp/(2 (Nc^2))), ((19 alp)/(80 Nc))+((47 alp Nc)/240), ((2 alp)/Nc)-((alp Nc)/2), 0, 0, 0, 0, 0, 0}, {((11 alp Nf)/60), -((11 alp Nf)/15), 1-((157 alp)/(60 Nc))+((157 alp Nc)/60), 0, (alp/8)-(alp/(8 (Nc^2))), -((13 alp)/(80 Nc))+((3 alp Nc)/80)-((7 alp Nf)/60), (alp/(2 Nc))-((alp Nc)/8), 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1+((11 alp Nc)/3)-((2 alp Nf)/3), 0, -((4 alp)/3), -6 alp, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1+((11 alp Nc)/3)-((2 alp Nf)/3), -((22 alp)/3), 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, -((3 alp)/2)+((3 alp)/(2 (Nc^2))), 1-((2 alp)/(3 Nc))+((13 alp Nc)/6)+((2 alp Nf)/3), -((6 alp)/Nc)+((3 alp Nc)/2), 0, 0, 0, 0, 0, 0}, {0, 0, 0, -((3 alp)/2)+((3 alp)/(2 (Nc^2))), 0, -((20 alp)/(3 Nc))+((3 alp Nc)/2), 1+((13 alp Nc)/6)-((2 alp Nf)/3), 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, ((2 alp)/3), 0, 1+((3 alp)/Nc)+((2 alp Nc)/3)-((2 alp Nf)/3), 0, 0, 0, 0, alp}, {0, 0, 0, 0, 0, -((2 alp)/3), 0, 0, 1+((3 alp)/Nc)+((2 alp Nc)/3)-((2 alp Nf)/3), 0, 0, 0, alp}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1-(alp/Nc)+((14 alp Nc)/3)-((2 alp Nf)/3), 24 alp, 24 alp, 0}, {0, 0, 0, 0, 0, (alp/(3 Nc)), 0, 0, 0, (alp/4)-(alp/(4 (Nc^2))), 1+((3 alp)/Nc)+((11 alp Nc)/3)-((2 alp Nf)/3), 0, (alp/Nc)-((alp Nc)/4)}, {0, 0, 0, 0, 0, -(alp/(3 Nc)), 0, 0, 0, (alp/4)-(alp/(4 (Nc^2))), 0, 1+((3 alp)/Nc)+((11 alp Nc)/3)-((2 alp Nf)/3), (alp/Nc)-((alp Nc)/4)}, {0, 0, 0, 0, 0, 0, 0, 6 alp-((6 alp)/(Nc^2)), 6 alp-((6 alp)/(Nc^2)), 0, ((24 alp)/Nc)-6 alp Nc, ((24 alp)/Nc)-6 alp Nc, 1-(alp/Nc)+((5 alp Nc)/3)-((2 alp Nf)/3)}}));massiveMixing=(({{((3 alp Nc)/2), 0, 0, 0, -((4 alp)/Nc)+4 alp Nc, 0, 0, 0, 0, 0, 0}, {((11 alp)/(240 Nc))+((79 alp Nc)/240), 0, 0, 0, -((149 alp)/(120 Nc))+((149 alp Nc)/120), 0, 0, 0, 0, 0, 0}, {((83 alp)/(240 Nc))+((37 alp Nc)/240), 0, alp, 0, -((197 alp)/(120 Nc))+((197 alp Nc)/120), 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 8 alp, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, -8 alp, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, ((4 alp)/Nc)-4 alp Nc, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, -((4 alp)/Nc)+4 alp Nc, 0, 0, 0, 0, 0, 0}, {alp, 0, 0, 0, 2 alp, 0, 0, 0, -8 alp Nc, 0, 0}, {alp, 0, 0, 0, 2 alp, 0, 0, 0, 0, 0, 0}, {-4 alp, 0, 0, 0, 24 alp, 0, 0, 0, 0, 0, 0}, {(alp/(2 Nc)), 0, 0, 0, (alp/Nc)-alp Nc, 0, 0, 0, 0, 0, 0}, {(alp/(2 Nc)), 0, 0, 0, (alp/Nc)-alp Nc, 0, 0, 0, 0, 0, 0}, {-((2 alp)/Nc), 4 alp, 0, 0, ((12 alp)/Nc)-12 alp Nc, 0, 0, 0, 0, 0, 0}}));massive=(({{1-((4 alp)/Nc)+2 alp Nc, 0, -4 alp, 0, ((6 alp)/Nc)-6 alp Nc, 0, 0, 0, 0, 0, 0}, {0, 1-((4 alp)/Nc)+2 alp Nc, 0, -4 alp, 0, ((6 alp)/Nc)-6 alp Nc, 0, 0, 0, 0, 0}, {0, 0, 1-((3 alp)/Nc)-((2 alp Nc)/3)+((2 alp Nf)/3), 0, 0, 0, -((3 alp)/Nc)+3 alp Nc, 0, 0, 0, 0}, {0, 0, 0, 1-((3 alp)/Nc)-((2 alp Nc)/3)+((2 alp Nf)/3), 0, 0, 0, -((3 alp)/Nc)+3 alp Nc, 0, 0, 0}, {0, 0, 0, 0, 1-((3 alp)/Nc)+3 alp Nc, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1-((3 alp)/Nc)+3 alp Nc, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 1-((3 alp)/Nc)+3 alp Nc, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 1-((3 alp)/Nc)+3 alp Nc, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1-((3 alp)/Nc)+3 alp Nc, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1-((3 alp)/Nc)+3 alp Nc, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1-((3 alp)/Nc)+3 alp Nc}}));

(* Normalise massive operators, such that operators equivalent in the mass-degenerate case carry then the same prefactor. *)
normaliseMassive=ConstantArray[1,{11}];
normaliseMassive[[{2,3,6,7,9}]]=1/Nf;
normaliseMassive[[{4,8,10}]]=1/Nf^2;
normaliseMassive[[11]]=1/Nf^3;
normaliseMassive=DiagonalMatrix[normaliseMassive];
fullMassiveMixing=ConstantArray[0,{24,24}];
fullMassiveMixing[[Range[13],Range[13]]]=massless;
fullMassiveMixing[[Range[13],13+Range[11]]]=massiveMixing.Inverse[normaliseMassive];
fullMassiveMixing[[13+Range[11],13+Range[11]]]=normaliseMassive.massive.Inverse[normaliseMassive];
pauliNorm=IdentityMatrix[24];
pauliNorm[[14,14]]=pauliNorm[[15,15]]=1/4;
fullMassiveMixing = pauliNorm.fullMassiveMixing.Inverse[pauliNorm]//Expand;

(* print out the full non-diagonal mixing matrix *)
fullMassiveMixing//MatrixForm

(* to compute the desired diagonal basis choose Nf, Nc and massless here *)
Nfvalue=4;Ncvalue=3;massless=False;GW=True;

range=Join[If[GW, Range[1,7],Range[1,13]],If[massless,{},Range[14,24]]];
jd=JordanDecomposition[D[fullMassiveMixing[[range,range]]/.Nf->Nfvalue/.Nc->Ncvalue,alp]//Expand];ord=Ordering[Diagonal[jd[[2]]]//N];
switch=IdentityMatrix[Length[jd[[2]]]][[ord]];
jnf=switch.jd[[2]].Inverse[switch];gammahat=Diagonal[jnf]/betaN/.betaN->(11/3Nc-2/3Nf)/.Nc->Ncvalue/.Nfq->Nfqvalue/.NfQ->NfQvalue/.Nf->Nfvalue;
dgammahat=(jnf-DiagonalMatrix[Diagonal[jnf]])//Expand;
invT=jd[[1]].Inverse[switch];T=Inverse[invT//N];norm={};
Do[AppendTo[norm, T[[i,Ordering[Abs[T[[i]]],-1][[1]]]]],{i,1,Length[range]}]
invT=invT.DiagonalMatrix[norm];
If[Not[massless],
isMassive={};
Do[AppendTo[isMassive,Not[AnyTrue[T[[i,Range[1,If[GW,7,13]]]]//Chop,TrueQ[# ≠0]&]]],{i,1,Length[invT]}]]
(* to express the diagonal basis terms of the non-diagonal set use: Inverse[T//N].basis[[range]]//MatrixForm *)

(* Choose matching coefficients for Overlap / Wilson QCD *)
If[GW,
(omega=ConstantArray[0,{1,18}];omega[[1,3]]=1/6;omega[[1,8]]=1;omega[[1,12]]=1/4;),
(omega=ConstantArray[0,{1,24}];
omega[[1,3]]=1/6;)]
omega.invT//N//Chop

isMassive

gammahat[[Ordering[isMassive][[Range[1,7+If[GW,0,6]]]]]]//N

gammahat[[Ordering[isMassive][[If[GW,0,6]+Range[8,18]]]]]//N

Map[TrueQ[Chop[#]==0]&,(omega.invT)[[1,Ordering[isMassive][[Range[1,7+If[GW,0,6]]]]]]]

Map[TrueQ[Chop[#]==0]&,(omega.invT)[[1,Ordering[isMassive][[If[GW,0,6]+Range[8,18]]]]]]