BeginPackage["REAP`SusyTC`", {"MixingParameterTools`MPT3x3`"}];

(* Public Interface *)

ClearAll[SCKMRotate];
SCKMRotate::usage = "SCKMRotate[x_] takes MSSM parameters contained in x and performs SCKM rotation (and some other minor stuff).";

ClearAll[CalcTreeValues];
CalcTreeValues::usage = "CalcTreeValues[x_] takes the SCKM rotated parameters contained in x and calculates the tree-level masses and mixings of sfermions, charginos, neutralinos and Higgs bosons including the \[Mu] parameter.";

ClearAll[CalcOneLoopValues]
CalcOneLoopValues::usage = "CalcOneLoopValues[x_] takes the tree-level values obtained via CalcTreeValues[..] contained in x and calculates the one-loop correction to the Higgs vev, the charged (CP-odd) Higgs mass for RGEcMSSM == True (False) and the \[Mu] parameter. Note: All other values are tree-level!";

ClearAll[OneLoopSMYd];
OneLoopSMYd::usage = "OneLoopSMYd[x] takes the one-loop values obtained via CalcOneLoopValues[..] contained in x and calculates SM Yukawa matrix Yd including its threshold corrections.";

ClearAll[OneLoopSMYe];
OneLoopSMYe::usage = "OneLoopSMYe[x] takes the one-loop values obtained via CalcOneLoopValues[..] contained in x and calculates SM Yukawa matrix Ye including its threshold corrections.";

ClearAll[OneLoopSMYu];
OneLoopSMYu::usage = "OneLoopSMYu[x] takes the one-loop values obtained via CalcOneLoopValues[..] contained in x and calculates SM Yukawa matrix Yu including its threshold corrections.";

ClearAll[OneLoopSM\[Kappa]];
OneLoopSM\[Kappa]::usage = "OneLoopSM\[Kappa][x] takes the one-loop values obtained via CalcOneLoopValues[..] contained in x and calculates wave-function renormalization threshold corrections to \[Kappa].";

ClearAll[hwave];
hwave::usage = "hwave[x] takes the one-loop values obtained via CalcOneLoopValues[..] contained in x and calculates wave-function renormalization threshold corrections to the SM-like Higgs field.";

ClearAll[Qwave];
Qwave::usage = "Qwave[x] takes the one-loop values obtained via CalcOneLoopValues[..] contained in x and calculates wave-function renormalization threshold corrections to Q.";

ClearAll[uwave];
uwave::usage = "uwave[x] takes the one-loop values obtained via CalcOneLoopValues[..] contained in x and calculates wave-function renormalization threshold corrections to u.";

ClearAll[dwave];
dwave::usage = "dwave[x] takes the one-loop values obtained via CalcOneLoopValues[..] contained in x and calculates wave-function renormalization threshold corrections to d.";

ClearAll[Lwave];
Lwave::usage = "Lwave[x] takes the one-loop values obtained via CalcOneLoopValues[..] contained in x and calculates wave-function renormalization threshold corrections to L.";

ClearAll[ewave];
ewave::usage = "ewave[x] takes the one-loop values obtained via CalcOneLoopValues[..] contained in x and calculates wave-function renormalization threshold corrections to e.";

ClearAll[GetSoftMassMixing];
GetSoftMassMixing::usage = "GetSoftMassMixing[m^2] diagonalizes the sfermion mass matrix m^2 in the EW unbroken phase and returns the sfermion masses and mixing matrix {msf, Wf}.";

ClearAll[FlavourBasis]; FlavourBasis; Protect[FlavourBasis];
FlavourBasis::usage = "Tag to signalise flavour basis.";

ClearAll[SCKMBasis]; SCKMBasis; Protect[SCKMBasis];
SCKMBasis::usage = "Tag to signalise SCKM basis.";

ClearAll[TachyonicMass]; TachyonicMass; Protect[TachyonicMass];
TachyonicMass::usage = "Exception thrown when tachyonic (or otherwise invalid) masses are encountered.";

ClearAll[NonHermitianMatrix]; NonHermitianMatrix; Protect[NonHermitianMatrix];
NonHermitianMatrix::usage = "Warning sowed when sfermion mass matrices become large non-Hermitian components in the numerics.";

ClearAll[SUSYScale];
SUSYScale::usage = "Calculates SUSY scale.";

ClearAll[TachyonFree];
TachyonFree::usage = "Checks for tachyonic states between l_high and l_low.";

ClearAll[STCGetQuarkMixing];
STCGetQuarkMixing::usage = "STCGetQuarkMixing[Yu,Yd] Diagonalizes Yu,Yd and returns {UuL,yu,UuR,UdL,yd,UdR}, s.t. Vckm=ConjugateTranspose[UuL].UdL is in PGD standard convention.";
ClearAll[STCGetLeptonMixing];
STCGetLeptonMixing::usage = "STCGetLeptonMixing[Mn,Ye] Diagonalizes Mn,Ye and returns {Un,mn,UeL,ye,UeR}, s.t. Vpmns=ConjugateTranspose[UeL].Um is in PGD standard convention.";

Begin["`Private`"];

(* Important constants *)

MZ::usage = "Z boson (pole) mass in GeV as used in the calculation."
Unprotect[MZ];
MZ = 91.1876;
Protect[MZ];

MW::usage = "W boson (pole) mass in GeV as used in the calculation."
Unprotect[MW];
MW = 80.385;
Protect[MW];

Unprotect[Loop];
Loop = 1. / (16 \[Pi]^2);
Protect[Loop];

MustBeReal[x_, name_:Null, tolerance_:10^-10] := With[{x$ = Chop[x]}, If[
  Abs[Im[x$]] > Abs[x$] tolerance,
  If[name =!= Null, Throw[{name, x}, TachyonicMass], Throw[x, TachyonicMass]],
  Re[x]
]]
SetAttributes[MustBeReal, Listable];

MustBeHermitian[x_, name_:Null, tol_:(10^-8)] := With[{x$ = 1/2 x + 1/2 ConjugateTranspose[x]},
If[
  Not[HermitianMatrixQ[x,Tolerance->tol]],
  If[name =!= Null, Sow[{name,x-ConjugateTranspose[x]//Norm},NonHermitianMatrix], Sow[x-ConjugateTranspose[x]//Norm,NonHermitianMatrix]]
  ];
  x$
]

(* Fields and field types *)

Unprotect[\[Psi]u, \[Psi]d, \[Psi]e, \[Psi]\[Nu]];
\[Psi]u = 1;
\[Psi]d = 2;
\[Psi]e = 3;
\[Psi]\[Nu] = 4;
Protect[\[Psi]u, \[Psi]d, \[Psi]e, \[Psi]\[Nu]];

Unprotect[stypes, ustypes, dstypes];
stypes = {\[Psi]u, \[Psi]d, \[Psi]e, \[Psi]\[Nu]};
ustypes = {\[Psi]u, \[Psi]\[Nu]};
dstypes = {\[Psi]d, \[Psi]e};
Protect[stypes, ustypes, dstypes];

Unprotect[\[Phi]Gp, \[Phi]Hp, \[Phi]h, \[Phi]H, \[Phi]G0, \[Phi]A];
\[Phi]Gp = 1;
\[Phi]Hp = 2;
\[Phi]h = 3;
\[Phi]H = 4;
\[Phi]G0 = 5;
\[Phi]A = 6;
Protect[\[Phi]Gp, \[Phi]Hp, \[Phi]h, \[Phi]H, \[Phi]G0, \[Phi]A];

Unprotect[Higgs0, Higgsp, Higgs0p];
Higgsp = {\[Phi]Gp, \[Phi]Hp};
Higgs0 = {\[Phi]h, \[Phi]H, \[Phi]G0, \[Phi]A};
Higgs0p = {\[Phi]h, \[Phi]H};
Protect[Higgs0, Higgsp, Higgs0p];

(* XXX for debugging *)
ClearAll[PrintAndReturn];
PrintAndReturn = Compile[{{x,_Complex}}, Print[x]; x];

(* Weak current couplings (without suffix L,R defaults to L) *)

I3[f_] := {1/2, -1/2, -1/2, 1/2}[[f]];
I3L[f_] := {1/2, -1/2, -1/2, 1/2}[[f]];
I3R[f_] := 0

Qe[f_] := {2/3, -1/3, -1, 0}[[f]];
QeL[f_] := {2/3, -1/3, -1, 0}[[f]];
QeR[f_] := {-2/3, 1/3, 1, 0}[[f]];

gfL[f_, \[Theta]W_] := I3L[f] - QeL[f] Sin[\[Theta]W]^2;
gfR[f_, \[Theta]W_] := I3R[f] - QeR[f] Sin[\[Theta]W]^2;

(* Color factors *)

NC[f_] := {3,3,1,1}[[f]];

(* Helpful definitions *)

Id3 = IdentityMatrix[3];
Zero3 = 0 Id3;
Complexify = Compile[{{x, _Complex}}, (1.0 + 0.0 I) x];
ComplexMatrix = Compile[{{x, _Complex, 2}}, (1.0 + 0.0 I) x];

(* Loop Functions *)

A0 = Compile[{{m, _Real}, {Q, _Real}}, If[Abs[m] > $MachineEpsilon, m^2 (1 - Log[m^2/Q^2]), 0]];

\[Epsilon] = $MachineEpsilon;
fb = Compile[{{x, _Complex}}, Log[1 - x] - x Log[1 - 1/x] - 1];
B0 = Compile[{{p, _Real}, {m1, _Real}, {m2, _Real}, {Q, _Real}},
  With[{s = p^2 - m2^2 + m1^2},
    If[p != 0,
      With[{
          xp = (s + Sqrt[s^2 - 4 p^2 (m1^2 - I \[Epsilon])])/(2 p^2),
          xm = (s - Sqrt[s^2 - 4 p^2 (m1^2 - I \[Epsilon])])/(2 p^2)
        },
        -Log[p^2/Q^2] - fb[xp] - fb[xm]
      ],
      (*else*)
      If[m1 != m2,
        1 + Log[Q^2/m2^2] + m1^2/(m1^2 - m2^2) Log[m2^2/m1^2],
        (*else*)
        Log[Q^2/m2^2]
      ]
    ]
  ]
];

B22 = Compile[{{p, _Real}, {m1, _Real}, {m2, _Real}, {Q, _Real}},
  1/6 * (
    1/2*(A0[m1, Q] + A0[m2, Q]) 
    + (m1^2 + m2^2 - 1/2 p^2) B0[p, m1, m2, Q] 
    + (m2^2 - m1^2)/(2 p^2)*(A0[m2, Q] - A0[m1, Q] - (m2^2 - m1^2) B0[p, m1, m2, Q]) 
    + m1^2 + m2^2 - 1/3 p^2
  )
];
B22t = Compile[{p, m1, m2, Q}, B22[p, m1, m2, Q] - 1/4 (A0[m1, Q] + A0[m2, Q])];
G = Compile[{p, m1, m2, Q}, (p^2 - m1^2 - m2^2) B0[p, m1, m2, Q] - (A0[m1, Q] + A0[m2, Q])];
H = Compile[{p, m1, m2, Q}, 4 B22[p, m1, m2, Q] + G[p, m1, m2, Q]];
F = Compile[{p, m1, m2, Q}, A0[m1, Q] - 2 A0[m2, Q] - (2 p^2 + 2 m1^2 - m2^2) B0[p, m1, m2, Q]];


(* Mass basis rotations *)

(* Takagi decomposition algorithm *)
Phase[x_] := If[Chop[x] =!= 0, x/Abs[x], 1];
UTakagi[A_, i_, j_, n_] := Block[{
    U = IdentityMatrix[n],
    c, t, \[CapitalDelta], D, ei\[CurlyPhi]
  },
  (* elementary decomposition step *)
  If[Chop[A[[i, j]]] === 0, Return[U]];
  ei\[CurlyPhi] = Phase[Conjugate[A[[i, i]]] A[[i, j]] + A[[j, j]] Conjugate[A[[i, j]]]];
  \[CapitalDelta] = 1/2 (ei\[CurlyPhi] A[[i, i]] - Conjugate[ei\[CurlyPhi]] A[[j, j]]);
  D = Sqrt[\[CapitalDelta]^2 + A[[i, j]]^2];
  If[Abs[\[CapitalDelta] + D] < Abs[\[CapitalDelta] - D], D = -D];
  t = Re[A[[i, j]]/(\[CapitalDelta] + D)];
  c = 1/Sqrt[1 + t^2];
  U[[i, i]] = c;
  U[[i, j]] = t c ei\[CurlyPhi];
  U[[j, i]] = -t c Conjugate[ei\[CurlyPhi]];
  U[[j, j]] = c;
  Conjugate[U]
];
  
TakagiDecompositionStep[u$_, m$_] := Block[{d = Length[m$], U, u, m},
  U = u$;
  m = m$;
  Do[
    Do[
      u = UTakagi[m, i, j, d];
      m = u.m.Transpose[u];
      U = u.U, 
      {j, i - 1, 1, -1}
    ],
    {i, d, 1, -1}
  ];
  {U, m}
];

TakagiDecomposition[M_] := Block[{u, d, \[CurlyPhi], order},
  {u, d} = FixedPoint[
    TakagiDecompositionStep[#[[1]], #[[2]]] &,
    {IdentityMatrix[Length[M]], M}
  ];
  \[CurlyPhi] = DiagonalMatrix[Exp[-(I/2) Arg[Diagonal[d]]]];
  u = \[CurlyPhi].u;
  d = Abs[Diagonal[d]];
  order = Ordering[d];
  {d[[order]], u[[order]]}
];

GetNeutralinoMixing[M1_, M2_, \[Mu]_, gp_, g2_, v_, \[Beta]_] := Block[
  {
    s2\[Theta]W, c2\[Theta]W, s\[Theta]W, c\[Theta]W, 
    c\[Beta] = Cos[\[Beta]], s\[Beta] = Sin[\[Beta]], 
    M, m\[Chi]0, N, Ni, ni,
    MZ = 1/2 v Sqrt[gp^2 + g2^2]
  },
  s2\[Theta]W = gp^2/(gp^2 + g2^2);
  c2\[Theta]W = 1 - s2\[Theta]W;
  s\[Theta]W = Sqrt[s2\[Theta]W];
  c\[Theta]W = Sqrt[c2\[Theta]W];
  M = {
    {M1, 0, -MZ c\[Beta] s\[Theta]W, MZ s\[Beta] s\[Theta]W}, 
    {0, M2, MZ c\[Beta] c\[Theta]W, -MZ s\[Beta] c\[Theta]W}, 
    {-MZ c\[Beta] s\[Theta]W, MZ c\[Beta] c\[Theta]W, 0, -\[Mu]}, 
    {MZ s\[Beta] s\[Theta]W, -MZ s\[Beta] c\[Theta]W, -\[Mu], 0}
  };
  {m\[Chi]0, N} = TakagiDecomposition[M];
  
  {m\[Chi]0, N}
];

GetCharginoMixing[M2_, \[Mu]_, g2_, v_, \[Beta]_] := Block[
  {
    c\[Beta] = Cos[\[Beta]], s\[Beta] = Sin[\[Beta]],
    M, m\[Chi]p,
    U, V, Ui, Vi, ni,
    MW = v/2 g2
  },
  M = {{M2, Sqrt[2] MW s\[Beta]}, {Sqrt[2] MW c\[Beta], \[Mu]}};
  {U, m\[Chi]p, V} = SingularValueDecomposition[M];
  U = Transpose[U];
  V = ConjugateTranspose[V];
  Reverse /@ {Diagonal[m\[Chi]p], U, V}
];

(*LR convention and SCKM input*)
GetSfermionMixing[ft_, ML2_?MatrixQ, MR2_?MatrixQ, \[Mu]_, A_?MatrixQ, y_?MatrixQ, gp_, g2_, v_, \[Beta]_] := Module[{
    m2l, m2r,
    s2\[Theta]W, c2\[Theta]W, s\[Theta]W, c\[Theta]W,
    c\[Beta] = Cos[\[Beta]], s\[Beta] = Sin[\[Beta]],
    M2, msf, Wf,
    sc\[Beta], cs\[Beta],
    MZ = 1/2 v Sqrt[gp^2 + g2^2]
  },
  s2\[Theta]W = 1 - g2^2/(gp^2 + g2^2);
  c2\[Theta]W = 1 - s2\[Theta]W;
  s\[Theta]W = Sqrt[s2\[Theta]W];
  c\[Theta]W = Sqrt[c2\[Theta]W];
  If[MemberQ[{\[Psi]u, \[Psi]\[Nu]}, ft],
    sc\[Beta] = s\[Beta];
    cs\[Beta] = c\[Beta],
    (* else *)
    sc\[Beta] = c\[Beta];
    cs\[Beta] = s\[Beta]
  ];
  m2l = (v/Sqrt[2] sc\[Beta])^2 Conjugate[y].Transpose[y];
  m2r = (v/Sqrt[2] sc\[Beta])^2 Transpose[y].Conjugate[y];
  
  If[ft == \[Psi]\[Nu],
    M2 = ML2 + MZ^2 (I3[ft] - Qe[ft] s2\[Theta]W) Cos[2 \[Beta]] Id3,
    (* else *)
    M2 = ArrayFlatten[
     {
      {
        ML2 + MZ^2 (I3[ft] - Qe[ft] s2\[Theta]W) Cos[2 \[Beta]] Id3 + m2l,
        v/Sqrt[2] sc\[Beta] ConjugateTranspose[A] - \[Mu] v / Sqrt[2] Conjugate[y] cs\[Beta]
      },
      {
        v/Sqrt[2] sc\[Beta] A - Conjugate[\[Mu]] v/Sqrt[2] Transpose[y] cs\[Beta],
        MR2 + MZ^2 Qe[ft] Cos[2 \[Beta]] s2\[Theta]W Id3 + m2r
      }
     }
    ]
  ];
  M2=MustBeHermitian[M2];
  {msf, Wf} = Eigensystem[M2];
  msf = MustBeReal[Reverse[Sqrt[msf]], "sfermion mass "<>Switch[ft,1,"u",2,"d",3,"e",4,"\[Nu]"]];
  Wf = Transpose[Reverse[Wf]];
  {msf, Wf, M2}
];

GetCPEvenHiggsMixing[MA_, gp_, g2_, v_, \[Beta]_] := Block[
  {
    MHsq, mhsq, MZ = 1/2 v Sqrt[gp^2 + g2^2],
    tan\[Alpha], X
  },
  M = {
    {
      MZ^2 Cos[\[Beta]]^2 + MA^2 Sin[\[Beta]]^2,
      -(MZ^2 + MA^2) Sin[\[Beta]] Cos[\[Beta]]
    }, 
    {
      -(MZ^2 + MA^2) Sin[\[Beta]] Cos[\[Beta]],
      MZ^2 Sin[\[Beta]]^2 + MA^2 Cos[\[Beta]]^2
    }
  };
  {{MHsq, mhsq}, X} = Eigensystem[M];
  tan\[Alpha] = -X[[2, 1]]/X[[2, 2]];
  {
    MustBeReal[Sqrt[mhsq], "mh"],
    MustBeReal[Sqrt[MHsq], "mH"],
    tan\[Alpha]
  }
];

(* Couplings *)

(* EWkino Z couplings *)

Off[Part::pspec, Part::pkspec1];

aZ0[g2_, \[Theta]W_] = DiagonalMatrix[{0, 0, g2/(2 Cos[\[Theta]W]), -g2/(2 Cos[\[Theta]W])}];
bZ0[g2_, \[Theta]W_] = -aZ0[g2, \[Theta]W];

fZ0[i_, j_, N_, g2_, \[Theta]W_] = (
  Abs[
    Sum[Conjugate[N[[i, k]]] N[[j, l]] aZ0[g2, \[Theta]W][[k, l]], {k, 4}, {l, 4}]
  ]^2
  + Abs[
    Sum[N[[i, k]] Conjugate[N[[j, l]]] bZ0[g2, \[Theta]W][[k, l]], {k, 4}, {l, 4}]
  ]^2
);
  
gZ0[i_, j_, N_, g2_, \[Theta]W_] = 2 Re[
  Conjugate[Sum[N[[i, k]] Conjugate[N[[j, l]]] bZ0[g2, \[Theta]W][[k, l]], {k, 4}, {l, 4}]] 
  *
  Sum[Conjugate[N[[i, k]]] N[[j, l]] aZ0[g2, \[Theta]W][[k, l]], {k, 4}, {l, 4}]
];

aZp[g2_, \[Theta]W_] := DiagonalMatrix[{ g2 Cos[\[Theta]W], g2 Cos[2 \[Theta]W]/(2 Cos[\[Theta]W])}];
bZp[g2_, \[Theta]W_] := aZp[g2, \[Theta]W];

fZp[i_, j_, U_, V_, g2_, \[Theta]W_] = (
  Abs[
    Sum[Conjugate[V[[i, k]]] V[[j, l]] aZp[g2, \[Theta]W][[k, l]], {k, 2}, {l, 2}]
  ]^2 
  + Abs[
    Sum[U[[i, k]] Conjugate[U[[j, l]]] bZp[g2, \[Theta]W][[k, l]], {k, 2}, {l, 2}]
  ]^2
);
gZp[i_, j_, U_, V_, g2_, \[Theta]W_] = 2 Re[
  Conjugate[Sum[U[[i, k]] Conjugate[U[[j, l]]] bZp[g2, \[Theta]W][[k, l]], {k, 2}, {l, 2}]]
  *
  Sum[Conjugate[V[[i, k]]] V[[j, l]] aZp[g2, \[Theta]W][[k, l]], {k, 2}, {l, 2}]
];

(* compile the highest level functions *)

fZ0 = Compile[{{i, _Integer}, {j, _Integer}, {N, _Complex, 2}, {g2, _Real}, {\[Theta]W, _Real}},
  Evaluate[fZ0[i, j, N, g2, \[Theta]W]]
];
gZ0 = Compile[{{i, _Integer}, {j, _Integer}, {N, _Complex, 2}, {g2, _Real}, {\[Theta]W, _Real}}, 
  Evaluate[gZ0[i, j, N, g2, \[Theta]W]]
];
fZp = Compile[{{i, _Integer}, {j, _Integer}, {U, _Complex, 2}, {V, _Complex, 2}, {g2, _Real}, {\[Theta]W, _Real}},
  Evaluate[fZp[i, j, U, V, g2, \[Theta]W]]
];
gZp = Compile[{{i, _Integer}, {j, _Integer}, {U, _Complex, 2}, {V, _Complex, 2}, {g2, _Real}, {\[Theta]W, _Real}},
  Evaluate[gZp[i, j, U, V, g2, \[Theta]W]]
]

On[Part::pspec, Part::pkspec1];

(* EWkino Higgs couplings *)

Off[Part::pspec, Part::pkspec1];

as0[gp_, g2_] = Block[{ai = ConstantArray[0, {4, 4, 2}]}, 
  ai[[1, 4, 2]] = gp/2;
  ai[[4, 1, 2]] = ai[[1, 4, 2]];
  ai[[1, 3, 1]] = -ai[[1, 4, 2]];
  ai[[3, 1, 1]] = -ai[[1, 4, 2]];
  ai[[2, 3, 1]] = g2/2;
  ai[[3, 2, 1]] = ai[[2, 3, 1]];
  ai[[2, 4, 2]] = -ai[[2, 3, 1]];
  ai[[4, 2, 2]] = -ai[[2, 3, 1]];
  
  ai
];
bs0[gp_, g2_] = as0[gp, g2];

ap0[gp_, g2_] = Block[{ai = ConstantArray[0, {4, 4, 2}]}, 
  ai[[1, 4, 2]] = -gp/2; 
  ai[[4, 1, 2]] = ai[[1, 4, 2]];
  ai[[1, 3, 1]] = ai[[1, 4, 2]];
  ai[[3, 1, 1]] = ai[[1, 4, 2]];
  ai[[2, 3, 1]] = g2/2;
  ai[[3, 2, 1]] = ai[[2, 3, 1]];
  ai[[2, 4, 2]] = ai[[2, 3, 1]];
  ai[[4, 2, 2]] = ai[[2, 3, 1]];
  
  ai
];
bp0[gp_, g2_] = -ap0[gp, g2];

fA0[i_, j_, N_, gp_, g2_, \[Beta]_] = (
  Abs[
    -Sin[\[Beta]] Sum[
      Conjugate[N[[i, k]]] Conjugate[N[[j, l]]] ap0[gp, g2][[k, l, 1]], {k, 4}, {l, 4}
    ] + 
    Cos[\[Beta]] Sum[
      Conjugate[N[[i, k]]] Conjugate[N[[j, l]]] ap0[gp, g2][[k, l, 2]], {k, 4}, {l, 4}]
  ]^2 + 
  Abs[
    -Sin[\[Beta]] Sum[
      N[[i, k]] N[[j, l]] bp0[gp, g2][[k, l, 1]], {k, 4}, {l, 4}] + 
    Cos[\[Beta]] Sum[
      N[[i, k]] N[[j, l]] bp0[gp, g2][[k, l, 2]], {k, 4}, {l, 4}]
  ]^2
);
gA0[i_, j_, N_, gp_, g2_, \[Beta]_] = 2 Re[
  Conjugate[
    -Sin[\[Beta]] Sum[
      N[[i, k]] N[[j, l]] bp0[gp, g2][[k, l, 1]], {k, 4}, {l, 4}
    ] 
    + Cos[\[Beta]] Sum[
      N[[i, k]] N[[j, l]] bp0[gp, g2][[k, l, 2]], {k, 4}, {l, 4}
    ]
  ]
  *
  (
    -Sin[\[Beta]] Sum[
      Conjugate[N[[i, k]]] Conjugate[N[[j, l]]] ap0[gp, g2][[k, l, 1]], {k, 4}, {l, 4}
    ] 
    + Cos[\[Beta]] Sum[
      Conjugate[N[[i, k]]] Conjugate[N[[j, l]]] ap0[gp, g2][[k, l, 2]], {k, 4}, {l, 4}
    ]
  )
];

asp[g2_] = Block[{ai = ConstantArray[0, {2, 2, 2}]},
  ai[[1, 2, 1]] = g2/Sqrt[2]; 
  ai[[2, 1, 2]] = ai[[1, 2, 1]]; 
  
  ai
];
bsp[g2_] = Block[{ai = ConstantArray[0, {2, 2, 2}]},
  ai[[2, 1, 1]] = g2/Sqrt[2]; 
  ai[[1, 2, 2]] = ai[[2, 1, 1]]; 
  
  ai
];

app[g2_] = Block[{ai = ConstantArray[0, {2, 2, 2}]},
  ai[[1, 2, 1]] = g2/Sqrt[2];
  ai[[2, 1, 2]] = -ai[[1, 2, 1]]; 
  
  ai
];
bpp[g2_] = Block[{ai = ConstantArray[0, {2, 2, 2}]},
  ai[[2, 1, 1]] = -g2/Sqrt[2];
  ai[[1, 2, 2]] = -ai[[2, 1, 1]];
  
  ai
];
   
fAp[i_, j_, U_, V_, g2_, \[Beta]_] = (
  Abs[
    -Sin[\[Beta]] Sum[
      Conjugate[V[[i, k]]] Conjugate[U[[j, l]]] app[g2][[k, l, 1]], {k, 2}, {l, 2}
    ] 
    + Cos[\[Beta]] Sum[
       Conjugate[V[[i, k]]] Conjugate[U[[j, l]]] app[g2][[k, l, 2]], {k, 2}, {l, 2}
    ]
  ]^2 +
  Abs[
    -Sin[\[Beta]] Sum[
      U[[i, k]] V[[j, l]] bpp[g2][[k, l, 1]], {k, 2}, {l, 2}
    ] 
    + Cos[\[Beta]] Sum[
      U[[i, k]] V[[j, l]] bpp[g2][[k, l, 2]], {k, 2}, {l, 2}
    ]
  ]^2
);
gAp[i_, j_, U_, V_, g2_, \[Beta]_] = 2 Re[
  Conjugate[
    -Sin[\[Beta]] Sum[
      U[[i, k]] V[[j, l]] bpp[g2][[k, l, 1]], {k, 2}, {l, 2}
    ] 
    + Cos[\[Beta]] Sum[
      U[[i, k]] V[[j, l]] bpp[g2][[k, l, 2]], {k, 2}, {l, 2}
    ]
  ]
  *
  (
    -Sin[\[Beta]] Sum[
      Conjugate[V[[i, k]]] Conjugate[U[[j, l]]] app[g2][[k, l, 1]], {k, 2}, {l, 2}
    ]
    + Cos[\[Beta]] Sum[
      Conjugate[V[[i, k]]] Conjugate[U[[j, l]]] app[g2][[k, l, 2]], {k, 2}, {l, 2}
    ]
  )
];

a0p[gp_,g2_]:=Block[{ai=ConstantArray[0,{4,2,2}]},
  ai[[1,2,1]]=gp/Sqrt[2];
  ai[[2,2,1]]=g2/Sqrt[2];
  ai[[3,1,1]]=-g2;

  ai
];
b0p[gp_,g2_]:=Block[{ai=ConstantArray[0,{4,2,2}]},
  ai[[1,2,2]]=gp/Sqrt[2];
  ai[[2,2,2]]=g2/Sqrt[2];
  ai[[4,1,2]]=g2;

  ai
];

fHp[i_, j_, N_, U_, V_, gp_, g2_, \[Beta]_] = (
  Abs[
    -Sin[\[Beta]] Sum[
      Conjugate[N[[i, k]]] Conjugate[U[[j, l]]] a0p[gp, g2][[k, l, 1]], {k, 4}, {l, 2}
    ]
    + Cos[\[Beta]] Sum[
      Conjugate[N[[i, k]]] Conjugate[U[[j, l]]] a0p[gp, g2][[k, l, 2]], {k, 4}, {l, 2}
    ]
  ]^2 +
  Abs[
    -Sin[\[Beta]] Sum[
      N[[i, k]] V[[j, l]] b0p[gp, g2][[k, l, 1]], {k, 4}, {l, 2}
    ]
    + Cos[\[Beta]] Sum[
      N[[i, k]] V[[j, l]] b0p[gp, g2][[k, l, 2]], {k, 4}, {l, 2}]
  ]^2
);
gHp[i_, j_, N_, U_, V_, gp_, g2_, \[Beta]_] = 2 Re[
  Conjugate[
    -Sin[\[Beta]] Sum[
      N[[i, k]] V[[j, l]] b0p[gp, g2][[k, l, 1]], {k, 4}, {l, 2}
    ]
    + Cos[\[Beta]] Sum[
      N[[i, k]] V[[j, l]] b0p[gp, g2][[k, l, 2]], {k, 4}, {l, 2}
    ]
  ]
  *
  (
    -Sin[\[Beta]] Sum[
      Conjugate[N[[i, k]]] Conjugate[U[[j, l]]] a0p[gp, g2][[k, l, 1]], {k, 4}, {l, 2}
    ]
    + Cos[\[Beta]] Sum[
      Conjugate[N[[i, k]]] Conjugate[U[[j, l]]] a0p[gp, g2][[k, l, 2]], {k, 4}, {l, 2}
    ]
  )
];

(* compile the functions *)
fA0 = Compile[{{i, _Integer}, {j, _Integer}, {N, _Complex, 2}, {gp, _Real}, {g2, _Real}, {\[Theta]W, _Real}},
  Evaluate[fA0[i, j, N, gp, g2, \[Theta]W]]
];
gA0 = Compile[{{i, _Integer}, {j, _Integer}, {N, _Complex, 2}, {gp, _Real}, {g2, _Real}, {\[Theta]W, _Real}},
  Evaluate[gA0[i, j, N, gp, g2, \[Theta]W]]
];
fAp = Compile[{{i, _Integer}, {j, _Integer}, {U, _Complex, 2}, {V, _Complex, 2}, {g2, _Real}, {\[Beta], _Real}},
  Evaluate[fAp[i, j, U, V, g2, \[Beta]]]
];
gAp = Compile[{{i, _Integer}, {j, _Integer}, {U, _Complex, 2}, {V, _Complex, 2}, {g2, _Real}, {\[Beta], _Real}},
  Evaluate[gAp[i, j, U, V, g2, \[Beta]]]
];
fHp = Compile[{{i, _Integer}, {j, _Integer}, {N, _Complex, 2}, {U, _Complex, 2}, {V, _Complex, 2}, {gp, _Real}, {g2, _Real}, {\[Beta], _Real}},
  Evaluate[fHp[i, j, N, U, V, gp, g2, \[Beta]]]
];
gHp = Compile[{{i, _Integer}, {j, _Integer}, {N, _Complex, 2}, {U, _Complex, 2}, {V, _Complex, 2}, {gp, _Real}, {g2, _Real}, {\[Beta], _Real}},
  Evaluate[gHp[i, j, N, U, V, gp, g2, \[Beta]]]
];

On[Part::pspec, Part::pkspec1];

(* Higgs Higgs couplings *)

\[Lambda]AA = With[{Rot\[Alpha] = {{Cos[\[Alpha]], -Sin[\[Alpha]]}, {Sin[\[Alpha]], Cos[\[Alpha]]}}},
  Compile[{
    {x, _Integer}, {g2, _Real}, {\[Theta]W, _Real}, {\[Alpha], _Real}, {\[Beta], _Real}
  }, 
  With[{
    M = Transpose[Rot\[Alpha]]
      .DiagonalMatrix[{-Cos[2 \[Beta]], Cos[2 \[Beta]]}]
      .Rot\[Alpha]
  },
  g2^2/(4 Cos[\[Theta]W]^2) *
    {
      (* Gp *)
      Cos[\[Theta]W]^2 (1 + Sin[2 \[Beta]]^2) - Sin[\[Theta]W]^2 Cos[2 \[Beta]]^2,
      (* Hp *)
      Cos[2 \[Beta]]^2,
      (* h *)
      M[[2, 2]],
      (* H *)
      M[[1, 1]],
      (* G0 *)
      3 Sin[2 \[Beta]]^2 - 1,
      (* A *)
      3 Cos[2 \[Beta]]^2
    }[[x]]
  ]
]];
  
\[Lambda]A = With[{Rot\[Alpha] = {{Cos[\[Alpha]], -Sin[\[Alpha]]}, {Sin[\[Alpha]], Cos[\[Alpha]]}}},
  Compile[{
    {x, _Integer}, {y, _Integer}, 
    {gp, _Real}, {g2, _Real}, {v, _Real}, {\[Theta]W, _Real},
    {\[Alpha], _Real}, {\[Beta], _Real}
  },
  Module[{
    MZ = 1/2 v Sqrt[gp^2 + g2^2],
    MA = Transpose[Rot\[Alpha]]
      .{
        {-Cos[2 \[Beta]] Cos[\[Beta]]},
        {Cos[2 \[Beta]] Sin[\[Beta]]}
       },
    MG = Transpose[Rot\[Alpha]]
      .{
        {-Sin[2 \[Beta]] Cos[\[Beta]]},
        {Sin[2 \[Beta]] Sin[\[Beta]]}
       }
  },
  (* order is \[Phi]Gp, \[Phi]Hp, \[Phi]h, \[Phi]H, \[Phi]G0, \[Phi]A *)
  g2 MZ/(2 Cos[\[Theta]W]) *
    If[x == 6 (* \[Phi]A *) || y == 6 (* \[Phi]A *),
      If[x == 4 (* \[Phi]H *) || y == 4 (* \[Phi]H *), 
        MA[[1, 1]],
        (*else x\[Equal]h||y\[Equal]h *)
        MA[[2, 1]]
      ],
      (*else x\[Equal]G0||y\[Equal]G0*)
    
      If[x == 4 (* \[Phi]H *) || y == 4 (* \[Phi]H *), 
        MG[[1, 1]],
        (*else x\[Equal]h||y\[Equal]h*)
        MG[[2, 1]]
      ]
    ]
  ]
]];

\[Lambda]HpH0Hm = With[{Rot\[Alpha] = {{Cos[\[Alpha]], -Sin[\[Alpha]]}, {Sin[\[Alpha]], Cos[\[Alpha]]}}},
  Compile[{
    {x, _Integer}, {y, _Integer}, {gp, _Real}, {g2, _Real}, {v, _Real},
    {\[Theta]W, _Real}, {\[Alpha], _Real}, {\[Beta], _Real}
  },
  Module[{
    MZ = 1/2 v Sqrt[gp^2 + g2^2], \[Lambda], \[Lambda]f
    },
    \[Lambda] = g2^2 MZ/(2 Cos[\[Theta]W]) * 
      {
        {(*GP*)
        {-Sin[2 \[Beta]] Cos[\[Beta]] + Cos[\[Theta]W]^2 Sin[\[Beta]]},
        {Sin[2 \[Beta]] Sin[\[Beta]] + Cos[\[Theta]W]^2 Cos[\[Beta]]}
        },
        {(* Hp *)
        {-Cos[2 \[Beta]] Cos[\[Beta]] + 2 Cos[\[Theta]W]^2 Cos[\[Beta]]},
        {Cos[2 \[Beta]] Sin[\[Beta]] + 2Cos[\[Theta]W]^2 Sin[\[Beta]]}
        }
      }[[y]];
    \[Lambda]f = Transpose[Rot\[Alpha]].\[Lambda];
    {
      (*Gp*) 0,
      (*Hp*) 0,
      (*h*) \[Lambda]f[[2,1]],
      (*H*) \[Lambda]f[[1,1]]
    }[[x]]
  ]
]];

\[Lambda]HpHmH0H0 = With[{Rot\[Alpha] = {{Cos[\[Alpha]], -Sin[\[Alpha]]}, {Sin[\[Alpha]], Cos[\[Alpha]]}}},
  Compile[{
    {x, _Integer}, {g2, _Real},
    {\[Theta]W, _Real}, {\[Alpha], _Real}, {\[Beta], _Real}
  },
  Module[{\[Lambda]},
    \[Lambda] = Transpose[Rot\[Alpha]].
      {
        {Cos[\[Theta]W]^2 - Sin[\[Theta]W]^2 Cos[2 \[Beta]], Cos[\[Theta]W]^2 Sin[2 \[Beta]]},
        {Cos[\[Theta]W]^2 Sin[2 \[Beta]], Cos[\[Theta]W]^2 + Sin[\[Theta]W]^2 Cos[2 \[Beta]]}
      }.Rot\[Alpha];
    g2^2/(4 Cos[\[Theta]W]^2) {
      (*Gp*) 0,
      (*Hp*) 0,
      (*h*) \[Lambda][[2,2]],
      (*H*) \[Lambda][[1,1]],
      (*G0*) Cos[\[Theta]W]^2 (1 + Sin[2 \[Beta]]^2) - Sin[\[Theta]W]^2 Cos[2 \[Beta]]^2,
      (*A*) Cos[2 \[Beta]]^2
    }[[x]]
  ]
]];

\[Lambda]HpHmHpHm = 
  Compile[{
    {x, _Integer}, {g2, _Real},
    {\[Theta]W, _Real}, {\[Beta], _Real}
  },
  g2^2/(4 Cos[\[Theta]W]^2) * 
  {
    (*Gp*) 2 Sin[2 \[Beta]]^2-1,
    (*Hp*) 2 Cos[2 \[Beta]]^2
  }[[x]]
]

(* Higgs sfermion Couplings *)

(*ComposeSfermion6x6MassMatrix = Compile[
  {{M, _Complex, 4}},
  Module[{M$},
    M$ = (1.0 + I 0.0) {
      {1, 0, 0, 0, 0, 0}, 
      {0, 1, 0, 0, 0, 0}, 
      {0, 0, 1, 0, 0, 0}, 
      {0, 0, 0, 1, 0, 0}, 
      {0, 0, 0, 0, 1, 0}, 
      {0, 0, 0, 0, 0, 1}
    };
    
    Do[
      Do[
        M$[[i, j]] = M[[1, 1]][[i, j]];
        M$[[i + 3, j]] = M[[2, 1]][[i, j]];
        M$[[i, j + 3]] = M[[1, 2]][[i, j]];
        M$[[i + 3, j + 3]] = M[[2, 2]][[i, j]];
      , {j, 3}
      ]
    , {i, 3}
    ];
    M$
  ]
];*)

ComposeSfermion6x6MassMatrix  = Compile[
  {{M, _Complex, 4}},
  {
    {
      M[[1, 1, 1, 1]], M[[1, 1, 1, 2]], M[[1, 1, 1, 3]], 
      M[[1, 2, 1, 1]], M[[1, 2, 1, 2]], M[[1, 2, 1, 3]]
    },
    {
      M[[1, 1, 2, 1]], M[[1, 1, 2, 2]], M[[1, 1, 2, 3]], 
      M[[1, 2, 2, 1]], M[[1, 2, 2, 2]], M[[1, 2, 2, 3]]
    }, 
    {
      M[[1, 1, 3, 1]], M[[1, 1, 3, 2]], M[[1, 1, 3, 3]], 
      M[[1, 2, 3, 1]], M[[1, 2, 3, 2]], M[[1, 2, 3, 3]]
    }, 
    {
      M[[2, 1, 1, 1]], M[[2, 1, 1, 2]], M[[2, 1, 1, 3]], 
      M[[2, 2, 1, 1]], M[[2, 2, 1, 2]], M[[2, 2, 1, 3]]
    }, 
    {
      M[[2, 1, 2, 1]], M[[2, 1, 2, 2]], M[[2, 1, 2, 3]], 
      M[[2, 2, 2, 1]], M[[2, 2, 2, 2]], M[[2, 2, 2, 3]]
    }, 
    {
      M[[2, 1, 3, 1]], M[[2, 1, 3, 2]], M[[2, 1, 3, 3]], 
      M[[2, 2, 3, 1]], M[[2, 2, 3, 2]], M[[2, 2, 3, 3]]
    }
  }
];

ComposeSfermion3x6MassMatrix  = Compile[
  {{M, _Complex, 4}},
  {
    {
      M[[1, 1, 1, 1]], M[[1, 1, 1, 2]], M[[1, 1, 1, 3]], 
      M[[1, 2, 1, 1]], M[[1, 2, 1, 2]], M[[1, 2, 1, 3]]
    },
    {
      M[[1, 1, 2, 1]], M[[1, 1, 2, 2]], M[[1, 1, 2, 3]], 
      M[[1, 2, 2, 1]], M[[1, 2, 2, 2]], M[[1, 2, 2, 3]]
    }, 
    {
      M[[1, 1, 3, 1]], M[[1, 1, 3, 2]], M[[1, 1, 3, 3]], 
      M[[1, 2, 3, 1]], M[[1, 2, 3, 2]], M[[1, 2, 3, 3]]
    }
  }
];
  
Diag3 = Compile[{{e, _Complex, 1}}, 
  {{e[[1]], 0, 0}, {0, e[[2]], 0}, {0, 0, e[[3]]}} (1.0 + I 0.0)
];

(* this pads \[Lambda]X matrices to 6 flavours for use with neutrinos *)
(* note that the formulae never use the zeros -- they are only there to make Mathematica happy *)
Pad\[Nu] = Compile[{{M, _Complex, 2}},
  With[{Zero3 = {
      {0. + 0. I, 0. + 0. I, 0. + 0. I}, 
      {0. + 0. I, 0. + 0. I, 0. + 0. I}, 
      {0. + 0. I, 0. + 0. I, 0. + 0. I}
    }},
    ComposeSfermion6x6MassMatrix[{{M, Zero3}, {Zero3, Zero3}}]
  ]
];

(* The same thing, but for padding (3 x 6) to (6 x 6) *)
Pad3x6\[Nu] = Compile[{{M, _Complex, 2}},
  Module[{ret = {
      {0. + 0. I, 0. + 0. I, 0. + 0. I,0. + 0. I, 0. + 0. I, 0. + 0. I}, 
      {0. + 0. I, 0. + 0. I, 0. + 0. I,0. + 0. I, 0. + 0. I, 0. + 0. I}, 
      {0. + 0. I, 0. + 0. I, 0. + 0. I,0. + 0. I, 0. + 0. I, 0. + 0. I},
      {0. + 0. I, 0. + 0. I, 0. + 0. I,0. + 0. I, 0. + 0. I, 0. + 0. I},
      {0. + 0. I, 0. + 0. I, 0. + 0. I,0. + 0. I, 0. + 0. I, 0. + 0. I},
      {0. + 0. I, 0. + 0. I, 0. + 0. I,0. + 0. I, 0. + 0. I, 0. + 0. I}
    }},
  ret[[1,1]] = M[[1,1]];
  ret[[1,2]] = M[[1,2]];
  ret[[1,3]] = M[[1,3]];
  ret[[1,4]] = M[[1,4]];
  ret[[1,5]] = M[[1,5]];
  ret[[1,6]] = M[[1,6]];
  ret[[2,1]] = M[[2,1]];
  ret[[2,2]] = M[[2,2]];
  ret[[2,3]] = M[[2,3]];
  ret[[2,4]] = M[[2,4]];
  ret[[2,5]] = M[[2,5]];
  ret[[2,6]] = M[[2,6]];
  ret[[3,1]] = M[[3,1]];
  ret[[3,2]] = M[[3,2]];
  ret[[3,3]] = M[[3,3]];
  ret[[3,4]] = M[[3,4]];
  ret[[3,5]] = M[[3,5]];
  ret[[3,6]] = M[[3,6]];
  ret
  ]
];

\[Lambda]d = Compile[{
    {Wude, _Complex, 3}, {W\[Nu], _Complex, 2},
    {gp, _Real}, {g2, _Real}, {v, _Real},
    {yf, _Real, 2}, {Af, _Complex, 3}, {mf, _Real, 2},
    {\[Theta]W, _Real}, {\[Beta], _Real}, {\[Mu], _Complex}
  }, Module[{
    Wu,Wd,We,yu,yd,ye,y\[Nu],Au,Ad,Ae,A\[Nu],mu,md,me,m\[Nu]D,f,ret,Id3, retu,retd,rete,ret\[Nu],
    MZ = 1/2 v Sqrt[gp^2 + g2^2]
  },
  Id3 = {{1. + 0. I, 0. + 0. I, 0. + 0. I}, {0. + 0. I, 1. + 0. I, 0. + 0. I}, {0. + 0. I, 0. + 0. I, 1. + 0. I}};
  {Wu,Wd,We} = {Id3, Id3, Id3};
  ret = {Id3, Id3, Id3};
  
  {Wu,Wd,We} = Wude;
  {yu,yd,ye,y\[Nu]} = yf;
  {Au,Ad,Ae,A\[Nu]} = Af;
  {mu,md,me,m\[Nu]D} = mf;
  
  retu = ret[[1]];
  ArrayQ[retu =
    ConjugateTranspose[Wu]
      .ComposeSfermion6x6MassMatrix[{
        {
          g2 MZ/Cos[\[Theta]W] gfL[\[Psi]u,\[Theta]W] Cos[\[Beta]] Id3,
          -Diag3[yu] \[Mu]/Sqrt[2]
        }, 
        {
          -Diag3[yu] Conjugate[\[Mu]]/Sqrt[2],
          g2 MZ/Cos[\[Theta]W] gfR[\[Psi]u,\[Theta]W] Cos[\[Beta]] Id3
        }
      }]
      .Wu];
  
  retd = ret[[2]];
  ArrayQ[retd =
    ConjugateTranspose[Wd]
      .ComposeSfermion6x6MassMatrix[{
        {
          g2 MZ/Cos[\[Theta]W] gfL[\[Psi]d, \[Theta]W] Cos[\[Beta]] Id3
          + Sqrt[2] Diag3[yd].Diag3[md],
          ConjugateTranspose[Ad]/Sqrt[2]
        },
        {
          Ad/Sqrt[2],
          g2 MZ/Cos[\[Theta]W] gfR[\[Psi]d, \[Theta]W] Cos[\[Beta]] Id3
          + Sqrt[2] Diag3[yd].Diag3[md]
        }
      }]
      .Wd];
  
  rete = ret[[3]];
  ArrayQ[rete =
    ConjugateTranspose[We]
      .ComposeSfermion6x6MassMatrix[{
        {
          g2 MZ/Cos[\[Theta]W] gfL[\[Psi]e, \[Theta]W] Cos[\[Beta]] Id3
          + Sqrt[2] Diag3[ye me], 
          ConjugateTranspose[Ae]/Sqrt[2]
        }, 
        {
          Ae/Sqrt[2],
          g2 MZ/Cos[\[Theta]W] gfR[\[Psi]e, \[Theta]W] Cos[\[Beta]] Id3
          + Sqrt[2] Diag3[ye me]
        }
      }]
      .We];
      
  ret\[Nu] =
    Pad\[Nu][
      g2 MZ/Cos[\[Theta]W] gfL[\[Psi]\[Nu], \[Theta]W] Cos[\[Beta]] Id3
    ];
  {retu,retd,rete,ret\[Nu]}
], {
  {gfL[__], _Complex}, {gfR[__], _Complex}, {Id3, _Complex, 2}, {Diag3[_], _Complex, 2},
  {ComposeSfermion6x6MassMatrix[__], _Complex, 2},
  {Pad\[Nu][__], _Complex, 2}, {ret, _Complex, 3},
  {Wu, _Complex, 2}, {Wd, _Complex, 2}, {We, _Complex, 2},
  {Au, _Complex, 2}, {Ad, _Complex, 2}, {Ae, _Complex, 2},
  {yu, _Complex, 1}, {yd, _Complex, 1}, {ye, _Complex, 1}
}];

\[Lambda]u = Compile[{
    {Wude, _Complex, 3}, {W\[Nu], _Complex, 2}, 
    {gp, _Real}, {g2, _Real}, {v, _Real},
    {yf, _Real, 2}, {Af, _Complex, 3}, {mf, _Real, 2},
    {\[Theta]W, _Real}, {\[Beta], _Real}, {\[Mu], _Complex}
  }, Module[{
    Wu,Wd,We,yu,yd,ye,y\[Nu],Au,Ad,Ae,A\[Nu],mu,md,me,m\[Nu]D,f,ret,Id3, retu,retd,rete,ret\[Nu],
    MZ = 1/2 v Sqrt[gp^2 + g2^2]
  },
  Id3 = {{1. + 0. I, 0. + 0. I, 0. + 0. I}, {0. + 0. I, 1. + 0. I, 0. + 0. I}, {0. + 0. I, 0. + 0. I, 1. + 0. I}};

  {Wu,Wd,We} = Wude;
  {yu,yd,ye,y\[Nu]} = yf;
  {Au,Ad,Ae,A\[Nu]} = Af;
  {mu,md,me,m\[Nu]D} = mf;
  ret = (1.0 + I 0.0) {Id3, Id3, Id3, Id3};
  
  retu = ret[[1]];
  ArrayQ[retu =
    ConjugateTranspose[ComplexMatrix[Wu]]
      .ComposeSfermion6x6MassMatrix[{
        {
          ComplexMatrix[-g2 MZ/Cos[\[Theta]W] gfL[\[Psi]u, \[Theta]W] Sin[\[Beta]] Id3
          + Sqrt[2] Diag3[yu].Diag3[mu]],
          ConjugateTranspose[Au]/Sqrt[2]
        },
        {
          Au/Sqrt[2],
          ComplexMatrix[-g2 MZ/Cos[\[Theta]W] gfR[\[Psi]u, \[Theta]W] Sin[\[Beta]] Id3
          + Sqrt[2] Diag3[yu].Diag3[mu]]
        }
      }]
      .ComplexMatrix[Wu]];

  retd = ret[[2]];
  ArrayQ[retd =
    ConjugateTranspose[ComplexMatrix[Wd]]
      .ComposeSfermion6x6MassMatrix[{
        {
          -g2 MZ/Cos[\[Theta]W] gfL[\[Psi]d, \[Theta]W] Sin[\[Beta]] Id3,
          -Diag3[yd] \[Mu]/Sqrt[2]
        }, 
        {
          -Diag3[yd] Conjugate[\[Mu]]/Sqrt[2], 
          -g2 MZ/Cos[\[Theta]W] gfR[\[Psi]d, \[Theta]W] Sin[\[Beta]] Id3 
        }
      }]
      .ComplexMatrix[Wd]];

  rete=ret[[3]];
  ArrayQ[rete =
    ConjugateTranspose[We]
      .ComposeSfermion6x6MassMatrix[{
        {
          -g2 MZ/Cos[\[Theta]W] gfL[\[Psi]e, \[Theta]W] Sin[\[Beta]] Id3, 
          -Diag3[ye] \[Mu]/Sqrt[2]
        },
        {
          -Diag3[ye] Conjugate[\[Mu]]/Sqrt[2], 
          -g2 MZ/Cos[\[Theta]W] gfR[\[Psi]e, \[Theta]W] Sin[\[Beta]] Id3
        }
      }]
      .We];

  ret\[Nu] =
    Pad\[Nu][
      -g2 MZ/Cos[\[Theta]W] gfL[\[Psi]\[Nu], \[Theta]W] Sin[\[Beta]] Id3
    ];
  {retu,retd,rete,ret\[Nu]}
], {
 {retu, _Complex, 2}, {retd, _Complex, 2}, {rete, _Complex, 2}, {ret\[Nu], _Complex, 2},
  {gfL[__], _Complex}, {gfR[__], _Complex}, {Id3, _Complex, 2}, {Diag3[_], _Complex, 2},
  {ComposeSfermion6x6MassMatrix[__], _Complex, 2},{ComplexMatrix[__], _Complex, 2},
  {Pad\[Nu][__], _Complex, 2}, {ret, _Complex, 3},
  {Wu, _Complex, 2}, {Wd, _Complex, 2}, {We, _Complex, 2},
  {Au, _Complex, 2}, {Ad, _Complex, 2}, {Ae, _Complex, 2},
  {yu, _Complex, 1}, {yd, _Complex, 1}, {ye, _Complex, 1}
}];

\[Lambda]As = Compile[{
    {Wude, _Complex, 3}, {W\[Nu], _Complex, 2}, 
    {yf, _Real, 2}, {Af, _Complex, 3},
    {\[Beta], _Real}, {\[Mu], _Complex}
  }, Module[{Wu,Wd,We,yu,yd,ye,y\[Nu],Au,Ad,Ae,A\[Nu],ret,Id3,Zero3, retu,retd,rete,ret\[Nu]},
  Id3 = {{1. + 0. I, 0. + 0. I, 0. + 0. I}, {0. + 0. I, 1. + 0. I, 0. + 0. I}, {0. + 0. I, 0. + 0. I, 1. + 0. I}};
  Zero3 = 0.0 Id3;
  {Wu,Wd,We} = Wude;
  {yu,yd,ye,y\[Nu]} = yf;
  {Au,Ad,Ae,A\[Nu]} = Af;
  ret = {Id3, Id3, Id3};

  retu =
    ConjugateTranspose[Wu]
      .ComposeSfermion6x6MassMatrix[{
        {
          Zero3, 
          -1/Sqrt[2] ( 
            ConjugateTranspose[Au] Cos[\[Beta]] 
            + \[Mu] Conjugate[DiagonalMatrix[yu]] Sin[\[Beta]]
          )
        }, {
          1/Sqrt[2] ( 
            Au Cos[\[Beta]] 
            + Conjugate[\[Mu]] Transpose[DiagonalMatrix[yu]] Sin[\[Beta]]
          ), 
          Zero3
        }
      }]
      .Wu;
      
  retd =
    ConjugateTranspose[Wd]
      .ComposeSfermion6x6MassMatrix[{
        {
          Zero3, 
          -1/Sqrt[2] (
            ConjugateTranspose[Ad] Sin[\[Beta]] 
            + \[Mu] Conjugate[DiagonalMatrix[yd]] Cos[\[Beta]]
          )
        }, 
        {
          1/Sqrt[2] (
            Ad Sin[\[Beta]] 
            + Conjugate[\[Mu]] Transpose[DiagonalMatrix[yd]] Cos[\[Beta]]
          ), 
          Zero3
        }
      }]
      .Wd;
      
  rete =
    ConjugateTranspose[We]
      .ComposeSfermion6x6MassMatrix[{
        {
          Zero3,
          -1/Sqrt[2] (
            ConjugateTranspose[Ae] Sin[\[Beta]] 
            + \[Mu] Conjugate[DiagonalMatrix[ye]] Cos[\[Beta]]
          )
        }, 
        {
          1/Sqrt[2] (
            Ae Sin[\[Beta]] 
            + Conjugate[\[Mu]] Transpose[DiagonalMatrix[ye]] Cos[\[Beta]]
          ), 
          Zero3
        }
      }]
      .We;
      
  ret\[Nu] =
    Pad\[Nu][
      Zero3
    ];
  
  {retu,retd,rete,ret\[Nu]}
]];

\[Lambda]AAs = Compile[{
    {Wude, _Complex, 3}, {W\[Nu], _Complex, 2}, 
    {yf, _Real, 2},
    {g2, _Real}, {\[Theta]W, _Real}, {\[Beta], _Real}
  }, Module[{Wu,Wd,We,yu,yd,ye,y\[Nu],ft,ret,Id3,Zero3, retu,retd,rete,ret\[Nu]},
  Id3 = {{1. + 0. I, 0. + 0. I, 0. + 0. I}, {0. + 0. I, 1. + 0. I, 0. + 0. I}, {0. + 0. I, 0. + 0. I, 1. + 0. I}};
  Zero3 = 0.0 Id3;
  {Wu,Wd,We} = Wude;
  {yu,yd,ye,y\[Nu]} = yf;
  ret =  {Id3, Id3, Id3, Id3};
  
  retu = ret[[1]];
  ArrayQ[retu =
    -2 ConjugateTranspose[Wu]
      .ComposeSfermion6x6MassMatrix[{
        {
        g2^2/4 Cos[2 \[Beta]] (1/Cos[\[Theta]W]^2 I3[\[Psi]u] - Tan[\[Theta]W]^2 Qe[\[Psi]u]) Id3 
        -1/2 Cos[\[Beta]]^2 DiagonalMatrix[yu]^2,
	Zero3
        },
	{
	Zero3,
	g2^2/4 Cos[2 \[Beta]] Tan[\[Theta]W]^2 Qe[\[Psi]u] Id3
	-1/2 Cos[\[Beta]]^2 DiagonalMatrix[yu]^2
	}
      }]
      .Wu];
  
  retd = ret[[2]];
  ArrayQ[retd =
    -2 ConjugateTranspose[Wd]
      .ComposeSfermion6x6MassMatrix[{
        {
        g2^2/4 Cos[2 \[Beta]] (1/Cos[\[Theta]W]^2 I3[\[Psi]d] - Tan[\[Theta]W]^2 Qe[\[Psi]d]) Id3 
        -1/2 Sin[\[Beta]]^2 DiagonalMatrix[yd]^2,
	Zero3
        },
	{
	Zero3,
	g2^2/4 Cos[2 \[Beta]] Tan[\[Theta]W]^2 Qe[\[Psi]d] Id3
	-1/2 Sin[\[Beta]]^2 DiagonalMatrix[yd]^2
	}
      }]
      .Wd];
  
  rete = ret[[3]];
  ArrayQ[rete =
    -2 ConjugateTranspose[We]
      .ComposeSfermion6x6MassMatrix[{
        {
        g2^2/4 Cos[2 \[Beta]] (1/Cos[\[Theta]W]^2 I3[\[Psi]e] - Tan[\[Theta]W]^2 Qe[\[Psi]e]) Id3 
        -1/2 Sin[\[Beta]]^2 DiagonalMatrix[ye]^2,
	Zero3
        },
	{
	Zero3,
	g2^2/4 Cos[2 \[Beta]] Tan[\[Theta]W]^2 Qe[\[Psi]e] Id3
	-1/2 Sin[\[Beta]]^2 DiagonalMatrix[ye]^2
	}
      }]
      .We];
  
  ret\[Nu] = ret[[4]];
  ArrayQ[ret\[Nu] =
   -2 Pad\[Nu][
    ConjugateTranspose[W\[Nu]]
      .(
        g2^2/4 Cos[2 \[Beta]](1/Cos[\[Theta]W]^2 I3[\[Psi]\[Nu]] - Tan[\[Theta]W]^2 Qe[\[Psi]\[Nu]]) Id3     - 1/2 Cos[\[Beta]]^2 DiagonalMatrix[y\[Nu]]^2
      )
      .W\[Nu]
    ]];
      
  {retu,retd,rete,ret\[Nu]}
]];

\[Lambda]Hps = Compile[{
    {Wude, _Complex, 3}, {W\[Nu], _Complex, 2}, 
    {yf, _Real, 2}, {Af, _Complex, 3},
    {g2, _Real}, {v, _Real}, {\[Beta], _Real}, {\[Mu], _Complex}
  }, Module[{Wu,Wd,We,yu,yd,ye,y\[Nu],Au,Ad,Ae,A\[Nu],ft,ret,Id3,Zero3,retu,ret\[Nu],MW},
  Id3 = {{1. + 0. I, 0. + 0. I, 0. + 0. I}, {0. + 0. I, 1. + 0. I, 0. + 0. I}, {0. + 0. I, 0. + 0. I, 1. + 0. I}};
  Zero3 = 0.0 Id3;
  {Wu,Wd,We} = Wude;
  {yu,yd,ye,y\[Nu]} = yf;
  {Au,Ad,Ae,A\[Nu]} = Af;
  ret =  {Id3, Id3,Id3, Id3};
  MW = g2 v/2;
  
  retu = ret[[1]];
  ArrayQ[retu =
    ConjugateTranspose[Wu]
      .ComposeSfermion6x6MassMatrix[{
        {
        g2/Sqrt[2] MW Sin[2 \[Beta]] Id3 - DiagonalMatrix[yu]^2 Cos[\[Beta]] v Sin[\[Beta]]/Sqrt[2] - DiagonalMatrix[yd]^2 Cos[\[Beta]] v Sin[\[Beta]]/Sqrt[2],
        -ConjugateTranspose[Ad] Sin[\[Beta]] - \[Mu] DiagonalMatrix[yd] Cos[\[Beta]]
        },
        {
        -Au Cos[\[Beta]] - Conjugate[\[Mu]] DiagonalMatrix[yu] Sin[\[Beta]],
        -DiagonalMatrix[yu].DiagonalMatrix[yd] v/Sqrt[2]
        }
      }]
      .Wd];

  
  ret\[Nu] = ret[[4]];
  ArrayQ[ret\[Nu] =
    Pad3x6\[Nu][
    ConjugateTranspose[W\[Nu]]
      .ComposeSfermion3x6MassMatrix[{
        {
        g2/Sqrt[2] MW Sin[2 \[Beta]] Id3 - DiagonalMatrix[ye]^2 Cos[\[Beta]] v Sin[\[Beta]]/Sqrt[2],
        - ConjugateTranspose[Ae] Sin[\[Beta]] - \[Mu] DiagonalMatrix[ye] Cos[\[Beta]]
        }
      }]
      .We
    ]];
   
  {retu,0 retu, 0 retu,ret\[Nu]}
]];

\[Lambda]HpHms = Compile[{
    {Wude, _Complex, 3}, {W\[Nu], _Complex, 2}, 
    {yf, _Real, 2},
    {g2, _Real}, {\[Theta]W, _Real}, {\[Beta], _Real}
  }, Module[{Wu,Wd,We,yu,yd,ye,y\[Nu],ft,ret,Id3,Zero3, retu,retd,rete,ret\[Nu]},
  Id3 = {{1. + 0. I, 0. + 0. I, 0. + 0. I}, {0. + 0. I, 1. + 0. I, 0. + 0. I}, {0. + 0. I, 0. + 0. I, 1. + 0. I}};
  Zero3 = 0.0 Id3;
  {Wu,Wd,We} = Wude;
  {yu,yd,ye,y\[Nu]} = yf;
  ret =  {Id3, Id3, Id3, Id3};
  
  retu = ret[[1]];
  ArrayQ[retu =
    -ConjugateTranspose[Wu]
      .ComposeSfermion6x6MassMatrix[{
        {
        -g2^2/2 Cos[2 \[Beta]] (Cos[2 \[Theta]W]/(Cos[\[Theta]W]^2) I3[\[Psi]u] + Tan[\[Theta]W]^2 Qe[\[Psi]u]) Id3 - Sin[\[Beta]]^2 DiagonalMatrix[yd]^2,
        Zero3
        },
        {
        Zero3,
        g2^2/2 Cos[2 \[Beta]] Tan[\[Theta]W]^2 Qe[\[Psi]u] Id3 - Cos[\[Beta]]^2 DiagonalMatrix[yu]^2
        }
      }]
      .Wu];

  retd = ret[[2]];
  ArrayQ[retd =
    -ConjugateTranspose[Wd]
      .ComposeSfermion6x6MassMatrix[{
        {
        -g2^2/2 Cos[2 \[Beta]] (Cos[2 \[Theta]W]/(Cos[\[Theta]W]^2) I3[\[Psi]d] + Tan[\[Theta]W]^2 Qe[\[Psi]d]) Id3 - Cos[\[Beta]]^2 DiagonalMatrix[yu]^2,
        Zero3
        },
        {
        Zero3,
        g2^2/2 Cos[2 \[Beta]] Tan[\[Theta]W]^2 Qe[\[Psi]d] Id3 - Sin[\[Beta]]^2 DiagonalMatrix[yd]^2
        }
      }]
      .Wd];

  rete = ret[[3]];
  ArrayQ[rete =
    -ConjugateTranspose[We]
      .ComposeSfermion6x6MassMatrix[{
        {
        -g2^2/2 Cos[2 \[Beta]] (Cos[2 \[Theta]W]/(Cos[\[Theta]W]^2) I3[\[Psi]e] + Tan[\[Theta]W]^2 Qe[\[Psi]e]) Id3,
        Zero3
        },
        {
        Zero3,
        g2^2/2 Cos[2 \[Beta]] Tan[\[Theta]W]^2 Qe[\[Psi]e] Id3 - Sin[\[Beta]]^2 DiagonalMatrix[ye]^2
        }
      }]
      .We];
  
  
  ret\[Nu] = ret[[4]];
  ArrayQ[ret\[Nu] =
    Pad\[Nu][
    -ConjugateTranspose[W\[Nu]]
      .(
        -g2^2/2 Cos[2 \[Beta]] (Cos[2 \[Theta]W]/(Cos[\[Theta]W]^2) I3[\[Psi]\[Nu]] + Tan[\[Theta]W]^2 Qe[\[Psi]\[Nu]]) Id3 - Sin[\[Beta]]^2 DiagonalMatrix[ye]^2
      )
      .W\[Nu]
    ]];
  
  {retu,retd,rete,ret\[Nu]}
]];

(* Z and A self-energies *)

\[CapitalPi]ZZT = Compile[{
    {p, _Real}, {\[Theta]W, _Real}, {gp, _Real}, {g2, _Real}, {v, _Real}, {\[Alpha], _Real}, {\[Beta], _Real}, 
    {MZ, _Real}, {MW, _Real}, {Mh, _Real}, {MH, _Real}, {MA, _Real}, {MHp, _Real}, 
    {mf, _Real, 2}, {Wude, _Complex, 3}, {W\[Nu], _Complex, 2}, {msude, _Real, 2}, {ms\[Nu], _Real, 1},
    {N, _Complex, 2}, {m\[Chi]0, _Real, 1}, {U, _Complex, 2}, {V, _Complex, 2}, {m\[Chi]p, _Real, 1}, {Q, _Real}
  }, Module[{ret},
  ret = 0.0 + I 0.0;
  
  ArrayQ[ret = Loop g2^2/Cos[\[Theta]W]^2 * With[{
      sab = Sin[\[Alpha] - \[Beta]],
      cab = Cos[\[Alpha] - \[Beta]]
    },
    - sab^2 (
      B22t[p, MA, MH, Q] + B22t[p, MZ, Mh, Q] - MZ^2 B0[p, MZ, Mh, Q]
    )
    - cab^2 (
      B22t[p, MZ, MH, Q] + B22t[p, MA, Mh, Q] - MZ^2 B0[p, MZ, MH, Q]
    )
    - 2 Cos[\[Theta]W]^4 (
      2 p^2 + MW^2 - MZ^2 Sin[\[Theta]W]^4/Cos[\[Theta]W]^2
    ) B0[p, MW, MW, Q]
    - (8 Cos[\[Theta]W]^4 + Cos[2 \[Theta]W]^2) B22t[p, MW, MW, Q]
    - Cos[2 \[Theta]W]^2 B22t[p, MHp, MHp, Q]
    - Sum[
      With[
        {
          g = If[ft == \[Psi]\[Nu], 3, 6],
          Wf = {Wude[[1]], Wude[[2]], Wude[[3]], W\[Nu]}[[ft]],
          msf = {msude[[1]], msude[[2]], msude[[3]], ms\[Nu]}[[ft]]
        },
        NC[ft] *
        Sum[
          (
            Sum[I3[ft] Conjugate[Wf[[i, s]]] Wf[[i, t]], {i, 3}] 
            - Qe[ft] Sin[\[Theta]W]^2 KroneckerDelta[s, t]
          ) * (
            Sum[I3[ft] Conjugate[Wf[[j, t]]] Wf[[j, s]], {j, 3}]
            - Qe[ft] Sin[\[Theta]W]^2 KroneckerDelta[s, t]
          ) 4 B22[p, msf[[s]], msf[[t]], Q]
          - 1/2 (
            Sum[
              (1 - 8 I3[ft] Qe[ft] Sin[\[Theta]W]^2) Conjugate[Wf[[i, s]]] Wf[[i, t]], 
              {i, 3}
            ]
            + 4 Qe[ft]^2 Sin[\[Theta]W]^4 KroneckerDelta[s, t]
          ) KroneckerDelta[s, t] A0[msf[[s]], Q], {s, g}, {t, g}
        ]
      ],
      {ft, 4}
    ]
    + Sum[
      NC[ft] * (
        (
          gfL[ft, \[Theta]W]^2
          + gfR[ft, \[Theta]W]^2
        ) H[p, mf[[ft, j]], mf[[ft, j]], Q] 
        - 4 gfL[ft, \[Theta]W] gfR[ft, \[Theta]W] mf[[ft, j]]^2
          B0[p, mf[[ft, j]], mf[[ft, j]], Q]
      ), 
      {ft, 4},
      {j, 3}
    ]
    + Cos[\[Theta]W]^2/(2 g2^2) Sum[
      fZ0[i, j, N, g2, \[Theta]W] H[p, m\[Chi]0[[i]], m\[Chi]0[[j]], Q] 
      + 2 gZ0[i, j, N, g2, \[Theta]W] m\[Chi]0[[i]] m\[Chi]0[[j]] B0[p, m\[Chi]0[[i]], m\[Chi]0[[j]], Q], 
      {i, 4}, {j, 4}
    ]
    + Cos[\[Theta]W]^2/g2^2 Sum[
      fZp[i, j, U, V, g2, \[Theta]W] H[p, m\[Chi]p[[i]], m\[Chi]p[[j]], Q] 
      + 2 gZp[i, j, U, V, g2, \[Theta]W] m\[Chi]p[[i]] m\[Chi]p[[j]] B0[p, m\[Chi]p[[i]], m\[Chi]p[[j]], Q], 
      {i, 2}, {j, 2}
    ]
  ]];
  ret
], {{ret, _Complex}, {fZ0[__], _Complex}, {gZ0[__], _Complex}, {fZp[__], _Complex}, {gZp[__], _Complex}}];

      
\[CapitalPi]AA = Compile[{
    {p, _Real}, {\[Theta]W, _Real}, {gp, _Real}, {g2, _Real}, {v, _Real}, {\[Alpha], _Real}, {\[Beta], _Real},
    {MZ, _Real}, {MW, _Real}, {Mh, _Real}, {MH, _Real}, {\[Mu], _Complex}, {MA, _Real}, {MHp, _Real},
    {y, _Real, 2}, {A, _Complex, 3}, {mf, _Real, 2}, 
    {Wude, _Complex, 3}, {W\[Nu], _Complex, 2}, {msude, _Real, 2}, {ms\[Nu], _Real, 1},
    {N, _Complex, 2}, {m\[Chi]0, _Real, 1}, {U, _Complex, 2}, {V, _Complex, 2}, {m\[Chi]p, _Real, 1}, {Q, _Real}
  }, Module[{\[Lambda]AAs$, \[Lambda]As$, ret},
  \[Lambda]As$ = \[Lambda]As[Wude, W\[Nu], y, A, \[Beta], \[Mu]];
  \[Lambda]AAs$ = \[Lambda]AAs[Wude, W\[Nu], y, g2, \[Theta]W, \[Beta]];
  ret = 0.0 + 0.0 I;
  
  ArrayQ[ret = Loop With[{
      sab = Sin[\[Alpha] - \[Beta]],
      cab = Cos[\[Alpha] - \[Beta]]
    },
    Cos[\[Beta]]^2 Sum[
      NC[ft] y[[ft, j]]^2 (
        p^2 B0[p, mf[[ft, j]], mf[[ft, j]], Q] 
        - 2 A0[mf[[ft, j]], Q]
      ), 
      {ft, {1 (*==\[Psi]u*), 4 (*==\[Psi]\[Nu]*)}}, {j, 3}
    ]
    + Sin[\[Beta]]^2 Sum[
      NC[ft] y[[ft, j]]^2 (
        p^2 B0[p, mf[[ft, j]], mf[[ft, j]], Q] 
        - 2 A0[mf[[ft, j]], Q]
      ), 
      {ft, {2 (*==\[Psi]d*), 3 (*==\[Psi]e*)}}, {j, 3}
    ]
    + Sum[
      With[{
          g = If[ft == \[Psi]\[Nu], 3, 6],
          msf = {msude[[1]], msude[[2]], msude[[3]], ms\[Nu]}[[ft]]
        },
        Sum[
          NC[ft] \[Lambda]AAs$[[ft, i, i]] A0[msf[[i]], Q],
          {i, g}
        ]
      ], 
      {ft, 4}
    ]
    + Sum[
      With[{
          g = If[ft == \[Psi]\[Nu], 3, 6],
          msf = {msude[[1]], msude[[2]], msude[[3]], ms\[Nu]}[[ft]]
        },
        Sum[
          NC[ft] \[Lambda]As$[[ft, i, j]]^2 B0[p, msf[[i]], msf[[j]], Q],
          {i, g}, {j, g}
        ]
      ],
      {ft, 4}
    ]
    + g2^2/4 (
      2 F[p, MHp, MW, Q] 
      + sab^2/Cos[\[Theta]W]^2 F[p, MH, MZ, Q] 
      + cab^2/Cos[\[Theta]W]^2 F[p, Mh, MZ, Q]
    )
    + 1/2 Sum[
      Sum[
        \[Lambda]A[n, m, gp, g2, v, \[Theta]W, \[Alpha], \[Beta]]^2 
          B0[
            p,
            (* order: \[Phi]Gp, \[Phi]Hp, \[Phi]h, \[Phi]H, \[Phi]G0, \[Phi]A *)
            {0, 0, Mh, MH, MZ, MA}[[n]],
            {0, 0, Mh, MH, MZ, MA}[[m]],
            Q
          ], 
        {m, {3,4,5,6} (* Higgs0 *)}
      ] 
      + \[Lambda]AA[n, g2, \[Theta]W, \[Alpha], \[Beta]] 
        * A0[{0, 0, Mh, MH, MZ, MA}[[n]], Q],
      {n, {3,4,5,6} (* Higgs0 *)}
    ]
    + g2^2 MW^2/2 B0[p, MW, MHp, Q]
    + Sum[
      \[Lambda]AA[x, g2, \[Theta]W, \[Alpha], \[Beta]]
      A0[{MW, MHp}[[x]], Q],
      {x, {1, 2} (* Higgsp *)}
    ]
    + g2^2 (2 A0[MW, Q] + A0[MZ, Q]/Cos[\[Theta]W]^2)
    + 1/2 Sum[
      fA0[i, j, N, gp, g2, \[Beta]] G[p, m\[Chi]0[[i]], m\[Chi]0[[j]], Q] 
      - 2 gA0[i, j, N, gp, g2, \[Beta]] m\[Chi]0[[i]] m\[Chi]0[[j]] B0[p, m\[Chi]0[[i]], m\[Chi]0[[j]], Q], 
      {i, 4}, {j, 4}
    ]
    + Sum[
      fAp[i, j, U, V, g2, \[Beta]] G[p, m\[Chi]p[[i]], m\[Chi]p[[j]], Q] 
      - 2 gAp[i, j, U, V, g2, \[Beta]] m\[Chi]p[[i]] m\[Chi]p[[j]] B0[p, m\[Chi]p[[i]], m\[Chi]p[[j]], Q], 
      {i, 2}, {j, 2}
    ]
  ]];
  ret
],
{{ret, _Complex}, {\[Lambda]AAs$, _Complex}, {\[Lambda]As$, _Complex},{fA0[__], _Complex}, {gA0[__], _Complex}, {fAp[__], _Complex}, {gAp[__], _Complex}}
];

\[CapitalPi]HpHm = Compile[{
    {p, _Real}, {\[Theta]W, _Real}, {gp, _Real}, {g2, _Real}, {v, _Real}, {\[Alpha], _Real}, {\[Beta], _Real},
    {MZ, _Real}, {MW, _Real}, {Mh, _Real}, {MH, _Real}, {\[Mu], _Complex}, {MA, _Real}, {MHp, _Real},
    {y, _Real, 2}, {A, _Complex, 3}, {mf, _Real, 2}, 
    {Wude, _Complex, 3}, {W\[Nu], _Complex, 2}, {msude, _Real, 2}, {ms\[Nu], _Real, 1},
    {N, _Complex, 2}, {m\[Chi]0, _Real, 1}, {U, _Complex, 2}, {V, _Complex, 2}, {m\[Chi]p, _Real, 1}, {Q, _Real}
  }, Module[{\[Lambda]HpHms$, \[Lambda]Hps$, ret},
  \[Lambda]Hps$ = \[Lambda]Hps[Wude, W\[Nu], y, A, g2, v, \[Beta], \[Mu]];
  \[Lambda]HpHms$ = \[Lambda]HpHms[Wude, W\[Nu], y, g2, \[Theta]W, \[Beta]]; 
  ret = 0.0 + 0.0 I;
  
  ArrayQ[ret = Loop With[{
      sab = Sin[\[Alpha] - \[Beta]],
      cab = Cos[\[Alpha] - \[Beta]],
      e = g2 Sin[\[Theta]W]
    },
    Sum[
      NC[ft] (Cos[\[Beta]]^2 y[[ft, j]]^2 + Sin[\[Beta]]^2 y[[Floor[(ft + 3)/2], j]]) G[p, mf[[ft, j]], mf[[Floor[(ft + 3)/2], j]], Q]
        -NC[ft] 2 Sin[2 \[Beta]] y[[ft, j]] y[[Floor[(ft + 3)/2], j]] mf[[ft, j]] mf[[Floor[(ft + 3)/2], j]] B0[p, mf[[ft, j]], mf[[Floor[(ft + 3)/2], j]], Q],
      {ft, {1 (*==\[Psi]u*), 4 (*==\[Psi]\[Nu]*)}}, {j, 3}   (*Floor[(1 + 3)/2] == \[Psi]d, Floor[(4 + 3)/2] == \[Psi]e*)
    ]
    + Sum[
      With[{
          g = If[ft == \[Psi]\[Nu], 3, 6],
          msf = {msude[[1]], msude[[2]], msude[[3]], ms\[Nu]}[[ft]],
          msfSU2partner = {msude[[1]], msude[[2]], msude[[3]], ms\[Nu]}[[Floor[(ft + 3)/2]]]
        },
        Sum[
          NC[ft] \[Lambda]Hps$[[ft, i, j]]^2 B0[p, msf[[i]], msfSU2partner[[j]], Q],
          {i, g}, {j, 6}
        ]
      ],
      {ft, {1 (*==\[Psi]u*), 4 (*==\[Psi]\[Nu]*)}}
    ]
    + Sum[
      With[{
          g = If[ft == \[Psi]\[Nu], 3, 6],
          msf = {msude[[1]], msude[[2]], msude[[3]], ms\[Nu]}[[ft]]
        },
        Sum[
          NC[ft] \[Lambda]HpHms$[[ft, i, i]] A0[msf[[i]], Q],
          {i,g}
        ]
      ],
      {ft, 4}
    ]
    + g2^2/4 (
      sab^2 F[p, MH, MW, Q]
      + cab^2 F[p, Mh, MW, Q] 
      + F[p, MA, MW, Q] 
      + Cos[2 \[Theta]W]^2/Cos[\[Theta]W]^2 F[p, MHp, MZ, Q]
    )
    + e^2 F[p, MHp, 0, Q]
    + 2 g2^2 A0[MW, Q]
    + g2^2 Cos[2 \[Theta]W]^2/Cos[\[Theta]W]^2 A0[MZ, Q]
    + Sum[
        Sum[
          \[Lambda]HpH0Hm[n, m, gp, g2, v, \[Theta]W, \[Alpha], \[Beta]]^2 
          B0[
            p,
            (* order: \[Phi]Gp, \[Phi]Hp, \[Phi]h, \[Phi]H, \[Phi]G0, \[Phi]A *)
            {MW, MHp, Mh, MH, 0, 0}[[n]],
            {MW, MHp, Mh, MH, 0, 0}[[m]],
            Q
          ],
          {n, {3, 4} (* Higgs0p *)}
        ]
        + \[Lambda]HpHmHpHm[m, g2, \[Theta]W, \[Beta]] A0[{MW, MHp, Mh, MH, 0, 0}[[m]], Q], 
        {m, {1, 2} (* Higgsp *)}
    ]
    + g2^2 MW^2/4 B0[p, MW, MA, Q]
    + 1/2 Sum[
      \[Lambda]HpHmH0H0[m, g2, \[Theta]W, \[Alpha], \[Beta]] A0[{0, 0, Mh, MH, MZ, MA}[[m]], Q],
      {m, {3,4,5,6} (* Higgs0 *)}
    ]
    +Sum[
      fHp[i, j, N, U, V, gp, g2, \[Beta]] G[p, m\[Chi]p[[j]], m\[Chi]0[[i]], Q],
      {j, 2}, {i, 4}
    ] 
    -2 Sum[m\[Chi]p[[j]] m\[Chi]0[[i]] gHp[i, j, N, U, V, gp, g2, \[Beta]] B0[p, m\[Chi]p[[j]], m\[Chi]0[[i]], Q],
      {j, 2}, {i, 4}
    ]
  ]];
  ret
]];

(* Higgs tadpoles *)
tadtd = Compile[{
    {\[Theta]W, _Real}, {gp, _Real}, {g2, _Real}, {v, _Real}, {\[Alpha], _Real}, {\[Beta], _Real},
    {MZ, _Real}, {MW, _Real}, {\[Mu], _Complex}, {Mh, _Real}, {MH, _Real}, {MA, _Real}, {MHp, _Real}, 
    {y, _Real, 2}, {A, _Complex, 3}, {mf, _Real, 2}, 
    {Wude, _Complex, 3}, {W\[Nu], _Complex, 2}, {msude, _Real, 2}, {ms\[Nu], _Real, 1},
    {N, _Complex, 2}, {m\[Chi]0, _Real, 1}, {U, _Complex, 2}, {V, _Complex, 2}, {m\[Chi]p, _Real, 1}, {Q, _Real}
  }, Module[{\[Lambda]d$, ret, Id3},
  Id3 = {{1. + 0. I, 0. + 0. I, 0. + 0. I}, {0. + 0. I, 1. + 0. I, 0. + 0. I}, {0. + 0. I, 0. + 0. I, 1. + 0. I}};
  \[Lambda]d$ = {Id3,Id3} (1.0 + I 0.0);
  \[Lambda]d$ = \[Lambda]d[Wude, W\[Nu], gp, g2, v, y, A, mf, \[Theta]W, \[Beta], \[Mu]];
  ret = 1.0 + I 0.0;
  ret = Loop * (
    - Sum[
      2 NC[ft] y[[ft, j]]^2 A0[mf[[ft, j]], Q],
      {ft, {2 (*==\[Psi]d*), 3 (*==\[Psi]e*)}}, 
      {j, 3}
    ]
    + Sum[With[{
          g = If[ft == \[Psi]\[Nu], 3, 6],
          msf = {msude[[1]], msude[[2]], msude[[3]], ms\[Nu]}[[ft]]
        },
        Complexify[Sum[
          NC[ft] g2^2/(2 MW Cos[\[Beta]])
          * \[Lambda]d$[[ft, i, i]]
          * A0[msf[[i]], Q],
          {i, g}
        ]]
      ],
      {ft, 4}
    ]
    - g2^2 Cos[2 \[Beta]]/(8 Cos[\[Theta]W]^2) (A0[MA, Q] + 2 A0[MHp, Q])
    + g2^2/2 A0[MHp, Q]
    + g2^2/(8 Cos[\[Theta]W]^2) (
      3 Sin[\[Alpha]]^2 - Cos[\[Alpha]]^2 
      + Sin[2 \[Alpha]] Tan[\[Beta]]
    ) A0[Mh, Q]
    + g2^2/(8 Cos[\[Theta]W]^2) (
      3 Cos[\[Alpha]]^2 - Sin[\[Alpha]]^2 
      - Sin[2 \[Alpha]] Tan[\[Beta]]
    ) A0[MH, Q]
    - Sum[
      g2^2 m\[Chi]0[[i]]/(MW Cos[\[Beta]]) Re[N[[i, 3]] (N[[i, 2]] - N[[i, 1]] Tan[\[Theta]W])]
      * A0[m\[Chi]0[[i]], Q],
      {i, 4}
    ]
    - Sum[
      Sqrt[2] g2^2 m\[Chi]p[[i]]/(MW Cos[\[Beta]]) Re[V[[i, 1]] U[[i, 2]]] A0[m\[Chi]p[[i]], Q], 
      {i, 2}
    ]
    + 3 g2^2/4 (2 A0[MW, Q] + A0[MZ, Q]/Cos[\[Theta]W]^2)
    + g2^2 Cos[2 \[Beta]]/(8 Cos[\[Theta]W]^2) (2 A0[MW, Q] + A0[MZ, Q])
  )
],{{ret, _Complex}, {\[Lambda]d$, _Complex}}
];

tadtu = Compile[{
    {\[Theta]W, _Real}, {gp, _Real}, {g2, _Real}, {v, _Real}, {\[Alpha], _Real}, {\[Beta], _Real}, 
    {MZ, _Complex}, {MW, _Complex}, {\[Mu], _Complex}, {Mh, _Real}, {MH, _Real}, {MA, _Real}, {MHp, _Real}, 
    {y, _Real, 2}, {A, _Complex, 3}, {mf, _Real, 2}, 
    {Wude, _Complex, 3}, {W\[Nu], _Complex, 2}, {msude, _Real, 2}, {ms\[Nu], _Real, 1},
    {N, _Complex, 2}, {m\[Chi]0, _Real, 1}, {U, _Complex, 2}, {V, _Complex, 2}, {m\[Chi]p, _Real, 1}, {Q, _Real}
  }, Module[{\[Lambda]u$, ret, Id3},
  Id3 = {{1. + 0. I, 0. + 0. I, 0. + 0. I}, {0. + 0. I, 1. + 0. I, 0. + 0. I}, {0. + 0. I, 0. + 0. I, 1. + 0. I}};
  \[Lambda]u$ = {Id3,Id3} (1.0 + I 0.0);
  \[Lambda]u$ = \[Lambda]u[Wude, W\[Nu], gp, g2, v, y, A, mf, \[Theta]W, \[Beta], \[Mu]];
  ret = 1.0 + I 0.0;
  
  ret = Loop * (
    - Sum[
      2. NC[ft] y[[ft, j]]^2 A0[mf[[ft, j]], Q], 
      {ft, {1 (*==\[Psi]u*), 4 (*==\[Psi]\[Nu]*)}}, 
      {j, 3}
    ]
    + Sum[With[{
          g = If[ft == \[Psi]\[Nu], 3, 6],
          msf = {msude[[1]], msude[[2]], msude[[3]], ms\[Nu]}[[ft]]
        },
        Complexify[Sum[
          NC[ft] g2^2/(2 MW Sin[\[Beta]])
          * \[Lambda]u$[[ft, i, i]]
          * A0[msf[[i]], Q],
          {i, g}
        ]]
      ],
      {ft, 4}
    ]
    + g2^2 Cos[2 \[Beta]]/(8 Cos[\[Theta]W]^2) (A0[MA, Q] + 2 A0[MHp, Q])
    + g2^2/2 A0[MHp, Q]
    + g2^2/(8 Cos[\[Theta]W]^2) (
      3 Cos[\[Alpha]]^2 - Sin[\[Alpha]]^2 + Sin[2 \[Alpha]] Cot[\[Beta]]
    ) A0[Mh, Q]
    + g2^2/(8 Cos[\[Theta]W]^2) (
      3 Sin[\[Alpha]]^2 - Cos[\[Alpha]]^2 - Sin[2 \[Alpha]] Cot[\[Beta]]
    ) A0[MH, Q]
    + Sum[
      g2^2 m\[Chi]0[[i]]/(MW Sin[\[Beta]]) Re[N[[i, 4]] (N[[i, 2]] - N[[i, 1]] Tan[\[Theta]W])]
      * A0[m\[Chi]0[[i]], Q], 
      {i, 4}
    ]
    - Sum[
      Sqrt[2] g2^2 m\[Chi]p[[i]]/(MW Sin[\[Beta]]) Re[V[[i, 2]] U[[i, 1]]] A0[m\[Chi]p[[i]], Q], 
      {i, 2}
    ]
    + 3 g2^2/4 (2 A0[MW, Q] + A0[MZ, Q]/Cos[\[Theta]W]^2)
    - g2^2 Cos[2 \[Beta]]/(8 Cos[\[Theta]W]^2) (2 A0[MW, Q] + A0[MZ, Q])
    
  ) 
],{{ret, _Complex}, {\[Lambda]u$, _Complex, 3}, {\[Lambda]u, _Complex, 3}}
];


STCGetQuarkMixing[Yu_,Yd_]:=Block[{UuL1,UuR1,yu,Yd2,UCKM,UR,yd,t12,t13,t23,d,d1,d2,d3,p1,p2,UdL1,UdR1,UuL,UuR,UdL,UdR},
   {UuL1, yu, UuR1} = SingularValueDecomposition[Conjugate[Yu]];
     yu=Reverse@Diagonal@yu;
       UuL1=Transpose@Reverse@Transpose@UuL1;
       UuR1=Transpose@Reverse@Transpose@UuR1;
     Yd2=ConjugateTranspose[UuL1].Conjugate[Yd].UuR1;
   {UCKM, yd, UR} = SingularValueDecomposition[Yd2];
     yd=Reverse@Diagonal@yd;
       UCKM=Transpose@Reverse@Transpose@UCKM;
       UR=Transpose@Reverse@Transpose@UR;
	UdL1=UuL1.UCKM;
	UdR1=UuR1.UR;
	{t12,t13,t23,d,d1,d2,d3,p1,p2}=MPT3x3MixingParameters[UCKM];
	UuL=UuL1.DiagonalMatrix[{E^(I d1),E^(I d2),E^(I d3)}];
	UuR=UuR1.DiagonalMatrix[{E^(I d1),E^(I d2),E^(I d3)}];
	UdL=UdL1.DiagonalMatrix[{E^(I/2 p1),E^(I/2 p2),1}];
	UdR=UdR1.DiagonalMatrix[{E^(I/2 p1),E^(I/2 p2),1}];
	{UuL,yu,UuR,UdL,yd,UdR}
];

STCGetLeptonMixing[Mn_,Ye_]:=Block[{UeL1,UeR1,ye,Mn2,UPMNS,Un,mn,t12,t13,t23,d,d1,d2,d3,p1,p2,UeL,UeR},
   {UeL1, ye, UeR1} = SingularValueDecomposition[Conjugate[Ye]];
     ye=Reverse@Diagonal@ye;
       UeL1=Transpose@Reverse@Transpose@UeL1;
       UeR1=Transpose@Reverse@Transpose@UeR1;
     Mn2=Transpose[UeL1].Mn.UeL1;
     UPMNS = MPT3x3NeutrinoMixingMatrix[Mn2];
	Un=UeL1.UPMNS;
	mn=Diagonal@(Transpose[UPMNS].Mn2.UPMNS);
	{t12,t13,t23,d,d1,d2,d3,p1,p2}=MPT3x3MixingParameters[UPMNS];
	UeL=UeL1.DiagonalMatrix[{E^(I d1),E^(I d2),E^(I d3)}];
	UeR=UeR1.DiagonalMatrix[{E^(I d1),E^(I d2),E^(I d3)}];
	{Un,mn,UeL,ye,UeR}
];


SCKMRotate[x_] := Block[{
    x$, Y$ = {0,0,0,0}, T$ = {0,0,0,0},
    UuL, UuR, UdL, UdR, UeL, UeR, U\[Nu]L
  },
  x$ = Join[
    x,
    {
      (* pole masses *) 
      "MZ" -> MZ,
      "MW" -> MW
    }
  ];
  
  (* Beware of relative minus signs and additional transpositions due to convention translation REAP -> SLHA *)
  {UuL, Y$[[1]], UuR, UdL, Y$[[2]], UdR} = STCGetQuarkMixing[Transpose["Yu" /. x$], Transpose["Yd" /. x$]];
  {U\[Nu]L, Y$[[4]], UeL, Y$[[3]], UeR } = STCGetLeptonMixing["M\[Nu]" /. x$, Transpose["Ye" /. x$]];
  Y$[[4]] = {0,0,0};
  
	(*SLHA SCKM defintion  Beware of relative minus signs and additional transpositions due to convention translation REAP -> SLHA *)

  T$[[1]] = ConjugateTranspose[UuR].("Tu" /. x$).UuL;
  T$[[2]] = ConjugateTranspose[UdR].("Td" /. x$).UdL;
  T$[[3]] = ConjugateTranspose[UeR].("Te" /. x$).UeL;
  T$[[4]] = Zero3;
  
  Join[
    FilterRules[
      x$,
      Except["Yu" | "Yd" | "Ye" | "Tu" | "Td" | "Te" | "m2Q" | "m2L" | "m2u" | "m2d" | "m2e" | "M\[Nu]"]
    ],
    {
      (* unrotated matrices *)
      FlavourBasis["Yu"] -> ("Yu" /. x$),
      FlavourBasis["Yd"] -> ("Yd" /. x$),
      FlavourBasis["Ye"] -> ("Ye" /. x$),
      
      FlavourBasis["Tu"] -> ("Tu" /. x$),
      FlavourBasis["Td"] -> ("Td" /. x$),
      FlavourBasis["Te"] -> ("Te" /. x$),
      
      FlavourBasis["m2Q"] -> ("m2Q" /. x$),
      FlavourBasis["m2L"] -> ("m2L" /. x$),
      FlavourBasis["m2u"] -> ("m2u" /. x$),
      FlavourBasis["m2d"] -> ("m2d" /. x$),
      FlavourBasis["m2e"] -> ("m2e" /. x$),

      "VCKM"->ConjugateTranspose[UuL].UdL,
      "VPMNS"->ConjugateTranspose[UeL].U\[Nu]L,

      (* SCKM rotated matrices *)
      SCKMBasis["m2Q"]      -> ConjugateTranspose[UdL].("m2Q" /. x$).UdL,
      SCKMBasis["m2L"]      -> ConjugateTranspose[UeL].("m2L" /. x$).UeL,
      
      SCKMBasis["m2u"]       -> ConjugateTranspose[UuR].("m2u" /. x$).UuR,
      SCKMBasis["m2d"]       -> ConjugateTranspose[UdR].("m2d" /. x$).UdR,
      SCKMBasis["m2e"]       -> ConjugateTranspose[UeR].("m2e" /. x$).UeR,
      
      SCKMBasis["Y"]         -> MustBeReal[Y$, "SCKMBasis[Y]"],
      SCKMBasis["T"]         -> T$
    }
  ]
];

TreeLevelVEV[MZ_, gp_, g2_] := 2 MZ / Sqrt[gp^2 + g2^2];

TreeLevel\[Mu][MZ_, m2Hd_, m2Hu_, \[Beta]_, signMu_] := signMu * MustBeReal[
  Sqrt[
    1/2 (Tan[2 \[Beta]] (m2Hu Tan[\[Beta]] - m2Hd Cot[\[Beta]]) - MZ^2)
  ],
  "|\[Mu]|"
];

TreeLevelm3[MZ_, m2Hd_, m2Hu_, \[Beta]_, signMu_] := MustBeReal[
  Sqrt[
    1/2 (Tan[2 \[Beta]](m2Hu - m2Hd) - MZ^2 Sin[2 \[Beta]])
  ],
  "|m3|"
];

TreeLevelMHp[MZ_, MW_, m2Hd_, m2Hu_, \[Beta]_, signMu_] := MustBeReal[
  Sqrt[1/Cos[2 \[Beta]] (m2Hu - m2Hd) - MZ^2 + MW^2],
  "MHp"
];

TreeLevelMA[MZ_, m2Hd_, m2Hu_, \[Beta]_, signMu_] := MustBeReal[
  Sqrt[1/Cos[2 \[Beta]] (m2Hu - m2Hd) - MZ^2],
  "MA"
];

CombineRules[dominant_, recessive_] := Join[dominant, FilterRules[recessive, Except[dominant]]]

CalcTreeValues[x_] := Block[{
    x$,
    Wude = {0,0,0},
    W\[Nu],
    msude = {0,0,0},
    ms\[Nu],
    mf = {0,0,0,0},
    M2\[Nu],
    M2ude = {0,0,0}
  },
  (* include simple Higgs parameters *)

x$ = CombineRules[
    x,
    {
      "vev" -> MustBeReal[TreeLevelVEV["MZ" /. x, "gp" /. x, "g2" /. x], "vev"],
      "\[Theta]W" -> MustBeReal[ArcSin[("gp" /. x) / Sqrt[("gp" /. x)^2 + ("g2" /. x)^2]], "\[Theta]W"],
      "\[Mu]" -> TreeLevel\[Mu]["MZ" /. x, "m2Hd" /. x, "m2Hu" /. x, ArcTan["tan\[Beta]" /. x], "sign\[Mu]" /. x],
      "m3" -> TreeLevelm3["MZ" /. x, "m2Hd" /. x, "m2Hu" /. x, ArcTan["tan\[Beta]" /. x], "sign\[Mu]" /. x]
    }
  ];


  x$ = CombineRules[
    x$,
    {
      If[      
         ("cMSSM" /. x) == True,
         "MHp" -> TreeLevelMHp["MZ" /. x$, ("g2" /. x) ("vev" /. x$)/2, "m2Hd" /. x$, "m2Hu" /. x$, ArcTan["tan\[Beta]" /. x$], "sign\[Mu]" /. x$],
         (*else "cMSSM" == False*)
         "MA" -> TreeLevelMA["MZ" /. x$, "m2Hd" /. x$, "m2Hu" /. x$, ArcTan["tan\[Beta]" /. x$], "sign\[Mu]" /. x$]
      ]
    }
  ];
  
  mf[[1]] = 1/Sqrt[2] ("vev" /. x$) Sin[ArcTan["tan\[Beta]" /. x$]] (SCKMBasis["Y"] /. x$)[[1]];
  mf[[2]] = 1/Sqrt[2] ("vev" /. x$) Cos[ArcTan["tan\[Beta]" /. x$]] (SCKMBasis["Y"] /. x$)[[2]];
  mf[[3]] = 1/Sqrt[2] ("vev" /. x$) Cos[ArcTan["tan\[Beta]" /. x$]] (SCKMBasis["Y"] /. x$)[[3]];
  mf[[4]] = {0,0,0};
  
  x$ = CombineRules[
    {
      If[      
        ("cMSSM" /. x) == True,
        "MA" -> MustBeReal[Sqrt[("MHp" /. x$)^2 - ("g2" /. x$)^2 ("vev" /. x$)^2 /4], "MA"],
        (*else cMSSM == False*)
        "MHp" -> MustBeReal[Sqrt[("MA" /. x$)^2 + ("g2" /. x$)^2 ("vev" /. x$)^2 /4], "MHp"]
      ],
      "mf" -> MustBeReal[mf, "mf"]
    },
    x$
  ];
    
  (* sfermion mass eigenbasis rotation *)
  {msude[[1]], Wude[[1]], M2ude[[1]]} = GetSfermionMixing[
    \[Psi]u,
    ("VCKM" /. x$).(SCKMBasis["m2Q"] /. x$).ConjugateTranspose["VCKM" /. x$], SCKMBasis["m2u"] /. x$, "\[Mu]" /. x$, 
    (SCKMBasis["T"] /. x$)[[1]], DiagonalMatrix[(SCKMBasis["Y"] /. x$)[[1]]], 
    "gp" /. x$, "g2" /. x$, "vev" /. x$,
    ArcTan["tan\[Beta]" /. x$]
  ];
  {msude[[2]], Wude[[2]], M2ude[[2]]} = GetSfermionMixing[
    \[Psi]d,
    SCKMBasis["m2Q"] /. x$, SCKMBasis["m2d"] /. x$, "\[Mu]" /. x$, 
    (SCKMBasis["T"] /. x$)[[2]], DiagonalMatrix[(SCKMBasis["Y"] /. x$)[[2]]], 
    "gp" /. x$, "g2" /. x$, "vev" /. x$,
    ArcTan["tan\[Beta]" /. x$]
  ];
  {msude[[3]], Wude[[3]], M2ude[[3]]} = GetSfermionMixing[
    \[Psi]e,
    SCKMBasis["m2L"] /. x$, SCKMBasis["m2e"] /. x$, "\[Mu]" /. x$, 
    (SCKMBasis["T"] /. x$)[[3]], DiagonalMatrix[(SCKMBasis["Y"] /. x$)[[3]]], 
    "gp" /. x$, "g2" /. x$, "vev" /. x$,
    ArcTan["tan\[Beta]" /. x$]
  ];
  {ms\[Nu], W\[Nu], M2\[Nu]} = GetSfermionMixing[
    \[Psi]\[Nu],
    ConjugateTranspose["VPMNS" /. x$].(SCKMBasis["m2L"] /. x$).("VPMNS" /. x$), (*RH \[Nu] integrated Out*) Zero3, "\[Mu]" /. x$,
    Zero3, Zero3,
    "gp" /. x$, "g2" /. x$, "vev" /. x$,
    ArcTan["tan\[Beta]" /. x$]
  ];
  
  
  CombineRules[
    Join[
      {
        "Wude" -> Wude,
        "msude" -> msude,
	"M2ude" -> M2ude,
        "W\[Nu]" -> W\[Nu],
        "ms\[Nu]" -> ms\[Nu],
	"M2\[Nu]" -> M2\[Nu]
      },
    
      (* neutralino, chargino and Higgs mass eigenbasis rotation *)
      Thread[
        {"m\[Chi]0", "N"} -> GetNeutralinoMixing[
          "M1" /. x$, "M2" /. x$, "\[Mu]" /. x$,
          "gp" /. x$, "g2" /. x$, "vev" /. x$, ArcTan["tan\[Beta]" /. x$]
        ]
      ],
      Thread[
        {"m\[Chi]p", "U", "V"} -> GetCharginoMixing[
          "M2" /. x$, "\[Mu]" /. x$, "g2" /. x$, "vev" /. x$, ArcTan["tan\[Beta]" /. x$]
        ]
      ],
      Thread[
        {"Mh", "MH", "tan\[Alpha]"} -> GetCPEvenHiggsMixing[
          "MA" /. x$, "gp" /. x$, "g2" /. x$, "vev" /. x$, ArcTan["tan\[Beta]" /. x$]
        ]
      ]
    ],
    x$
  ]
];

OneLoop\[Mu][
  MZ_, MW_, m2Hd_, m2Hu_, \[Mu]_, Mh_, MH_, MA_, MHp_, \[Alpha]_, \[Beta]_, signMu_, 
  \[Theta]W_, gp_, g2_, v_, 
  y_, A_, mf_, Wude_, W\[Nu]_, msude_, ms\[Nu]_,
  N_, m\[Chi]0_, U_, V_, m\[Chi]p_, Q_
] := signMu * MustBeReal[
  Sqrt[ 1/2 (
    Tan[2 \[Beta]] (
      ( 
        m2Hu
        - tadtu[ \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
          MZ, MW, \[Mu], Mh, MH, MA, MHp, 
          y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
          N, m\[Chi]0, U, V, m\[Chi]p, Q
        ]
      ) Tan[\[Beta]]
      - (
        m2Hd 
        - tadtd[ \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
          MZ, MW, \[Mu], Mh, MH, MA, MHp, 
          y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
          N, m\[Chi]0, U, V, m\[Chi]p, Q]
      ) Cot[\[Beta]]
    ) 
    - MZ^2 
    - Re[\[CapitalPi]ZZT[
      MZ, \[Theta]W, gp, g2, v, \[Alpha], \[Beta], MZ, MW, Mh, MH, MA, MHp, 
      mf, Wude, W\[Nu], msude, ms\[Nu],
      N, m\[Chi]0, U, V, m\[Chi]p, Q
    ]]
    
  )],
  "|\[Mu]|"
];

OneLoopm3[
  MZ_, MW_, m2Hd_, m2Hu_, \[Mu]_, Mh_, MH_, MA_, MHp_, \[Alpha]_, \[Beta]_, signMu_, 
  \[Theta]W_, gp_, g2_, v_, 
  y_, A_, mf_, Wude_, W\[Nu]_, msude_, ms\[Nu]_,
  N_, m\[Chi]0_, U_, V_, m\[Chi]p_, Q_
] := MustBeReal[
  Sqrt[ 1/2 (
    Tan[2 \[Beta]] (
      m2Hu
        - tadtu[ \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
          MZ, MW, \[Mu], Mh, MH, MA, MHp, 
          y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
          N, m\[Chi]0, U, V, m\[Chi]p, Q
        ]
      - m2Hd 
        - tadtd[ \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
          MZ, MW, \[Mu], Mh, MH, MA, MHp, 
          y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
          N, m\[Chi]0, U, V, m\[Chi]p, Q]
      
    ) 
    - (
	MZ^2 
	+ Re[\[CapitalPi]ZZT[
      	   MZ, \[Theta]W, gp, g2, v, \[Alpha], \[Beta], MZ, MW, Mh, MH, MA, MHp, 
	   mf, Wude, W\[Nu], msude, ms\[Nu], N, m\[Chi]0, U, V, m\[Chi]p, Q
    	   ]
	]
      ) Sin[2 \[Beta]]
    )
  ],
  "|m3|"
];

OneLoopMA[
  MZ_, MW_, m2Hd_, m2Hu_, \[Mu]_, Mh_, MH_, MA_, MHp_, \[Alpha]_, \[Beta]_, signMu_, \[Theta]W_, 
  gp_, g2_, v_,
  y_, A_, mf_, Wude_, W\[Nu]_, msude_, ms\[Nu]_,
  N_, m\[Chi]0_, U_, V_, m\[Chi]p_, Q_
] := MustBeReal[
  Sqrt[
    1/Cos[2 \[Beta]] (
      m2Hu 
      - Re[tadtu[
        \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
        MZ, MW, \[Mu], Mh, MH, MA, MHp, 
        y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
        N, m\[Chi]0, U, V, m\[Chi]p, Q
      ]]
      - m2Hd 
      + Re[tadtd[
        \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
        MZ, MW, \[Mu], Mh, MH, MA, MHp, 
        y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
        N, m\[Chi]0, U, V, m\[Chi]p, Q
      ]]
    ) 
    - MZ^2
    - Re[\[CapitalPi]ZZT[
      MZ, \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
      MZ, MW, Mh, MH, MA, MHp,
      mf, Wude, W\[Nu], msude, ms\[Nu],
      N, m\[Chi]0, U, V, m\[Chi]p, Q
    ]]
    - Re[\[CapitalPi]AA[
      MA, \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
      MZ, MW, Mh, MH, \[Mu], MA, MHp, 
      y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
      N, m\[Chi]0, U, V, m\[Chi]p, Q
    ]]
    + Sin[\[Beta]]^2 Re[tadtd[ 
      \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
      MZ, MW, \[Mu], Mh, MH, MA, MHp, 
      y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
      N, m\[Chi]0, U, V, m\[Chi]p, Q
    ]]
    + Cos[\[Beta]]^2 Re[tadtu[
      \[Theta]W, gp, g2, v, \[Alpha], \[Beta],
      MZ, MW, \[Mu], Mh, MH, MA, MHp, 
      y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
      N, m\[Chi]0, U, V, m\[Chi]p, Q
    ]]
  ],
  "MA"
];

OneLoopMHp[
  MZ_, MW_, m2Hd_, m2Hu_, \[Mu]_, Mh_, MH_, MA_, MHp_, \[Alpha]_, \[Beta]_, signMu_, \[Theta]W_,
  gp_, g2_, v_,
  y_, A_, mf_, Wude_, W\[Nu]_, msude_, ms\[Nu]_,
  N_, m\[Chi]0_, U_, V_, m\[Chi]p_, Q_
] := With[{$MW = g2 v/2},
  MustBeReal[
    Sqrt[
      1/Cos[2 \[Beta]] (
      m2Hu 
      - Re[tadtu[
        \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
        MZ, MW, \[Mu], Mh, MH, MA, MHp, 
        y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
        N, m\[Chi]0, U, V, m\[Chi]p, Q
      ]]
      - m2Hd 
      + Re[tadtd[
        \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
        MZ, MW, \[Mu], Mh, MH, MA, MHp, 
        y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
        N, m\[Chi]0, U, V, m\[Chi]p, Q
      ]]
    ) 
    - MZ^2
    - Re[\[CapitalPi]ZZT[
      MZ, \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
      MZ, MW, Mh, MH, MA, MHp,
      mf, Wude, W\[Nu], msude, ms\[Nu],
      N, m\[Chi]0, U, V, m\[Chi]p, Q
    ]]
    + $MW^2
    -Re[\[CapitalPi]HpHm[
     MHp, \[Theta]W, gp, g2, v, \[Alpha], \[Beta],
     MZ, MW, Mh, MH, \[Mu], MA, MHp,
     y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
     N, m\[Chi]0, U, V, m\[Chi]p, Q
    ]]
    + Sin[\[Beta]]^2 Re[tadtd[ 
      \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
      MZ, MW, \[Mu], Mh, MH, MA, MHp, 
      y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
      N, m\[Chi]0, U, V, m\[Chi]p, Q
    ]]
    + Cos[\[Beta]]^2 Re[tadtu[
      \[Theta]W, gp, g2, v, \[Alpha], \[Beta],
      MZ, MW, \[Mu], Mh, MH, MA, MHp, 
      y, A, mf, Wude, W\[Nu], msude, ms\[Nu],
      N, m\[Chi]0, U, V, m\[Chi]p, Q
    ]]
  ],
  "MHp"
]];
  
OneLoopVEV[MZ_, MW_, m2Hd_, m2Hu_, \[Mu]_, Mh_, MH_, MA_, MHp_, \[Alpha]_, \[Beta]_, \[Theta]W_, 
  gp_, g2_, v_, 
  y_, A_, mf_, Wude_, W\[Nu]_, msude_, ms\[Nu]_,
  N_, m\[Chi]0_, U_, V_, m\[Chi]p_, Q_
] := MustBeReal[
  2 Sqrt[
    (
      MZ^2
      + Re[\[CapitalPi]ZZT[
        MZ, \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
        MZ, MW, Mh, MH, MA, MHp, 
        mf, Wude, W\[Nu], msude, ms\[Nu], 
        N, m\[Chi]0, U, V, m\[Chi]p, Q
      ]]
    ) / (gp^2 + g2^2)
  ],
  "vev"
];

OneLoopMZ[MZ_, MW_, Mh_, MH_, MA_, MHp_, \[Alpha]_, \[Beta]_, \[Theta]W_, 
  gp_, g2_, v_,
  mf_, Wude_, W\[Nu]_, msude_, ms\[Nu]_,
  N_, m\[Chi]0_, U_, V_, m\[Chi]p_, Q_
] := MustBeReal[
  Sqrt[
    MZ^2
    + Re[\[CapitalPi]ZZT[
      MZ, \[Theta]W, gp, g2, v, \[Alpha], \[Beta], 
      MZ, MW, Mh, MH, MA, MHp, 
      mf, Wude, W\[Nu], msude, ms\[Nu],
      N, m\[Chi]0, U, V, m\[Chi]p, Q
    ]]
  ],
  "MZ"
];

CalcOneLoopValues[x_] := Block[{iterate},

  (*Print["tree \[Mu]" -> ("\[Mu]" /. x)];*)
  iterate = Function[{xin}, Block[{xout},

    (*Print["iteration"];*)
    xout = CalcTreeValues[xin];
	
    (*Print["Calculate One Loop Values..."];*)
    xout = Join[
      FilterRules[xout, If[("cMSSM" /. xout) == True, Except["vev" | "MHp" | "\[Mu]" | "m3"], Except["vev" | "MA" | "\[Mu]" | "m3"]]],
      {
        "vev" -> OneLoopVEV[
          Re["MZ" /. xout], Re["MW" /. xout],
          Re["m2Hd" /. xout], Re["m2Hu" /. xout], "\[Mu]" /. xout,
          Re["Mh" /. xout], Re["MH" /. xout], Re["MA" /. xout], Re["MHp" /. xout], 
          Re[ArcTan["tan\[Alpha]" /. xout]], Re[ArcTan["tan\[Beta]" /. xout]], 
          Re["\[Theta]W" /. xout], Re["gp" /. xout], Re["g2" /. xout], Re["vev" /. xout], 
          SCKMBasis["Y"] /. xout, SCKMBasis["T"] /. xout, Re["mf" /. xout], 
          "Wude" /. xout, "W\[Nu]" /. xout, Re["msude" /. xout], Re["ms\[Nu]" /. xout],
          "N" /. xout, Re["m\[Chi]0" /. xout], "U" /. xout, "V" /. xout, Re["m\[Chi]p" /. xout], Re["Q" /. xout]
        ],
        If[      
          ("cMSSM" /. xout) == True,
          "MHp" -> OneLoopMHp[
            Re["MZ" /. xout], Re["MW" /. xout], 
            Re["m2Hd" /. xout], Re["m2Hu" /. xout], "\[Mu]" /. xout,
            Re["Mh" /. xout], Re["MH" /. xout], Re["MA" /. xout], Re["MHp" /. xout],
            Re[ArcTan["tan\[Alpha]" /. xout]], Re[ArcTan["tan\[Beta]" /. xout]], "sign\[Mu]" /. xout,
            Re["\[Theta]W" /. xout], Re["gp" /. xout], Re["g2" /. xout], Re["vev" /. xout],
            SCKMBasis["Y"] /. xout, SCKMBasis["T"] /. xout, 
            Re["mf" /. xout], "Wude" /. xout, "W\[Nu]" /. xout, Re["msude" /. xout], Re["ms\[Nu]" /. xout],
            "N" /. xout, Re["m\[Chi]0" /. xout], "U" /. xout, "V" /. xout, Re["m\[Chi]p" /. xout], Re["Q" /. xout]
          ],
          (*else cMSSM == False*)
	  "MA" -> OneLoopMA[
            Re["MZ" /. xout], Re["MW" /. xout], 
            Re["m2Hd" /. xout], Re["m2Hu" /. xout], "\[Mu]" /. xout,
            Re["Mh" /. xout], Re["MH" /. xout], Re["MA" /. xout], Re["MHp" /. xout],
            Re[ArcTan["tan\[Alpha]" /. xout]], Re[ArcTan["tan\[Beta]" /. xout]], "sign\[Mu]" /. xout,
            Re["\[Theta]W" /. xout], Re["gp" /. xout], Re["g2" /. xout], Re["vev" /. xout],
            SCKMBasis["Y"] /. xout, SCKMBasis["T"] /. xout, 
            Re["mf" /. xout], "Wude" /. xout, "W\[Nu]" /. xout, Re["msude" /. xout], Re["ms\[Nu]" /. xout],
            "N" /. xout, Re["m\[Chi]0" /. xout], "U" /. xout, "V" /. xout, Re["m\[Chi]p" /. xout], Re["Q" /. xout]
          ]
        ],
	"\[Mu]" -> OneLoop\[Mu][
          Re["MZ" /. xout], Re["MW" /. xout], 
          Re["m2Hd" /. xout], Re["m2Hu" /. xout], "\[Mu]" /. xout, 
          Re["Mh" /. xout], Re["MH" /. xout], Re["MA" /. xout], Re["MHp" /. xout], 
          Re[ArcTan["tan\[Alpha]" /. xout]], Re[ArcTan["tan\[Beta]" /. xout]], "sign\[Mu]" /. xout, 
          Re["\[Theta]W" /. xout], Re["gp" /. xout], Re["g2" /. xout], Re["vev" /. xout], 
          SCKMBasis["Y"] /. xout, SCKMBasis["T"] /. xout,
          Re["mf" /. xout], "Wude" /. xout, "W\[Nu]" /. xout, Re["msude" /. xout], Re["ms\[Nu]" /. xout],
          "N" /. xout, Re["m\[Chi]0" /. xout], "U" /. xout, "V" /. xout, Re["m\[Chi]p" /. xout], Re["Q" /. xout]
        ],
	"m3" -> OneLoopm3[
          Re["MZ" /. xout], Re["MW" /. xout], 
          Re["m2Hd" /. xout], Re["m2Hu" /. xout], "\[Mu]" /. xout, 
          Re["Mh" /. xout], Re["MH" /. xout], Re["MA" /. xout], Re["MHp" /. xout], 
          Re[ArcTan["tan\[Alpha]" /. xout]], Re[ArcTan["tan\[Beta]" /. xout]], "sign\[Mu]" /. xout, 
          Re["\[Theta]W" /. xout], Re["gp" /. xout], Re["g2" /. xout], Re["vev" /. xout], 
          SCKMBasis["Y"] /. xout, SCKMBasis["T"] /. xout,
          Re["mf" /. xout], "Wude" /. xout, "W\[Nu]" /. xout, Re["msude" /. xout], Re["ms\[Nu]" /. xout],
          "N" /. xout, Re["m\[Chi]0" /. xout], "U" /. xout, "V" /. xout, Re["m\[Chi]p" /. xout], Re["Q" /. xout]
        ]
      }
    ]
  ]];
  
  FixedPoint[
    iterate, x, 1 (*If[("cMSSM" /. x) == True, 1, 100]*),
    SameTest -> (
      Abs[("vev" /. #2) / ("vev" /. #1) - 1] < 1*^-4 &&
      Abs[("MA" /. #2) / ("MA" /. #1) - 1] < 1*^-4 &&
      Abs[("MHp" /. #2) / ("MHp" /. #1) - 1] < 1*^-4 &&
      Abs[("Mh" /. #2) / ("Mh" /. #1) - 1] < 1*^-4 &&
      Abs[("MH" /. #2) / ("MH" /. #1) - 1] < 1*^-4 &&
      Abs[("\[Mu]" /. #2) / ("\[Mu]" /. #1) - 1] < 1*^-4 &&
      Abs[("m3" /. #2) / ("m3" /. #1) - 1] < 1*^-4 &
    )
  ]
];

(*Tachyon Part*)

TachyonFree[m2Q_InterpolatingFunction, m2u_InterpolatingFunction,  m2d_InterpolatingFunction, m2L_InterpolatingFunction, m2e_InterpolatingFunction, tmin_, tmax_]:=Block[{sol,t,tab},
tab=Table[{t,Det[m2Q[t]//Chop]//Re},{t,tmin,tmax,(tmax-tmin)/1000}];
sol=Select[tab,#[[2]]<0&];
If[sol!={},Throw[{"m2Q tachyonic at {Log[scale],Det[m2Q]}", sol[[-1]]}, TachyonicMass]];
tab=Table[{t,Det[m2u[t]//Chop]//Re},{t,tmin,tmax,(tmax-tmin)/1000}];
sol=Select[tab,#[[2]]<0&];
If[sol!={},Throw[{"m2u tachyonic at {Log[scale],Det[m2u]}", sol[[-1]]}, TachyonicMass]];
tab=Table[{t,Det[m2d[t]//Chop]//Re},{t,tmin,tmax,(tmax-tmin)/1000}];
sol=Select[tab,#[[2]]<0&];
If[sol!={},Throw[{"m2d tachyonic at {Log[scale],Det[m2d]}", sol[[-1]]}, TachyonicMass]];
tab=Table[{t,Det[m2L[t]//Chop]//Re},{t,tmin,tmax,(tmax-tmin)/1000}];
sol=Select[tab,#[[2]]<0&];
If[sol!={},Throw[{"m2L tachyonic at {Log[scale],Det[m2L]}", sol[[-1]]}, TachyonicMass]];
tab=Table[{t,Det[m2e[t]//Chop]//Re},{t,tmin,tmax,(tmax-tmin)/1000}];
sol=Select[tab,#[[2]]<0&];
If[sol!={},Throw[{"m2e tachyonic at {Log[scale],Det[m2e]}", sol[[-1]]}, TachyonicMass]];
];

(*SUSYScale Part*)

SUSYScale[Wu_?MatrixQ, msu_?VectorQ, fromStops_?BooleanQ] := If[fromStops,
	GeometricMean[
 	  (* stop1 and stop2 correspond to the up squark i with Wu[[{3,6}, i]] maximal *)
	  msu[[SortBy[Range[6], -Norm[Wu[[{3, 6}, #]]] &][[{1, 2}]]]]
	],
	GeometricMean[
 	  (* take lightest and heaviest up-type squarks *)
	  {msu[[1]],msu[[6]]}
	]
];

SUSYScale[x_,fromStops_?BooleanQ] := With[
  {
    Wu=("Wude"/.x)[[1]], msu = ("msude"/.x)[[1]]
  },
  SUSYScale[Wu, msu, fromStops]
]

SUSYScale[m2Q_?MatrixQ, m2u_?MatrixQ, Tu_?MatrixQ, Yu_?MatrixQ, Yd_?MatrixQ, g1_, g2_, \[Beta]_, m2Hd_, m2Hu_, signMu_, fromStops_?BooleanQ] := Block[{Wu, msu, M2u, UuL, yu, 
UuR, UdL, yd, UdR,m2Qscm,m2usckm,Tusckm},
  (* Beware of convention conversions, cf. SCKMRotate! *)
  {UuL, yu, UuR, UdL, yd, UdR} = STCGetQuarkMixing[Transpose[Yu],Transpose[Yd]];
  m2Qsckm = ConjugateTranspose[UdL].m2Q.UdL;
  m2usckm = ConjugateTranspose[UuR].m2u.UuR;
  Tusckm = ConjugateTranspose[UuR].Transpose[Tu].UuL;
  {msu, Wu, M2u} = GetSfermionMixing[
    \[Psi]u,
    m2Qsckm, m2usckm,
    TreeLevel\[Mu][MZ, m2Hd, m2Hu, \[Beta], signMu],
    Tusckm, DiagonalMatrix[yu], 
    g1 Sqrt[3/5], g2, TreeLevelVEV[MZ, g1 Sqrt[3/5], g2],
    \[Beta]
  ];
  SUSYScale[Wu, msu, fromStops]
];


(*Because of numerical stability, improve FindRoot:*)
SUSYScale[
  m2Q_InterpolatingFunction, m2u_InterpolatingFunction, Au_InterpolatingFunction, Yu_InterpolatingFunction, Yd_InterpolatingFunction,
  g1_InterpolatingFunction, g2_InterpolatingFunction, \[Beta]_, 
  m2Hd_InterpolatingFunction, m2Hu_InterpolatingFunction, signMu_, fromStops_
] := Block[
  {t},
  E^t /. Quiet[
    FindRoot[
      Log[SUSYScale[m2Q[t], m2u[t], Au[t], Yu[t], Yd[t], g1[t], g2[t], \[Beta], m2Hd[t], m2Hu[t], signMu, fromStops]] - t, 
      {t, Log[1000]}
    ],
    InterpolatingFunction::dmval
  ]
]
  
(* SUSY threshold correction part *)

H2[x$_, y$_] := With[{x=Abs[x$],y=Abs[y$]},x^2 Log[x^2]/((1 - x^2) (x^2 - y^2)) + y^2 Log[y^2]/((1 - y^2) (y^2 - x^2))];
C00[q$_,x$_,y$_] := With[{q=Abs[q$],x=Abs[x$],y=Abs[y$]},  1/4 (3/2 - Log[q^2] + x^4 Log[x^2]/((1 - x^2) (x^2 - y^2))+ y^4 Log[y^2]/((1 - y^2) (y^2 - x^2)))];

B1[q$_] := With[{q=Abs[q$]}, -3/4 + 1/2 Log[q^2]];
B1t[q$_] := B1t[q$] = With[{q=Abs[q$]}, -1/4 + 1/2 Log[q^2]];
B2[m1$_,m2$_] := With[{m1=Abs[m1$],m2=Abs[m2$]},1/2 (m2^4-m1^4)/(m1^2-m2^2)^3 + m1^2 m2^2 /(m1^2-m2^2)^3 Log[m1^2/m2^2]];
B2t[Q$_,m1$_,m2$_] := With[{Q=Abs[Q$],m1=Abs[m1$],m2=Abs[m2$]}, 5 (m1^6 - m2^6) + 27 m1^2 m2^2 (m2^2 - m1^2) + 6 m2^4 (m2^2 - 3 m1^2) Log[m2^2/m1^2] + 6 (m2^2 - m1^2)^3 Log[m1^2/Q^2]];
B3[Q$_,m1$_,m2$_] := With[{Q=Abs[Q$],m1=Abs[m1$],m2=Abs[m2$]},-1/4 - 1/2 m1^2/(m1^2-m2^2) + 1/2 m1^4 /(m1^2-m2^2)^2 Log[m1^2/m2^2] + 1/2 Log[m2^2/Q^2]];



GetSoftMassMixing[Msq$_] := GetSoftMassMixing[Msq$] = Block[{msf, Wf, Msq},
  Msq=MustBeHermitian[Msq$];
  {msf, Wf} = Eigensystem[Msq];
  msf = MustBeReal[Reverse[Sqrt[msf]], "scalar soft mass"]; (*Note: Yields the scalar masses ms and NOT ms^2*)
  Wf = Transpose[Reverse[Wf]];
  {msf, Wf}
];

GetHiggsDoubletMass[x_] := GetHiggsDoubletMass[x] = Block[{mt$, Wt$, M$, gspec=REAP`RGEMSSMsoftbroken`Private`gSUSYSpectrum},
	M$ = {{Abs["\[Mu]" /.gspec]^2 + "m2Hd" /.x //Chop, -("m3" /.gspec)^2},
	      {-Conjugate["m3" /.gspec]^2, Abs["\[Mu]" /.gspec]^2 + "m2Hu" /.x//Chop}};
	{mt$, Wt$} = Eigensystem[M$];
	{mt$, Wt$} = {Reverse[mt$], Reverse[Wt$]};
 	Wt$ = RotationMatrix[ArcTan[Wt$[[1, 2]]/Wt$[[2, 2]]]];
	mt$ = mt$[[2]]//Sqrt;
	{mt$, Wt$}
];

GetCN[K2$_]:=Block[{K2diag, U, K2},
   K2=MustBeHermitian[K2$];
   {K2diag,U}=Eigensystem[K2];
   U = Conjugate[U];
   (*{Reverse[K2diag], Reverse[U]}*)
   {K2diag,U}
]

hwave[x_] := hwave[x] = Block[{
(* Convert to LR *)
    Yd = Transpose[FlavourBasis["Yd"] /. x],
    Td = Transpose[FlavourBasis["Td"] /. x],
    Yu = Transpose[FlavourBasis["Yu"] /. x],
    Tu = Transpose[FlavourBasis["Tu"] /. x],
    Ye = Transpose[FlavourBasis["Ye"] /. x],
    Te = Transpose[FlavourBasis["Te"] /. x], 
    rotYd, rotYu, rotTu, rotTd, rotYe, rotTe, msq, msu, msd, Wq, Wu, Wd,
    msl, mse, Wl, We, mH, U, X11, X12, X22, X22t, Y11, Y11t, Y22, Z12, Z12t, XA11, XB11, XA12, XB12, XA22, XB22,
    gp, g2, g3, M1, M2, M3, \[Mu], Q, \[Delta]
  },
  {msl, Wl} = GetSoftMassMixing[FlavourBasis["m2L"] /. x];
  {mse, We} = GetSoftMassMixing[FlavourBasis["m2e"] /. x];
  {msq, Wq} = GetSoftMassMixing[FlavourBasis["m2Q"] /. x];
  {msu, Wu} = GetSoftMassMixing[FlavourBasis["m2u"] /. x];
  {msd, Wd} = GetSoftMassMixing[FlavourBasis["m2d"] /. x];
  rotYe = Transpose[Wl].Ye.Conjugate[We];
  rotTe = Transpose[Wl].Te.Conjugate[We];
  rotYu = Transpose[Wq].Yu.Conjugate[Wu];
  rotTu = Transpose[Wq].Tu.Conjugate[Wu];
  rotYd = Transpose[Wq].Yd.Conjugate[Wd];
  rotTd = Transpose[Wq].Td.Conjugate[Wd];
  {mH, U} = GetHiggsDoubletMass[x];
  {gp, g2, g3, M1, M2, M3, \[Mu], Q, NonTan\[Beta]} = {"gp", "g2", "g3", "M1", "M2", "M3", "\[Mu]", "Q", "NonTan\[Beta]"} /. x;
	X11 = Loop Abs[\[Mu]]^2 Sum[Abs[rotYu[[i,j]]]^2 B2[msq[[i]], msu[[j]]],{i,3},{j,3}]//Chop;
	X12 = - Loop \[Mu] Sum[rotTu[[i,j]] Conjugate[rotYu[[i,j]]] B2[msq[[i]], msu[[j]]],{i,3},{j,3}]//Chop;
	Y11 = Loop Sum[Abs[rotTd[[i,j]]]^2 B2[msq[[i]], msd[[j]]],{i,3},{j,3}]//Chop;
	Z12 = - Loop \[Mu] Sum[rotTd[[i,j]] Conjugate[rotYd[[i,j]]] B2[msq[[i]], msd[[j]]],{i,3},{j,3}]//Chop;
	X22 = Loop Abs[\[Mu]]^2 Sum[Abs[rotYd[[i,j]]]^2 B2[msq[[i]], msd[[j]]],{i,3},{j,3}]//Chop;
	Y22 = Loop Sum[Abs[rotTu[[i,j]]]^2 B2[msq[[i]], msu[[j]]],{i,3},{j,3}]//Chop;
	Y11t = Loop Sum[Abs[rotTe[[i,j]]]^2 B2[msl[[i]], mse[[j]]],{i,3},{j,3}]//Chop;
	Z12t = - Loop \[Mu] Sum[rotTe[[i,j]] Conjugate[rotYe[[i,j]]] B2[msl[[i]], mse[[j]]],{i,3},{j,3}]//Chop;
	X22t = Loop Abs[\[Mu]]^2 Sum[Abs[rotYe[[i,j]]]^2 B2[msl[[i]], mse[[j]]],{i,3},{j,3}]//Chop;
	XA11 = XA22 = Loop gp^2 1/12 1/(Abs[\[Mu]]^2-Abs[M1]^2)^3 B2t[Q, M1, \[Mu]]//Chop;
	XB11 = XB22 = Loop g2^2 1/4 1/(Abs[\[Mu]]^2-Abs[M2]^2)^3 B2t[Q, M2, \[Mu]]//Chop;
	XA12 = Loop gp^2 \[Mu] M1 B2[\[Mu],M1]//Chop;
	XB12 = Loop g2^2 3 \[Mu] M2 B2[\[Mu],M2]//Chop;
  \[Delta] = 	{{1 + X11 + Y11 + Y11t + XA11 + XB11, X12 + Z12 + Z12t + XA12 + XB12},
		{Conjugate[X12 + Z12 + Z12t +  XA12 + XB12], 1 + X22 + X22t + Y22 + XA22 + XB22}};
  If[NonTan\[Beta],Sqrt[1/(U.\[Delta].ConjugateTranspose[U])[[1,1]]//Chop//N],1]
];

Qwave[x_] := Qwave[x] =Block[{
(* Convert to LR *)
    Yd = Transpose[FlavourBasis["Yd"] /. x],
    Td = Transpose[FlavourBasis["Td"] /. x],
    Yu = Transpose[FlavourBasis["Yu"] /. x],
    Tu = Transpose[FlavourBasis["Tu"] /. x],
    rotYd, rotYu, rotTu, rotTd, msq, msu, msd, Wq, Wu, Wd,
    mH, U, r1, r2, r3, rd, rdt, ru, rut,
    gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta], \[Delta], Kq, Uq
  },
  {msq, Wq} = GetSoftMassMixing[FlavourBasis["m2Q"] /. x];
  {msu, Wu} = GetSoftMassMixing[FlavourBasis["m2u"] /. x];
  {msd, Wd} = GetSoftMassMixing[FlavourBasis["m2d"] /. x];
  rotYu = Transpose[Wq].Yu.Conjugate[Wu];
  rotTu = Transpose[Wq].Tu.Conjugate[Wu];
  rotYd = Transpose[Wq].Yd.Conjugate[Wd];
  rotTd = Transpose[Wq].Td.Conjugate[Wd];
  {mH, U} = GetHiggsDoubletMass[x];
  {gp, g2, g3, M1, M2, M3, \[Mu], Q,tan\[Beta], NonTan\[Beta]} = {"gp", "g2", "g3", "M1", "M2", "M3", "\[Mu]", "Q", "tan\[Beta]" , "NonTan\[Beta]"} /. x;
	r3 = DiagonalMatrix[Table[8 /3 g3^2 Loop B3[Q,M3,msq[[i]]],{i,3}]]//Chop;
	r2 = DiagonalMatrix[Table[3 /2 g2^2 Loop B3[Q,M2,msq[[i]]],{i,3}]]//Chop;
	r1 = DiagonalMatrix[Table[1 /18 gp^2 Loop B3[Q,M1,msq[[i]]],{i,3}]]//Chop;
	rd = Loop Table[Sum[rotYd[[i,k]] ConjugateTranspose[rotYd][[k,j]] B3[Q,\[Mu],msd[[k]]] ,{k,3}],{i,3},{j,3}]//Chop;
	rdt = Loop Sin[ArcTan[tan\[Beta]]]^2 Table[Sum[rotYd[[i,k]] ConjugateTranspose[rotYd][[k,j]] B1t[mH/Q] ,{k,3}],{i,3},{j,3}]//Chop;
	ru =  Loop Table[Sum[rotYu[[i,k]] ConjugateTranspose[rotYu][[k,j]] B3[Q,\[Mu],msu[[k]]] ,{k,3}],{i,3},{j,3}]//Chop;
	rut = Loop Cos[ArcTan[tan\[Beta]]]^2 Table[Sum[rotYu[[i,k]] ConjugateTranspose[rotYu][[k,j]] B1t[mH/Q] ,{k,3}],{i,3},{j,3}]//Chop;
  \[Delta] = IdentityMatrix[3] + r3 + r2 + r1 + rd + rdt + ru + rut;
  {Kq,Uq} = GetCN[\[Delta]];
  If[NonTan\[Beta],Sqrt[Inverse[DiagonalMatrix[Kq]]].Uq//Chop//N,IdentityMatrix[3]]
];

dwave[x_] := dwave[x] =Block[{
(* Convert to LR *)
    Yd = Transpose[FlavourBasis["Yd"] /. x],
    Td = Transpose[FlavourBasis["Td"] /. x],
    rotYd, rotTd, msq, msu, msd, Wq, Wu, Wd,
    mH, U, r1,  r3, rd, rdt,
    gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta], \[Delta], Kd, Ud
  },
  {msq, Wq} = GetSoftMassMixing[FlavourBasis["m2Q"] /. x];
  {msu, Wu} = GetSoftMassMixing[FlavourBasis["m2u"] /. x];
  {msd, Wd} = GetSoftMassMixing[FlavourBasis["m2d"] /. x];
  rotYd = Transpose[Wq].Yd.Conjugate[Wd];
  rotTd = Transpose[Wq].Td.Conjugate[Wd];
  {mH, U} = GetHiggsDoubletMass[x];
  {gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta], NonTan\[Beta]} = {"gp", "g2", "g3", "M1", "M2", "M3", "\[Mu]", "Q", "tan\[Beta]", "NonTan\[Beta]"} /. x;
	r3 = DiagonalMatrix[Table[8 /3 g3^2 Loop B3[Q,M3,msd[[j]]],{j,3}]]//Chop;
	r1 = DiagonalMatrix[Table[2 /9 gp^2 Loop B3[Q,M1,msd[[j]]],{j,3}]]//Chop;
	rd = 2 Loop Table[Sum[ConjugateTranspose[rotYd][[j,k]] rotYd[[k,i]] B3[Q,\[Mu],msq[[k]]] ,{k,3}],{i,3},{j,3}]//Chop;
	rdt = 2 Loop Sin[ArcTan[tan\[Beta]]]^2 Table[Sum[ConjugateTranspose[rotYd][[j,k]] rotYd[[k,i]] B1t[mH/Q] ,{k,3}],{i,3},{j,3}]//Chop;
  \[Delta] = IdentityMatrix[3]+ r3 + r1 + rd + rdt;
  {Kd,Ud} = GetCN[\[Delta]];
  If[NonTan\[Beta],Sqrt[Inverse[DiagonalMatrix[Kd]]].Ud//Chop//N,IdentityMatrix[3]]
];

uwave[x_] := uwave[x] =Block[{
(* Convert to LR *)
    Yu = Transpose[FlavourBasis["Yu"] /. x],
    Tu = Transpose[FlavourBasis["Tu"] /. x],
    rotYu, rotTu, msq, msu, msd, Wq, Wu, Wd,
    mH, U, r1,  r3, ru, rut,
    gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta], \[Delta], Ku, Uu
  },
  {msq, Wq} = GetSoftMassMixing[FlavourBasis["m2Q"] /. x];
  {msu, Wu} = GetSoftMassMixing[FlavourBasis["m2u"] /. x];
  {msd, Wd} = GetSoftMassMixing[FlavourBasis["m2d"] /. x];
  rotYu = Transpose[Wq].Yu.Conjugate[Wu];
  rotTu = Transpose[Wq].Tu.Conjugate[Wu];
  {mH, U} = GetHiggsDoubletMass[x];
  {gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta],NonTan\[Beta]} = {"gp", "g2", "g3", "M1", "M2", "M3", "\[Mu]", "Q", "tan\[Beta]","NonTan\[Beta]"} /. x;
	r3 = DiagonalMatrix[Table[8 /3 g3^2 Loop B3[Q,M3,msu[[j]]],{j,3}]]//Chop;
	r1 = DiagonalMatrix[Table[8 /9 gp^2 Loop B3[Q,M1,msu[[j]]],{j,3}]]//Chop;
	ru = 2 Loop Table[Sum[ConjugateTranspose[rotYu][[j,k]] rotYu[[k,i]] B3[Q,\[Mu],msq[[k]]] ,{k,3}],{i,3},{j,3}]//Chop;
	rut = 2 Loop Cos[ArcTan[tan\[Beta]]]^2 Table[Sum[ConjugateTranspose[rotYu][[j,k]] rotYu[[k,i]] B1t[mH/Q] ,{k,3}],{i,3},{j,3}]//Chop;
  \[Delta] = IdentityMatrix[3]+ r3 + r1 + ru + rut;
  {Ku,Uu} = GetCN[\[Delta]];
  If[NonTan\[Beta],Sqrt[Inverse[DiagonalMatrix[Ku]]].Uu//Chop//N,IdentityMatrix[3]]
];

Lwave[x_] := Lwave[x] =Block[{
(* Convert to LR *)
    Ye = Transpose[FlavourBasis["Ye"] /. x],
    Te = Transpose[FlavourBasis["Te"] /. x],
    rotYe, rotTe, msl, mse, Wl, We,
    mH, U, r1, r2, re, ret,
    gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta], \[Delta], Kl, Ul
  },
  {msl, Wl} = GetSoftMassMixing[FlavourBasis["m2L"] /. x];
  {mse, We} = GetSoftMassMixing[FlavourBasis["m2e"] /. x];
  rotYe = Transpose[Wl].Ye.Conjugate[We];
  rotTe = Transpose[Wl].Te.Conjugate[We];
  {mH, U} = GetHiggsDoubletMass[x];
  {gp, g2, g3, M1, M2, M3, \[Mu], Q,tan\[Beta],NonTan\[Beta]} = {"gp", "g2", "g3", "M1", "M2", "M3", "\[Mu]", "Q", "tan\[Beta]","NonTan\[Beta]"} /. x;
	r2 = DiagonalMatrix[Table[3 /2 g2^2 Loop B3[Q,M2,msl[[i]]],{i,3}]]//Chop;
	r1 = DiagonalMatrix[Table[1 /2 gp^2 Loop B3[Q,M1,msl[[i]]],{i,3}]]//Chop;
	re = Loop Table[Sum[rotYe[[i,k]] ConjugateTranspose[rotYe][[k,j]] B3[Q,\[Mu],mse[[k]]] ,{k,3}],{i,3},{j,3}]//Chop;
	ret = Loop Sin[ArcTan[tan\[Beta]]]^2 Table[Sum[rotYe[[i,k]] ConjugateTranspose[rotYe][[k,j]] B1t[mH/Q] ,{k,3}],{i,3},{j,3}]//Chop;
  \[Delta] = IdentityMatrix[3]+ r2 + r1 + re + ret;
  {Kl,Ul} = GetCN[\[Delta]];
  If[NonTan\[Beta],Sqrt[Inverse[DiagonalMatrix[Kl]]].Ul//Chop//N,IdentityMatrix[3]]
];

ewave[x_] := ewave[x] =Block[{
(* Convert to LR *)
    Ye = Transpose[FlavourBasis["Ye"] /. x],
    Te = Transpose[FlavourBasis["Te"] /. x],
    rotYe, rotTe, msl, mse, Wl, We,
    mH, U, r1, r2, re, ret,
    gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta], \[Delta], Ke, Ue
  },
  {msl, Wl} = GetSoftMassMixing[FlavourBasis["m2L"] /. x];
  {mse, We} = GetSoftMassMixing[FlavourBasis["m2e"] /. x];
  rotYe = Transpose[Wl].Ye.Conjugate[We];
  rotTe = Transpose[Wl].Te.Conjugate[We];
  {mH, U} = GetHiggsDoubletMass[x];
  {gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta],NonTan\[Beta]} = {"gp", "g2", "g3", "M1", "M2", "M3", "\[Mu]", "Q", "tan\[Beta]","NonTan\[Beta]"} /. x;
	r1 = DiagonalMatrix[Table[2 gp^2 Loop B3[Q,M1,mse[[j]]],{j,3}]]//Chop;
	re = 2 Loop Table[Sum[ConjugateTranspose[rotYe][[j,k]] rotYe[[k,i]] B3[Q,\[Mu],msl[[k]]] ,{k,3}],{i,3},{j,3}]//Chop;
	ret = 2 Loop Sin[ArcTan[tan\[Beta]]]^2 Table[Sum[ConjugateTranspose[rotYe][[j,k]] rotYe[[k,i]] B1t[mH/Q] ,{k,3}],{i,3},{j,3}]//Chop;
  \[Delta] = IdentityMatrix[3] + r1 + re + ret;
  {Ke,Ue} = GetCN[\[Delta]];
  If[NonTan\[Beta],Sqrt[Inverse[DiagonalMatrix[Ke]]].Ue//Chop//N,IdentityMatrix[3]]
];

OneLoopSMYd[x_] := Block[{
    (* Convert to LR *)
    Yd = Transpose[FlavourBasis["Yd"] /. x],
    Td = Transpose[FlavourBasis["Td"] /. x],
    Yu = Transpose[FlavourBasis["Yu"] /. x],
    Tu = Transpose[FlavourBasis["Tu"] /. x],
    SMYd, rotYd, rotYu, rotTu, rotTd, msq, msu, msd, Wq, Wu, Wd,
    gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta], mH, U
  },
  {msq, Wq} = GetSoftMassMixing[FlavourBasis["m2Q"] /. x];
  {msu, Wu} = GetSoftMassMixing[FlavourBasis["m2u"] /. x];
  {msd, Wd} = GetSoftMassMixing[FlavourBasis["m2d"] /. x];
  rotYu = Transpose[Wq].Yu.Conjugate[Wu];
  rotTu = Transpose[Wq].Tu.Conjugate[Wu];
  rotYd = Transpose[Wq].Yd.Conjugate[Wd];
  rotTd = Transpose[Wq].Td.Conjugate[Wd];
  {mH, U} = GetHiggsDoubletMass[x];
  {gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta], NonTan\[Beta]} = {"gp", "g2", "g3", "M1", "M2", "M3", "\[Mu]", "Q", "tan\[Beta]","NonTan\[Beta]"} /. x;
  NonTan\[Beta] = Boole[NonTan\[Beta]];
  SMYd=Table[
	 rotYd[[i, j]] Cos[ArcTan[tan\[Beta]]] (
        	1
        	- Loop tan\[Beta] 8/3 g3^2 Conjugate[\[Mu]]/M3 H2[ msq[[i]]/M3, msd[[j]]/M3 ]  
        	+ Loop tan\[Beta] 3/2 g2^2 Conjugate[M2]/\[Mu] H2[M2/\[Mu], msq[[i]]/\[Mu]]
        	+ Loop tan\[Beta] gp^2 (
        	  1/9 Conjugate[\[Mu]]/M1 H2[msq[[i]]/M1, msd[[j]]/M1]
        	  + 1/6 Conjugate[M1]/\[Mu] H2[M1/\[Mu], msq[[i]]/\[Mu]]
        	  + 1/3 Conjugate[M1]/\[Mu] H2[M1/\[Mu], msd[[j]]/\[Mu]]
        	)
        	+ NonTan\[Beta] * Loop 12 g2^2 C00[\[Mu]/Q, msq[[i]]/\[Mu], M2/\[Mu]]
		+ NonTan\[Beta] * Loop gp^2 (
		  2/3 C00[\[Mu]/Q, msq[[i]]/\[Mu], M1/\[Mu]]
		  + 4/3 C00[\[Mu]/Q, msd[[j]]/\[Mu], M1/\[Mu]]
		)
      	)
	 - Cos[ArcTan[tan\[Beta]]] Loop tan\[Beta] 1/\[Mu] Sum[
        	  rotYu[[i, n]] 
        	  * ConjugateTranspose[rotTu][[n, m]] 
        	  * rotYd[[m, j]] H2[msu[[n]]/\[Mu], msq[[m]]/\[Mu]], 
        	  {n, 3}, {m, 3}
        	]
	 + NonTan\[Beta] rotTd[[i,j]] Cos[ArcTan[tan\[Beta]]] (
		+ Loop  8/3 g3^2 / M3 H2[msq[[i]]/M3, msd[[j]]/M3]
		- Loop  1/9 gp^2 / M1 H2[msq[[i]]/M1, msd[[j]]/M1]
	)
	 + NonTan\[Beta] Loop Cos[ArcTan[tan\[Beta]]] Sum[
        	  rotYu[[i, n]] 
        	  * ConjugateTranspose[rotYu][[n, m]] 
        	  * rotYd[[m, j]] H2[msu[[n]]/\[Mu], msq[[m]]/\[Mu]], 
        	  {n, 3}, {m, 3}
        	]
	+ NonTan\[Beta] Loop Cos[ArcTan[tan\[Beta]]] Sin[ArcTan[tan\[Beta]]]^2 2 Sum[
        	  rotYu[[i, n]] 
        	  * ConjugateTranspose[rotYu][[n, m]] 
        	  * rotYd[[m, j]] B1[mH/Q], 
        	  {n, 3}, {m, 3}
        	],
      {i, 3}, {j, 3}
  ];
  Transpose[Conjugate[Wq].Qwave[x].SMYd.Transpose[dwave[x]].Transpose[Wd]]*hwave[x]
];

OneLoopSMYe[x_] := Block[{
    (*Convert to SLHA2*)
    Ye = Transpose[FlavourBasis["Ye"] /. x],
    Te = Transpose[FlavourBasis["Te"] /. x], 
    SMYe, rotYe, rotTe, msl, mse, Wl, We,
    gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta]
  },
  {msl, Wl} = GetSoftMassMixing[FlavourBasis["m2L"] /. x];
  {mse, We} = GetSoftMassMixing[FlavourBasis["m2e"] /. x];
  rotYe = Transpose[Wl].Ye.Conjugate[We];
  rotTe = Transpose[Wl].Te.Conjugate[We];
  {gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta] , NonTan\[Beta]} = {"gp", "g2", "g3", "M1", "M2", "M3", "\[Mu]", "Q", "tan\[Beta]" , "NonTan\[Beta]"} /. x;
  NonTan\[Beta] = Boole[NonTan\[Beta]];
  SMYe=Table[
      rotYe[[i, j]] Cos[ArcTan[tan\[Beta]]] (
        1
        + Loop tan\[Beta] 3/2 g2^2 Conjugate[M2]/\[Mu] H2[M2/\[Mu], msl[[i]]/\[Mu]] 
        + Loop tan\[Beta] gp^2 (
          - Conjugate[\[Mu]]/M1 H2[msl[[i]]/M1, mse[[j]]/M1]
          - 1/2 Conjugate[M1]/\[Mu] H2[M1/\[Mu], msl[[i]]/\[Mu]]
          + Conjugate[M1]/\[Mu] H2[M1/\[Mu], mse[[j]]/\[Mu]]
        )
	+ NonTan\[Beta] Loop 12 g2^2 C00[\[Mu]/Q, msl[[i]]/\[Mu], M2/\[Mu]]
	+ NonTan\[Beta] Loop gp^2 (
	  -2 C00[\[Mu]/Q, msl[[i]]/\[Mu], M1/\[Mu]]
	  + 4 C00[\[Mu]/Q, mse[[j]]/\[Mu], M1/\[Mu]]
	)
      )
      + NonTan\[Beta] rotTe[[i, j]] Cos[ArcTan[tan\[Beta]]] Loop gp^2 / M1 H2[msl[[i]]/M1, mse[[j]]/M1],
      {i, 3}, {j, 3}
  ];
  Transpose[Conjugate[Wl].Lwave[x].SMYe.Transpose[ewave[x]].Transpose[We]]*hwave[x]
];

OneLoopSMYu[x_] := Block[{
    (* Convert to SLHA2 *)
    Yd = Transpose[FlavourBasis["Yd"] /. x],
    Td = Transpose[FlavourBasis["Td"] /. x],
    Yu = Transpose[FlavourBasis["Yu"] /. x],
    Tu = Transpose[FlavourBasis["Tu"] /. x],
    SMYu, rotYd, rotYu, rotTu, rotTd, msq, msu, msd, Wq, Wu, Wd,
    gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta], mH, U
  },
  {msq, Wq} = GetSoftMassMixing[FlavourBasis["m2Q"] /. x];
  {msu, Wu} = GetSoftMassMixing[FlavourBasis["m2u"] /. x];
  {msd, Wd} = GetSoftMassMixing[FlavourBasis["m2d"] /. x];
  rotYu = Transpose[Wq].Yu.Conjugate[Wu];
  rotTu = Transpose[Wq].Tu.Conjugate[Wu];
  rotYd = Transpose[Wq].Yd.Conjugate[Wd];
  rotTd = Transpose[Wq].Td.Conjugate[Wd];
  {mH, U} = GetHiggsDoubletMass[x];
  {gp, g2, g3, M1, M2, M3, \[Mu], Q, tan\[Beta], NonTan\[Beta]} = {"gp", "g2", "g3", "M1", "M2", "M3", "\[Mu]", "Q", "tan\[Beta]", "NonTan\[Beta]"} /. x;
  If[NonTan\[Beta],SMYu = Table[
	 rotYu[[i, j]] Sin[ArcTan[tan\[Beta]]] (
        	1
        	- Loop 1/tan\[Beta] 8/3 g3^2 Conjugate[\[Mu]]/M3 H2[ msq[[i]]/M3, msu[[j]]/M3]  
        	+ Loop 1/tan\[Beta] 3/2 g2^2 Conjugate[M2]/\[Mu] H2[M2/\[Mu], msq[[i]]/\[Mu]]
        	+ Loop 1/tan\[Beta] gp^2 (
        	  -2/9 Conjugate[\[Mu]]/M1 H2[msq[[i]]/M1, msu[[j]]/M1]
        	  - 1/6 Conjugate[M1]/\[Mu] H2[M1/\[Mu], msq[[i]]/\[Mu]]
        	  + 2/3 Conjugate[M1]/\[Mu] H2[M1/\[Mu], msu[[j]]/\[Mu]]
        	)
        	+ Loop 12 g2^2 C00[\[Mu]/Q, msq[[i]]/\[Mu], M2/\[Mu]]
		+ Loop gp^2 (
		  - 2/3 C00[\[Mu]/Q, msq[[i]]/\[Mu], M1/\[Mu]]
		  + 8/3 C00[\[Mu]/Q, msu[[j]]/\[Mu], M1/\[Mu]]
		)
      	)
	 - Sin[ArcTan[tan\[Beta]]] Loop 1/tan\[Beta] 1/\[Mu] Sum[
        	  rotYd[[i, n]] 
        	  * ConjugateTranspose[rotTd][[n, m]] 
        	  * rotYu[[m, j]] H2[msd[[n]]/\[Mu], msq[[m]]/\[Mu]], 
        	  {n, 3}, {m, 3}
        	]
	 + rotTu[[i,j]] Sin[ArcTan[tan\[Beta]]] (
		+ Loop  8/3 g3^2 / M3 H2[msq[[i]]/M3, msu[[j]]/M3]
		+ Loop  2/9 gp^2 / M1 H2[msq[[i]]/M1, msu[[j]]/M1]
	)
	 + Loop Sin[ArcTan[tan\[Beta]]] Sum[
        	  rotYd[[i, n]] 
        	  * ConjugateTranspose[rotYd][[n, m]] 
        	  * rotYu[[m, j]] H2[msd[[n]]/\[Mu], msq[[m]]/\[Mu]], 
        	  {n, 3}, {m, 3}
        	]
	+ Loop Sin[ArcTan[tan\[Beta]]] Cos[ArcTan[tan\[Beta]]]^2 2 Sum[
        	  rotYd[[i, n]] 
        	  * ConjugateTranspose[rotYd][[n, m]] 
        	  * rotYu[[m, j]] B1[mH/Q], 
        	  {n, 3}, {m, 3}
        	],
      {i, 3}, {j, 3}
    ],(*else NonTan\[Beta]==False*)
  SMYu = rotYu Sin[ArcTan[tan\[Beta]]]
  ];
  Transpose[Conjugate[Wq].Qwave[x].SMYu.Transpose[uwave[x]].Transpose[Wu]]*hwave[x]
];

OneLoopSM\[Kappa][x_] := Block[{
  \[Kappa], tan\[Beta], SM\[Kappa], rot\[Kappa], msl, Wl
  },
  {\[Kappa], tan\[Beta]} = {"\[Kappa]", "tan\[Beta]"} /. x;
  {msl, Wl} = GetSoftMassMixing[FlavourBasis["m2L"] /. x];
  \[Kappa] = \[Kappa] * 10^16; (*for numeric stability change units*)
  rot\[Kappa] = Transpose[Wl].\[Kappa].Wl;
  SM\[Kappa] = rot\[Kappa] * Sin[ArcTan[tan\[Beta]]]^2;
  Conjugate[Wl].Lwave[x].SM\[Kappa].Transpose[Lwave[x]].ConjugateTranspose[Wl]*(hwave[x]^2)*10^(-16)(*restore original units*)
];

End[];

EndPackage[];
