(* ::Package:: *)

(*
run
math -initfile dim6top.m
having set
$Path = Join[{ToFileName[$HomeDirectory, ".../feynrules-current"]}, $Path]
$FeynRulesPath = ToFileName[$HomeDirectory, ".../feynrules-current"]
in ~/.Mathematica/Kernel/init.m
*)


Print["--- Set Feynman gauge"]
FeynmanGauge = True;

Print["--- Load FeynRules"];
FR$Parallel = True;
<<FeynRules`

Print["--- Load models"];
(*LoadModel[ ToFileName[$FeynRulesPath,"Models/SM/SM.fr"],"dim6top_each_coupling_order.fr"];*)
(**)LoadModel[ ToFileName[$FeynRulesPath,"Models/SM/SM.fr"],"dim6top.fr"];(**)

M$Information = {
	Authors      -> {"G. Durieux", "C. Zhang"}, 
	Version      -> DateString[],
	Date         -> DateString[],
	Institutions -> {},
	Emails       -> {},
	URLs         -> {"http://arxiv.org/abs/1802.07237", "https://feynrules.irmp.ucl.ac.be/wiki/dim6top"}
};

Print["--- Load restriction"];
LoadRestriction[ToFileName[$FeynRulesPath,"Models/SM/DiagonalCKM.rst"]]
Print["--- Load restriction"];
LoadRestriction[ToFileName[$FeynRulesPath,"Models/SM/Massless.rst"]]
(*Print["--- Load restriction"];
LoadRestriction["dim6top.rst"];*)

(*
Print["--- Get vertices"]
verttot =<<saved_verttot ;

Print["--- Write UFO model with those"]
WriteUFO[
	{},Input->verttot,
	AddDecays->False,
	(*Output->"dim6top_LO_UFO"*)
	Output->"dim6top_LO_UFO_each_coupling_order"
];
Exit[];
*)


Print["--- Define operators"];
(* Sigma[mu_,nu_, sp1_,sp2_] = I/2 (Ga[mu,sp1,sp5] Ga[nu,sp5,sp2] - Ga[nu,sp1,sp5] Ga[mu,sp5,sp2]); *)
   Sigma[mu_,nu_, sp1_,sp2_] = I   (Ga[mu,sp1,sp5] Ga[nu,sp5,sp2] - ME[mu,nu] IndexDelta[sp1,sp2] );

(** 2q Hermitian **)
Ophiq1  = Cphiq1[f1,f2]  I (Phibar[k1]DC[Phi[k1],mu]-DC[Phibar[k1],mu]Phi[k1])               (QLbar[sp1,k3,f1,c1].QL[sp2,k3,f2,c1] Ga[mu,sp1,sp2]) ;
Ophiq3  = Cphiq3[f1,f2]  I (Phibar[k1]DC[Phi[k2],mu]-DC[Phibar[k1],mu]Phi[k2]) PS[j1,k1,k2]  (QLbar[sp1,k3,f1,c1].QL[sp2,k4,f2,c1] Ga[mu,sp1,sp2]) PS[j1,k3,k4] ;
Ophiu   = Cphiu [f1,f2]  I (Phibar[k1]DC[Phi[k1],mu]-DC[Phibar[k1],mu]Phi[k1])               (uRbar[sp1,   f1,c1].uR[sp2,   f2,c1] Ga[mu,sp1,sp2]) ;
Ophid   = Cphid [f1,f2]  I (Phibar[k1]DC[Phi[k1],mu]-DC[Phibar[k1],mu]Phi[k1])               (dRbar[sp1,   f1,c1].dR[sp2,   f2,c1] Ga[mu,sp1,sp2]) ;
(** 2q non-Hermitian **)
Ophiud = Cphiud[f1,f2] I Phi[k1] DC[Phi[k2],mu]        Eps[k1,k2]             (uRbar[sp1,   f1,c1].dR[sp2,f2,c1]) Ga[mu,sp1,sp2] ;
Ouphi  = Cuphi [f1,f2] (Phibar[k3] Phi[k3]-vev^2/2)  Phibar[k2] Eps[k1,k2]  (QLbar[sp1,k1,f1,c1].uR[sp1,f2,c1]) ;
OdW    = CdW   [f1,f2] FS[Wi,mu,nu,j1]    Phi[k2]              PS[j1,k1,k2] (QLbar[sp1,k1,f1,c1].dR[sp2,f2,c1]) Sigma[mu,nu,sp1,sp2] ;
OuW    = CuW   [f1,f2] FS[Wi,mu,nu,j1]    Phibar[k2] Eps[k3,k2]PS[j1,k1,k3] (QLbar[sp1,k1,f1,c1].uR[sp2,f2,c1]) Sigma[mu,nu,sp1,sp2] ;
OuB    = CuB   [f1,f2] FS[B ,mu,nu   ]    Phibar[k2] Eps[k1,k2]             (QLbar[sp1,k1,f1,c1].uR[sp2,f2,c1]) Sigma[mu,nu,sp1,sp2] ;
OuG    = CuG   [f1,f2] FS[G ,mu,nu,a1]    Phibar[k2] Eps[k1,k2]             (QLbar[sp1,k1,f1,c1].uR[sp2,f2,c2]) Sigma[mu,nu,sp1,sp2] T[a1,c1,c2];
(* 2l2t *)
Olq1 = Clq1[f3,f4,f1,f2] (QLbar[sp1,k1,f1,c1].QL[sp2,k1,f2,c1] Ga[mu,sp1,sp2]) (LLbar[sp3,k3,f3].LL[sp4,k3,f4] Ga[mu,sp3,sp4]) ;
Olq3 = Clq3[f3,f4,f1,f2] (QLbar[sp1,k1,f1,c1].QL[sp2,k2,f2,c1] Ga[mu,sp1,sp2]) (LLbar[sp3,k3,f3].LL[sp4,k4,f4] Ga[mu,sp3,sp4])  PS[j1,k1,k2] PS[j1,k3,k4] ;
Olu  = Clu [f3,f4,f1,f2] (uRbar[sp1,   f1,c1].uR[sp2,   f2,c1] Ga[mu,sp1,sp2]) (LLbar[sp3,k3,f3].LL[sp4,k3,f4] Ga[mu,sp3,sp4]) ;
Oeq  = Ceq [f3,f4,f1,f2] (QLbar[sp1,k1,f1,c1].QL[sp2,k1,f2,c1] Ga[mu,sp1,sp2]) (lRbar[sp3,   f3].lR[sp4,   f4] Ga[mu,sp3,sp4]) ;
Oeu  = Ceu [f3,f4,f1,f2] (uRbar[sp1,   f1,c1].uR[sp2,   f2,c1] Ga[mu,sp1,sp2]) (lRbar[sp3,   f3].lR[sp4,   f4] Ga[mu,sp3,sp4]) ;
(* 2l2q scalar and tensor *)
Olequ1 = Clequ1[f3,f4,f1,f2]  (LLbar[sp3,k3,f3].lR[sp3,f4]) (QLbar[sp1,k1,f1,c1].uR[sp1,f2,c1]) Eps[k3,k1] ;
Olequ3 = Clequ3[f3,f4,f1,f2]  (LLbar[sp3,k3,f3].lR[sp4,f4]) (QLbar[sp1,k1,f1,c1].uR[sp2,f2,c1]) Eps[k3,k1] Sigma[mu,nu,sp1,sp2] Sigma[mu,nu,sp3,sp4] ;
Oledq  = Cledq[f3,f4,f1,f2]   (LLbar[sp3,k2,f3].lR[sp3,f4]) (dRbar[sp1,f1,c1].QL[sp1,k2,f2,c1]) ;

(* 4q *)
Oqq83 = Cqq83[f5,f6,f1,f2] (QLbar[sp1,k1,f1,c1].QL[sp2,k2,f2,c2] Ga[mu,sp1,sp2] T[a,c1,c2] PS[j1,k1,k2]) (QLbar[sp3,k3,f5,c3].QL[sp4,k4,f6,c4] Ga[mu,sp3,sp4] T[a,c3,c4] PS[j1,k3,k4]);
Oqq81 = Cqq81[f5,f6,f1,f2] (QLbar[sp1,k1,f1,c1].QL[sp2,k1,f2,c2] Ga[mu,sp1,sp2] T[a,c1,c2]) (QLbar[sp3,k3,f5,c3].QL[sp4,k3,f6,c4] Ga[mu,sp3,sp4] T[a,c3,c4]);
Oqu8  = Cqu8 [f5,f6,f1,f2] (QLbar[sp1,k1,f1,c1].QL[sp2,k1,f2,c2] Ga[mu,sp1,sp2] T[a,c1,c2]) (uRbar[sp3,   f5,c3].uR[sp4,   f6,c4] Ga[mu,sp3,sp4] T[a,c3,c4]);
Oqd8  = Cqd8 [f5,f6,f1,f2] (QLbar[sp1,k1,f1,c1].QL[sp2,k1,f2,c2] Ga[mu,sp1,sp2] T[a,c1,c2]) (dRbar[sp3,   f5,c3].dR[sp4,   f6,c4] Ga[mu,sp3,sp4] T[a,c3,c4]);
Ouu8  = Cuu8 [f5,f6,f1,f2] (uRbar[sp1,   f1,c1].uR[sp2,   f2,c2] Ga[mu,sp1,sp2] T[a,c1,c2]) (uRbar[sp3,   f5,c3].uR[sp4,   f6,c4] Ga[mu,sp3,sp4] T[a,c3,c4]);
Oud8  = Cud8 [f5,f6,f1,f2] (uRbar[sp1,   f1,c1].uR[sp2,   f2,c2] Ga[mu,sp1,sp2] T[a,c1,c2]) (dRbar[sp3,   f5,c3].dR[sp4,   f6,c4] Ga[mu,sp3,sp4] T[a,c3,c4]);
(* Ouq8  = Cuq8 [f5,f6,f1,f2] (uRbar[sp1,   f1,c1].uR[sp2,   f2,c2] Ga[mu,sp1,sp2] T[a,c1,c2]) (QLbar[sp3,k3,f5,c3].QL[sp4,k3,f6,c4] Ga[mu,sp3,sp4] T[a,c3,c4]); *)

Oqq13 = Cqq13[f5,f6,f1,f2] (QLbar[sp1,k1,f1,c1].QL[sp2,k2,f2,c1] Ga[mu,sp1,sp2] PS[j1,k1,k2]) (QLbar[sp3,k3,f5,c3].QL[sp4,k4,f6,c3] Ga[mu,sp3,sp4] PS[j1,k3,k4]);
Oqq11 = Cqq11[f5,f6,f1,f2] (QLbar[sp1,k1,f1,c1].QL[sp2,k1,f2,c1] Ga[mu,sp1,sp2]) (QLbar[sp3,k3,f5,c3].QL[sp4,k3,f6,c3] Ga[mu,sp3,sp4]);
Oqu1  = Cqu1 [f5,f6,f1,f2] (QLbar[sp1,k1,f1,c1].QL[sp2,k1,f2,c1] Ga[mu,sp1,sp2]) (uRbar[sp3,   f5,c3].uR[sp4,   f6,c3] Ga[mu,sp3,sp4]);
Oqd1  = Cqd1 [f5,f6,f1,f2] (QLbar[sp1,k1,f1,c1].QL[sp2,k1,f2,c1] Ga[mu,sp1,sp2]) (dRbar[sp3,   f5,c3].dR[sp4,   f6,c3] Ga[mu,sp3,sp4]);
Ouu1  = Cuu1 [f5,f6,f1,f2] (uRbar[sp1,   f1,c1].uR[sp2,   f2,c1] Ga[mu,sp1,sp2]) (uRbar[sp3,   f5,c3].uR[sp4,   f6,c3] Ga[mu,sp3,sp4]);
Oud1  = Cud1 [f5,f6,f1,f2] (uRbar[sp1,   f1,c1].uR[sp2,   f2,c1] Ga[mu,sp1,sp2]) (dRbar[sp3,   f5,c3].dR[sp4,   f6,c3] Ga[mu,sp3,sp4]);
(* Ouq1  = Cuq1 [f5,f6,f1,f2] (uRbar[sp1,   f1,c1].uR[sp2,   f2,c1] Ga[mu,sp1,sp2]) (QLbar[sp3,k3,f5,c3].QL[sp4,k3,f6,c3] Ga[mu,sp3,sp4]); *)

Oquqd1  = Cquqd1[f5,f6,f1,f2]  (QLbar[sp1,k1,f1,c1].uR[sp1,   f2,c1]) (QLbar[sp3,k3,f5,c3].dR[sp3,   f6,c3]) Eps[k1,k3];
Oquqd8  = Cquqd8[f5,f6,f1,f2]  (QLbar[sp1,k1,f1,c1].uR[sp1,   f2,c2]) (QLbar[sp3,k3,f5,c3].dR[sp3,   f6,c4]) Eps[k1,k3] T[a,c1,c2] T[a,c3,c4];
Oquqd1T = Cquqd1T[f5,f6,f1,f2] (QLbar[sp1,k1,f1,c1].uR[sp2,   f2,c1]) (QLbar[sp3,k3,f5,c3].dR[sp4,   f6,c3]) Eps[k1,k3] Sigma[mu,nu,sp1,sp2] Sigma[mu,nu,sp3,sp4];
Oquqd8T = Cquqd8T[f5,f6,f1,f2] (QLbar[sp1,k1,f1,c1].uR[sp2,   f2,c2]) (QLbar[sp3,k3,f5,c3].dR[sp4,   f6,c4]) Eps[k1,k3] T[a,c1,c2] T[a,c3,c4] Sigma[mu,nu,sp1,sp2] Sigma[mu,nu,sp3,sp4];

Ouqqu1 = Cuqqu1[f5,f6,f1,f2] (uRbar[sp1,   f1,c1].QL[sp1,k2,f2,c1]) (QLbar[sp3,k2,f5,c3].uR[sp3,   f6,c3]);
Ouqqu8 = Cuqqu8[f5,f6,f1,f2] (uRbar[sp1,   f1,c1].QL[sp1,k2,f2,c2]) (QLbar[sp3,k2,f5,c3].uR[sp3,   f6,c4]) T[a,c1,c2] T[a,c3,c4];
Odqqd1 = Cdqqd1[f5,f6,f1,f2] (dRbar[sp1,   f1,c1].QL[sp1,k2,f2,c1]) (QLbar[sp3,k2,f5,c3].dR[sp3,   f6,c3]);
Odqqd8 = Cdqqd8[f5,f6,f1,f2] (dRbar[sp1,   f1,c1].QL[sp1,k2,f2,c2]) (QLbar[sp3,k2,f5,c3].dR[sp3,   f6,c4]) T[a,c1,c2] T[a,c3,c4];
Oduud1 = Cduud1[f5,f6,f1,f2] (dRbar[sp1,   f1,c1].uR[sp2,   f2,c1] Ga[mu,sp1,sp2]) (uRbar[sp3,   f5,c3].dR[sp4,   f6,c3] Ga[mu,sp3,sp4]);
Oduud8 = Cduud8[f5,f6,f1,f2] (dRbar[sp1,   f1,c1].uR[sp2,   f2,c2] Ga[mu,sp1,sp2]) (uRbar[sp3,   f5,c3].dR[sp4,   f6,c4] Ga[mu,sp3,sp4]) T[a,c1,c2] T[a,c3,c4];

FierzPauli = { PS[j1_,k1_,k2_]*PS[j1_,k3_,k4_] -> 2*(IndexDelta[k1,k4]*IndexDelta[k2,k3] - 1/2*IndexDelta[k1,k2]*IndexDelta[k3,k4]) };
DefinePS   = { PauliSigma[j1_,k1_,k2_] -> PS[j1,k1,k2] };
Gen33      = { f1|f2 -> 3 };
Gen31      = { f1 -> 3, f2 -> 1 };
Gen32      = { f1 -> 3, f2 -> 2 };
Gen13      = { f1 -> 1, f2 -> 3 };
Gen23      = { f1 -> 2, f2 -> 3 };
Gen3311    = { f1|f2 -> 3, f5|f6 -> 1 };
Gen3322    = { f1|f2 -> 3, f5|f6 -> 2 };
Gen1133    = { f1|f2 -> 1, f5|f6 -> 3 };
Gen2233    = { f1|f2 -> 2, f5|f6 -> 3 };

Gen3111    = { f1 -> 3, f2 -> 1, f5|f6 -> 1 };
Gen3122    = { f1 -> 3, f2 -> 1, f5|f6 -> 2 };
Gen3211    = { f1 -> 3, f2 -> 2, f5|f6 -> 1 };
Gen3222    = { f1 -> 3, f2 -> 2, f5|f6 -> 2 };
Gen1131    = { f1|f2 -> 1, f5 -> 3, f6 -> 1 };
Gen2231    = { f1|f2 -> 2, f5 -> 3, f6 -> 1 };
Gen1132    = { f1|f2 -> 1, f5 -> 3, f6 -> 2 };
Gen2232    = { f1|f2 -> 2, f5 -> 3, f6 -> 2 };

Gen3331    = { f1|f2|f5 -> 3, f6 -> 1};
Gen3332    = { f1|f2|f5 -> 3, f6 -> 2};
Gen3313    = { f1|f2|f6 -> 3, f5 -> 1};
Gen3323    = { f1|f2|f6 -> 3, f5 -> 2};
Gen3133    = { f1|f5|f6 -> 3, f2 -> 1};
Gen3233    = { f1|f5|f6 -> 3, f2 -> 2};
Gen1333    = { f2|f5|f6 -> 3, f1 -> 1};
Gen2333    = { f2|f5|f6 -> 3, f1 -> 2};

Gen3333    = { f1|f2|f5|f6 -> 3};

Electron     = { f3|f4 -> 1 };
Muon         = { f3|f4 -> 2 };
Tau          = { f3|f4 -> 3 };
RemovePart = { GP|GPbar|G0 -> 0};
ckm        = { Conjugate[CKM[i_,j_]] -> IndexDelta[i,j], CKM[i_,j_] -> IndexDelta[i,j] };
Print["--- !!!! Goldstones will be removed"]
Print["--- !!!! CKM will be set to I"]


Print["--- Expand indices"]
(*FR$Parallel=True;*)

O2fHerm    = {Ophiq1, Ophiq3, Ophiu, Ophid};
O2fNonHerm = {Ophiud, Ouphi, OdW, OuW, OuB, OuG};
O4fHerm    = {Olq1, Olq3, Olu, Oeq, Oeu};
O4fNonHerm = {Olequ1, Olequ3, Oledq};

(* Expand two-fermion *)
Print["--- Expand two fermion operators"]
L2fHerm=ParallelMap[
	ExpandIndices[
	ExpandIndices[
		((#)		/. FierzPauli ) /Lambda^2 
	,FlavorExpand->SU2W]	/. DefinePS
	,FlavorExpand->SU2D]    /. RemovePart
&, O2fHerm ] ;
L2fNonHerm=ParallelMap[
	ExpandIndices[
	ExpandIndices[
		((#)		/. FierzPauli ) /Lambda^2 
	,FlavorExpand->SU2W]	/. DefinePS
	,FlavorExpand->SU2D]    /. RemovePart
&, O2fNonHerm ] ;


Print["--- Expand 2f FCNC"]
L2fFCNC = Join[
	ParallelMap[ExpandIndices[ (# /. Gen31) + (# /. Gen32)                              , FlavorExpand->Generation] /. ckm &, L2fHerm ],
	ParallelMap[ExpandIndices[ (# /. Gen31) + (# /. Gen32) + (# /. Gen13) + (# /. Gen23), FlavorExpand->Generation] /. ckm &, L2fNonHerm ] 
	];


L2fFCNC = ParallelMap[ (# + HC[#]) &, L2fFCNC ] ;
Print["--- Expand 2f (flavour conserving)"]
L2fHerm    = ParallelMap[ExpandIndices[  # /. Gen33, FlavorExpand->Generation] /. ckm &, L2fHerm ] ;
L2fNonHerm = ParallelMap[ExpandIndices[  # /. Gen33, FlavorExpand->Generation] /. ckm &, L2fNonHerm ] ;
L2f = Join[L2fHerm, ParallelMap[ (# + HC[#]) &, L2fNonHerm ]];


(* Expand four-fermion 2q2l *)
Print["--- Expand four-fermion operators"]
L4fHerm=ParallelMap[
	ExpandIndices[
	ExpandIndices[
		(# /Lambda^2 )	/. FierzPauli 
	,FlavorExpand->SU2W]	/. DefinePS
	,FlavorExpand->SU2D]    /. RemovePart
&, O4fHerm ] ;
L4fNonHerm=ParallelMap[
	ExpandIndices[
	ExpandIndices[
		(# /Lambda^2 )	/. FierzPauli 
	,FlavorExpand->SU2W]	/. DefinePS
	,FlavorExpand->SU2D]    /. RemovePart
&, O4fNonHerm ] ;


Print["--- Expand 2q2l FCNC"]
L4fFCNC = Join[
	ParallelMap[ExpandIndices[
		   (( (# /. Electron ) + (# /. Muon ) + (# /. Tau ) ) /. Gen31)
		 + (( (# /. Electron ) + (# /. Muon ) + (# /. Tau ) ) /. Gen32)
		 , FlavorExpand->Generation] /. ckm &, L4fHerm ],
	ParallelMap[ExpandIndices[
		   (( (# /. Electron ) + (# /. Muon ) + (# /. Tau ) ) /. Gen31)
		 + (( (# /. Electron ) + (# /. Muon ) + (# /. Tau ) ) /. Gen32)
		 + (( (# /. Electron ) + (# /. Muon ) + (# /. Tau ) ) /. Gen13)
		 + (( (# /. Electron ) + (# /. Muon ) + (# /. Tau ) ) /. Gen23)
		 , FlavorExpand->Generation] /. ckm &, L4fNonHerm ] 
	];
L4fFCNC = ParallelMap[ (# + HC[#]) &, L4fFCNC ] ;


Print["--- Expand 2q2l (flavour conserving)"]
L4fHerm    = ParallelMap[ExpandIndices[ ((# /. Electron ) + (# /. Muon ) + (# /. Tau )  /. Gen33), FlavorExpand->Generation] /. ckm &, L4fHerm ] ;
L4fNonHerm = ParallelMap[ExpandIndices[ ((# /. Electron ) + (# /. Muon ) + (# /. Tau )  /. Gen33), FlavorExpand->Generation] /. ckm &, L4fNonHerm ] ;
L4f = Join[ L4fHerm, ParallelMap[ (# + HC[#]) &, L4fNonHerm ]] ;


(* Expand 4q *)
Print["--- Expand 4q operators"]
O4q = {Oqq83, Oqq81, Oqu8, Oqd8, Ouu8, Oud8,
       Oqq13, Oqq11, Oqu1, Oqd1, Ouu1, Oud1,
       Oquqd1, Oquqd8, Oquqd1T, Oquqd8T,
       Ouqqu1, Ouqqu8, Odqqd1, Odqqd8,
       Oduud1, Oduud8 };
L4q = ParallelMap[
	ExpandIndices[
	ExpandIndices[
		((#) /Lambda^2 )/. FierzPauli 
	,FlavorExpand->SU2W]	/. DefinePS
	,FlavorExpand->SU2D]    /. RemovePart
&, O4q ] ;
Oqq83  = L4q[[ 1]];
Oqq81  = L4q[[ 2]];
Oqu8   = L4q[[ 3]];
Oqd8   = L4q[[ 4]];
Ouu8   = L4q[[ 5]];
Oud8   = L4q[[ 6]];
Oqq13  = L4q[[ 7]];
Oqq11  = L4q[[ 8]];
Oqu1   = L4q[[ 9]];
Oqd1   = L4q[[10]];
Ouu1   = L4q[[11]];
Oud1   = L4q[[12]];
Oquqd1 = L4q[[13]];
Oquqd8 = L4q[[14]];
Oquqd1T= L4q[[15]];
Oquqd8T= L4q[[16]];
Ouqqu1 = L4q[[17]];
Ouqqu8 = L4q[[18]];
Odqqd1 = L4q[[19]];
Odqqd8 = L4q[[20]];
Oduud1 = L4q[[21]];
Oduud8 = L4q[[22]];



Print["--- Expand qqtt"]
Lqqtt = Join[ ParallelMap[ ExpandIndices[ (# /. Gen3311) + (# /. Gen3322) , FlavorExpand->Generation]	/. ckm &,
		{Oqq83, Oqq81, Oqu8, Oqd8, Ouu8, Oud8,
		 Oqq13, Oqq11, Oqu1, Oqd1, Ouu1, Oud1} ],
	      ParallelMap[ ExpandIndices[ (# /. Gen1133) + (# /. Gen2233) , FlavorExpand->Generation]	/. ckm &,
		{Oqu8, Oqu1} ]
	] ;


Print["--- Expand tttt"]
Ltttt = ParallelMap[ ExpandIndices[ (# /. Gen3333)]	/. ckm &, 
	{Oqq13,  Oqq11,  Oqu1, Oqd1, Ouu1,  Oud1,
	                  Oqu8, Oqd8,        Oud8,
	 Oquqd1, Oquqd8}] ;
Ltttt = Join[ Ltttt, Table[ HC[Ltttt[[i]]],{i,Length[Ltttt]-1,Length[Ltttt]}] ];


Print["--- Expand qqttS"]
LqqttS= Join[ ParallelMap[ ExpandIndices[ (# /. Gen3311) + (# /. Gen3322) ]	/. ckm &,
		{Ouqqu1, Ouqqu8,
		 Odqqd1, Odqqd8,
		 Oquqd1, Oquqd8, Oquqd1T, Oquqd8T,
		 Oduud1, Oduud8} ],
	      ParallelMap[ ExpandIndices[ (# /. Gen1133) + (# /. Gen2233) ]	/. ckm &,
		{Oquqd1, Oquqd8, Oquqd1T, Oquqd8T} ]             (* Should there be a sign flip here from epsilon? *)
	] ;
LqqttS = ParallelMap[ (# + HC[#]) &, LqqttS];


Print["--- Expand 3h1l FCNC"]
L3h1lFCNC = Join[ ParallelMap[ ExpandIndices[ (# /. Gen3331) + (# /. Gen3332) ]	/. ckm &,
			{Oqq13, Oqq11, Ouu1,
			 Oqu1, Oqu8, Oqd1, Oqd8, Oud1, Oud8,
			 Oquqd1, Oquqd8} ],
		  ParallelMap[ ExpandIndices[ (# /. Gen3133) + (# /. Gen3233) ]	/. ckm &,
			{Oqu1, Oqu8, Oqd1, Oqd8, Oud1, Oud8,
			 Oquqd1, Oquqd8} ],
		  ParallelMap[ ExpandIndices[ (# /. Gen3313) + (# /. Gen3323) ]	/. ckm &,
			{Oquqd1, Oquqd8} ],
		  ParallelMap[ ExpandIndices[ (# /. Gen1333) + (# /. Gen2333) ]	/. ckm &,
			{Oquqd1, Oquqd8} ]
	] ;
L3h1lFCNC = ParallelMap[ (# + HC[#]) &, L3h1lFCNC];


Print["--- Expand 1l3h FCNC"]
L1h3lFCNC = Join[ ParallelMap[ ExpandIndices[ 
			 (# /. Gen3111) + (# /. Gen3122)
			+(# /. Gen3211) + (# /. Gen3222) ]	/. ckm &,
			{Oqq11, Oqq13, Oqq81, Oqq83, Ouu1, Ouu8, Oud1, Oud8,  Oqd1, Oqd8,
			 Oqu1, Oqu8} ],
		  ParallelMap[ ExpandIndices[ 
			 (# /. Gen1131) + (# /. Gen2231)
			+(# /. Gen1132) + (# /. Gen2232) ]	/. ckm &,
			{Oqu1, Oqu8} ]
	] ;
L1h3lFCNC = ParallelMap[ (# + HC[#]) &, L1h3lFCNC];


Print["--- Total Lagrangian"]
(* Add hermitian conjugate to form Ldim6tot *)
Ldim6tot = Join[ L2f, L4f, Lqqtt, Ltttt, LqqttS, L2fFCNC, L4fFCNC, L3h1lFCNC, L1h3lFCNC ];
Ldim6totSum = Total[ Ldim6tot ];


Print["--- Save Lagrangian"]
Ldim6tot >> saved_Ldim6tot ;


Lsm = ExpandIndices[
ExpandIndices[
	ExpandIndices[
		(LSM		/. FierzPauli )
	,FlavorExpand->SU2W]	/. DefinePS
	,FlavorExpand->SU2D]    /. RemovePart /. ckm ];


Print["--- Define total Lagrangian and rename operator coefficients"]
Ltot = Lsm + Total[Ldim6tot] /. {
	Cuphi[3,3]   -> ctp + I ctpI,
	Cphiq1[3,3]  -> cpQM + cpQ3,
	Cphiq3[3,3]  -> cpQ3,
	Cphiu[3,3]   -> cpt,
	Cphid[3,3]   -> cpb,
	Cphiud[3,3]  -> cptb + I cptbI,
	CuB[3,3]     -> (cw*ctW - ctZ)/sw + I (cw*ctWI - ctZI)/sw,
	CuW[3,3]     -> ctW + I ctWI,
	CdW[3,3]     -> cbW + I cbWI,
	CuG[3,3]     -> (ctG + I ctGI)*Gstrong,
	
	Clq1[1,1,3,3]   -> cQlM1 + cQl31,
	Clq1[2,2,3,3]   -> cQlM2 + cQl32,
	Clq1[3,3,3,3]   -> cQlM3 + cQl33,
	Clq3[1,1,3,3]   -> cQl31,
	Clq3[2,2,3,3]   -> cQl32,
	Clq3[3,3,3,3]   -> cQl33,
	Ceq[1,1,3,3]    -> cQe1,
	Ceq[2,2,3,3]    -> cQe2,
	Ceq[3,3,3,3]    -> cQe3,
	Clu[1,1,3,3]    -> ctl1,
	Clu[2,2,3,3]    -> ctl2,
	Clu[3,3,3,3]    -> ctl3,
	Ceu[1,1,3,3]    -> cte1,
	Ceu[2,2,3,3]    -> cte2,
	Ceu[3,3,3,3]    -> cte3,
	Clequ1[1,1,3,3] -> ctlS1 + I ctlSI1,
	Clequ1[2,2,3,3] -> ctlS2 + I ctlSI2,
	Clequ1[3,3,3,3] -> ctlS3 + I ctlSI3,
	Clequ3[1,1,3,3] -> ctlT1 + I ctlTI1,
	Clequ3[2,2,3,3] -> ctlT2 + I ctlTI2,
	Clequ3[3,3,3,3] -> ctlT3 + I ctlTI3,
	Cledq[1,1,3,3]  -> cblS1 + I cblSI1,
	Cledq[2,2,3,3]  -> cblS2 + I cblSI2,
	Cledq[3,3,3,3]  -> cblS3 + I cblSI3,

	Cqq83[1,1,3,3] -> cQq83,
	Cqq81[1,1,3,3] -> cQq81,
	Cqu8 [1,1,3,3] -> cQu8,
	Cqu8 [3,3,1,1] -> ctq8,
	Cqd8 [1,1,3,3] -> cQd8,
	Cuu8 [1,1,3,3] -> ctu8,
	Cud8 [1,1,3,3] -> ctd8,
	Cqq13[1,1,3,3] -> cQq13,
	Cqq11[1,1,3,3] -> cQq11,
	Cqu1 [1,1,3,3] -> cQu1,
	Cqu1 [3,3,1,1] -> ctq1,
	Cqd1 [1,1,3,3] -> cQd1,
	Cuu1 [1,1,3,3] -> ctu1,
	Cud1 [1,1,3,3] -> ctd1,

	Cqq83[2,2,3,3] -> cQq83,
	Cqq81[2,2,3,3] -> cQq81,
	Cqu8 [2,2,3,3] -> cQu8,
	Cqu8 [3,3,2,2] -> ctq8,
	Cqd8 [2,2,3,3] -> cQd8,
	Cuu8 [2,2,3,3] -> ctu8,
	Cud8 [2,2,3,3] -> ctd8,
	Cqq13[2,2,3,3] -> cQq13,
	Cqq11[2,2,3,3] -> cQq11,
	Cqu1 [2,2,3,3] -> cQu1,
	Cqu1 [3,3,2,2] -> ctq1,
	Cqd1 [2,2,3,3] -> cQd1,
	Cuu1 [2,2,3,3] -> ctu1,
	Cud1 [2,2,3,3] -> ctd1,

	Cqq13[3,3,3,3] -> cQQ8/8,
	Cqq11[3,3,3,3] -> cQQ1/2 + cQQ8/24,
	Cqu1 [3,3,3,3] -> cQt1,
	Cqd1 [3,3,3,3] -> cQb1,
	Cuu1 [3,3,3,3] -> ctt1,
	Cud1 [3,3,3,3] -> ctb1,
	Cqu8 [3,3,3,3] -> cQt8,
	Cqd8 [3,3,3,3] -> cQb8,
	Cud8 [3,3,3,3] -> ctb8,
	
	Cquqd1[3,3,3,3] -> cQtQb1 + I cQtQb1I,
	Cquqd8[3,3,3,3] -> cQtQb8 + I cQtQb8I,

	Cuqqu1[1,1,3,3] -> ctQqu1 + I*ctQqu1I,
	Cuqqu8[1,1,3,3] -> ctQqu8 + I*ctQqu8I,
	Cdqqd1[1,1,3,3] -> cbQqd1 + I*cbQqd1I,
	Cdqqd8[1,1,3,3] -> cbQqd8 + I*cbQqd8I,
	Cquqd1[1,1,3,3] -> cQtqd1 + I*cQtqd1I,
	Cquqd8[1,1,3,3] -> cQtqd8 + I*cQtqd8I,
	Cquqd1[3,3,1,1] -> cQbqu1 + I*cQbqu1I,
	Cquqd8[3,3,1,1] -> cQbqu8 + I*cQbqu8I,
	Cquqd1T[1,1,3,3] -> cQtqd1T + I*cQtqd1TI,
	Cquqd8T[1,1,3,3] -> cQtqd8T + I*cQtqd8TI,
	Cquqd1T[3,3,1,1] -> cQbqu1T + I*cQbqu1TI,
	Cquqd8T[3,3,1,1] -> cQbqu8T + I*cQbqu8TI,
	Cduud1[1,1,3,3] -> cbtud1 + I*cbtud1I,
	Cduud8[1,1,3,3] -> cbtud8 + I*cbtud8I,

	Cuqqu1[2,2,3,3] -> ctQqu1 + I*ctQqu1I,
	Cuqqu8[2,2,3,3] -> ctQqu8 + I*ctQqu8I,
	Cdqqd1[2,2,3,3] -> cbQqd1 + I*cbQqd1I,
	Cdqqd8[2,2,3,3] -> cbQqd8 + I*cbQqd8I,
	Cquqd1[2,2,3,3] -> cQtqd1 + I*cQtqd1I,
	Cquqd8[2,2,3,3] -> cQtqd8 + I*cQtqd8I,
	Cquqd1[3,3,2,2] -> cQbqu1 + I*cQbqu1I,
	Cquqd8[3,3,2,2] -> cQbqu8 + I*cQbqu8I,
	Cquqd1T[2,2,3,3] -> cQtqd1T + I*cQtqd1TI,
	Cquqd8T[2,2,3,3] -> cQtqd8T + I*cQtqd8TI,
	Cquqd1T[3,3,2,2] -> cQbqu1T + I*cQbqu1TI,
	Cquqd8T[3,3,2,2] -> cQbqu8T + I*cQbqu8TI,
	Cduud1[2,2,3,3] -> cbtud1 + I*cbtud1I,
	Cduud8[2,2,3,3] -> cbtud8 + I*cbtud8I,
	
	(* FCNC up *)
	Cuphi [3,1]  -> ctpx31 + I ctpIx31,
	Cphiud[3,1]  -> cptbx31 + I cptbIx31,
	CuB   [3,1]  -> (- sw*ctZx31 + cw*ctAx31) + I*(- sw*ctZIx31 + cw*ctAIx31),
	CuW   [3,1]  -> (+ cw*ctZx31 + sw*ctAx31) + I*(+ cw*ctZIx31 + sw*ctAIx31),
	CdW   [3,1]  -> cbWx31 + I cbWIx31,
	CuG   [3,1]  -> (ctGx31 + I ctGIx31)*Gstrong,
	
	Cuphi [1,3]  -> ctpx13 + I ctpIx13,
	Cphiud[1,3]  -> cptbx13 + I cptbIx13,
	CuB   [1,3]  -> (- sw*ctZx13 + cw*ctAx13) + I*(- sw*ctZIx13 + cw*ctAIx13),
	CuW   [1,3]  -> (+ cw*ctZx13 + sw*ctAx13) + I*(+ cw*ctZIx13 + sw*ctAIx13),
	CdW   [1,3]  -> cbWx13 + I cbWIx13,
	CuG   [1,3]  -> (ctGx13 + I ctGIx13)*Gstrong,
	
	Cphiq1[3,1]  -> (cpQMx31 + cpQ3x31) + I*(cpQMIx31 + cpQ3Ix31),
	Cphiq3[3,1]  -> cpQ3x31 + I*cpQ3Ix31,
	Cphiu [3,1]  -> cptx31 + I*cptIx31,
	Cphid [3,1]  -> cpbx31 + I*cpbIx31,
	
	Clq1[1,1,3,1]   -> (cQlMx1x31 + cQl3x1x31) + I*(cQlMIx1x31 + cQl3Ix1x31),
	Clq1[2,2,3,1]   -> (cQlMx2x31 + cQl3x2x31) + I*(cQlMIx2x31 + cQl3Ix2x31),
	Clq1[3,3,3,1]   -> (cQlMx3x31 + cQl3x3x31) + I*(cQlMIx3x31 + cQl3Ix3x31),
	Clq3[1,1,3,1]   -> cQl3x1x31 + I*cQl3Ix1x31,
	Clq3[2,2,3,1]   -> cQl3x2x31 + I*cQl3Ix2x31,
	Clq3[3,3,3,1]   -> cQl3x3x31 + I*cQl3Ix3x31,
	Ceq [1,1,3,1]   -> cQex1x31 + I*cQeIx1x31,
	Ceq [2,2,3,1]   -> cQex2x31 + I*cQeIx2x31,
	Ceq [3,3,3,1]   -> cQex3x31 + I*cQeIx3x31,
	Clu [1,1,3,1]   -> ctlx1x31 + I*ctlIx1x31,
	Clu [2,2,3,1]   -> ctlx2x31 + I*ctlIx2x31,
	Clu [3,3,3,1]   -> ctlx3x31 + I*ctlIx3x31,
	Ceu [1,1,3,1]   -> ctex1x31 + I*cteIx1x31,
	Ceu [2,2,3,1]   -> ctex2x31 + I*cteIx2x31,
	Ceu [3,3,3,1]   -> ctex3x31 + I*cteIx3x31,
	
	Clequ1[1,1,3,1] -> ctlSx1x31 + I ctlSIx1x31,
	Clequ1[2,2,3,1] -> ctlSx2x31 + I ctlSIx2x31,
	Clequ1[3,3,3,1] -> ctlSx3x31 + I ctlSIx3x31,
	Clequ3[1,1,3,1] -> ctlTx1x31 + I ctlTIx1x31,
	Clequ3[2,2,3,1] -> ctlTx2x31 + I ctlTIx2x31,
	Clequ3[3,3,3,1] -> ctlTx3x31 + I ctlTIx3x31,
	
	Clequ1[1,1,1,3] -> ctlSx1x13 + I ctlSIx1x13,
	Clequ1[2,2,1,3] -> ctlSx2x13 + I ctlSIx2x13,
	Clequ1[3,3,1,3] -> ctlSx3x13 + I ctlSIx3x13,
	Clequ3[1,1,1,3] -> ctlTx1x13 + I ctlTIx1x13,
	Clequ3[2,2,1,3] -> ctlTx2x13 + I ctlTIx2x13,
	Clequ3[3,3,1,3] -> ctlTx3x13 + I ctlTIx3x13,

	Cledq[1,1,3,1] -> cblSx1x31 + I cblSIx1x31,
	Cledq[2,2,3,1] -> cblSx2x31 + I cblSIx2x31,
	Cledq[3,3,3,1] -> cblSx3x31 + I cblSIx3x31,
	Cledq[1,1,1,3] -> cblSx1x13 + I cblSIx1x13,
	Cledq[2,2,1,3] -> cblSx2x13 + I cblSIx2x13,
	Cledq[3,3,1,3] -> cblSx3x13 + I cblSIx3x13,
	
	(* FCNC charm *)
	Cuphi [3,2]  -> ctpx32 + I ctpIx32,
	Cphiud[3,2]  -> cptbx32 + I cptbIx32,
	CuB   [3,2]  -> (- sw*ctZx32 + cw*ctAx32) + I*(- sw*ctZIx32 + cw*ctAIx32),
	CuW   [3,2]  -> (+ cw*ctZx32 + sw*ctAx32) + I*(+ cw*ctZIx32 + sw*ctAIx32),
	CdW   [3,2]  -> cbWx32 + I cbWIx32,
	CuG   [3,2]  -> (ctGx32 + I ctGIx32)*Gstrong,
	
	Cuphi [2,3]  -> ctpx23 + I ctpIx23,
	Cphiud[2,3]  -> cptbx23 + I cptbIx23,
	CuB   [2,3]  -> (- sw*ctZx23 + cw*ctAx23) + I*(- sw*ctZIx23 + cw*ctAIx23),
	CuW   [2,3]  -> (+ cw*ctZx23 + sw*ctAx23) + I*(+ cw*ctZIx23 + sw*ctAIx23),
	CdW   [2,3]  -> cbWx23 + I cbWIx23,
	CuG   [2,3]  -> (ctGx23 + I ctGIx23)*Gstrong,
	
	Cphiq1[3,2]  -> (cpQMx32 + cpQ3x32) + I*(cpQMIx32 + cpQ3Ix32),
	Cphiq3[3,2]  -> cpQ3x32 + I*cpQ3Ix32,
	Cphiu [3,2]  -> cptx32 + I*cptIx32,
	Cphid [3,2]  -> cpbx32 + I*cpbIx32,
	
	Clq1[1,1,3,2]   -> (cQlMx1x32 + cQl3x1x32) + I*(cQlMIx1x32 + cQl3Ix1x32),
	Clq1[2,2,3,2]   -> (cQlMx2x32 + cQl3x2x32) + I*(cQlMIx2x32 + cQl3Ix2x32),
	Clq1[3,3,3,2]   -> (cQlMx3x32 + cQl3x3x32) + I*(cQlMIx3x32 + cQl3Ix3x32),
	Clq3[1,1,3,2]   -> cQl3x1x32 + I*cQl3Ix1x32,
	Clq3[2,2,3,2]   -> cQl3x2x32 + I*cQl3Ix2x32,
	Clq3[3,3,3,2]   -> cQl3x3x32 + I*cQl3Ix3x32,
	Ceq [1,1,3,2]   -> cQex1x32 + I*cQeIx1x32,
	Ceq [2,2,3,2]   -> cQex2x32 + I*cQeIx2x32,
	Ceq [3,3,3,2]   -> cQex3x32 + I*cQeIx3x32,
	Clu [1,1,3,2]   -> ctlx1x32 + I*ctlIx1x32,
	Clu [2,2,3,2]   -> ctlx2x32 + I*ctlIx2x32,
	Clu [3,3,3,2]   -> ctlx3x32 + I*ctlIx3x32,
	Ceu [1,1,3,2]   -> ctex1x32 + I*cteIx1x32,
	Ceu [2,2,3,2]   -> ctex2x32 + I*cteIx2x32,
	Ceu [3,3,3,2]   -> ctex3x32 + I*cteIx3x32,
	
	Clequ1[1,1,3,2] -> ctlSx1x32 + I ctlSIx1x32,
	Clequ1[2,2,3,2] -> ctlSx2x32 + I ctlSIx2x32,
	Clequ1[3,3,3,2] -> ctlSx3x32 + I ctlSIx3x32,
	Clequ3[1,1,3,2] -> ctlTx1x32 + I ctlTIx1x32,
	Clequ3[2,2,3,2] -> ctlTx2x32 + I ctlTIx2x32,
	Clequ3[3,3,3,2] -> ctlTx3x32 + I ctlTIx3x32,
	
	Clequ1[1,1,2,3] -> ctlSx1x23 + I ctlSIx1x23,
	Clequ1[2,2,2,3] -> ctlSx2x23 + I ctlSIx2x23,
	Clequ1[3,3,2,3] -> ctlSx3x23 + I ctlSIx3x23,
	Clequ3[1,1,2,3] -> ctlTx1x23 + I ctlTIx1x23,
	Clequ3[2,2,2,3] -> ctlTx2x23 + I ctlTIx2x23,
	Clequ3[3,3,2,3] -> ctlTx3x23 + I ctlTIx3x23,

	Cledq[1,1,3,2] -> cblSx1x32 + I cblSIx1x32,
	Cledq[2,2,3,2] -> cblSx2x32 + I cblSIx2x32,
	Cledq[3,3,3,2] -> cblSx3x32 + I cblSIx3x32,
	Cledq[1,1,2,3] -> cblSx1x23 + I cblSIx1x23,
	Cledq[2,2,2,3] -> cblSx2x23 + I cblSIx2x23,
	Cledq[3,3,2,3] -> cblSx3x23 + I cblSIx3x23,
	
	(* FCNC 1l3h up *)
(*
for i in 1 2; do
for lab in qq13 qq11 uu1 qu1 qu8 qu8 qd1 qd8 ud1 ud8 quqd1 quqd8; do
	echo -e "\tC${lab}[3,${i},3,3] -> c${lab}x333${i} + I c${lab}Ix333${i},"
done
echo ""
for lab in qu1 qu8 qu8 qd1 qd8 ud1 ud8 quqd1 quqd8; do
	echo -e "\tC${lab}[3,3,3,${i}] -> c${lab}x3${i}33 + I c${lab}Ix3${i}33,"
done
echo ""
for lab in quqd1 quqd8; do
	echo -e "\tC${lab}[${i},3,3,3] -> c${lab}x33${i}3 + I c${lab}Ix33${i}3,"
done
echo ""
for lab in quqd1 quqd8; do
	echo -e "\tC${lab}[3,3,${i},3] -> c${lab}x${i}333 + I c${lab}Ix${i}333,"
done
echo ""
done
*)
	Cqq13[3,1,3,3] -> cqq13x3331 + I cqq13Ix3331,
	Cqq11[3,1,3,3] -> cqq11x3331 + I cqq11Ix3331,
	Cuu1[3,1,3,3] -> cuu1x3331 + I cuu1Ix3331,
	Cqu1[3,1,3,3] -> cqu1x3331 + I cqu1Ix3331,
	Cqu8[3,1,3,3] -> cqu8x3331 + I cqu8Ix3331,
	Cqu8[3,1,3,3] -> cqu8x3331 + I cqu8Ix3331,
	Cqd1[3,1,3,3] -> cqd1x3331 + I cqd1Ix3331,
	Cqd8[3,1,3,3] -> cqd8x3331 + I cqd8Ix3331,
	Cud1[3,1,3,3] -> cud1x3331 + I cud1Ix3331,
	Cud8[3,1,3,3] -> cud8x3331 + I cud8Ix3331,
	Cquqd1[3,1,3,3] -> cquqd1x3331 + I cquqd1Ix3331,
	Cquqd8[3,1,3,3] -> cquqd8x3331 + I cquqd8Ix3331,

	Cqu1[3,3,3,1] -> cqu1x3133 + I cqu1Ix3133,
	Cqu8[3,3,3,1] -> cqu8x3133 + I cqu8Ix3133,
	Cqu8[3,3,3,1] -> cqu8x3133 + I cqu8Ix3133,
	Cqd1[3,3,3,1] -> cqd1x3133 + I cqd1Ix3133,
	Cqd8[3,3,3,1] -> cqd8x3133 + I cqd8Ix3133,
	Cud1[3,3,3,1] -> cud1x3133 + I cud1Ix3133,
	Cud8[3,3,3,1] -> cud8x3133 + I cud8Ix3133,
	Cquqd1[3,3,3,1] -> cquqd1x3133 + I cquqd1Ix3133,
	Cquqd8[3,3,3,1] -> cquqd8x3133 + I cquqd8Ix3133,

	Cquqd1[1,3,3,3] -> cquqd1x3313 + I cquqd1Ix3313,
	Cquqd8[1,3,3,3] -> cquqd8x3313 + I cquqd8Ix3313,

	Cquqd1[3,3,1,3] -> cquqd1x1333 + I cquqd1Ix1333,
	Cquqd8[3,3,1,3] -> cquqd8x1333 + I cquqd8Ix1333,

	Cqq13[3,2,3,3] -> cqq13x3332 + I cqq13Ix3332,
	Cqq11[3,2,3,3] -> cqq11x3332 + I cqq11Ix3332,
	Cuu1[3,2,3,3] -> cuu1x3332 + I cuu1Ix3332,
	Cqu1[3,2,3,3] -> cqu1x3332 + I cqu1Ix3332,
	Cqu8[3,2,3,3] -> cqu8x3332 + I cqu8Ix3332,
	Cqu8[3,2,3,3] -> cqu8x3332 + I cqu8Ix3332,
	Cqd1[3,2,3,3] -> cqd1x3332 + I cqd1Ix3332,
	Cqd8[3,2,3,3] -> cqd8x3332 + I cqd8Ix3332,
	Cud1[3,2,3,3] -> cud1x3332 + I cud1Ix3332,
	Cud8[3,2,3,3] -> cud8x3332 + I cud8Ix3332,
	Cquqd1[3,2,3,3] -> cquqd1x3332 + I cquqd1Ix3332,
	Cquqd8[3,2,3,3] -> cquqd8x3332 + I cquqd8Ix3332,

	Cqu1[3,3,3,2] -> cqu1x3233 + I cqu1Ix3233,
	Cqu8[3,3,3,2] -> cqu8x3233 + I cqu8Ix3233,
	Cqu8[3,3,3,2] -> cqu8x3233 + I cqu8Ix3233,
	Cqd1[3,3,3,2] -> cqd1x3233 + I cqd1Ix3233,
	Cqd8[3,3,3,2] -> cqd8x3233 + I cqd8Ix3233,
	Cud1[3,3,3,2] -> cud1x3233 + I cud1Ix3233,
	Cud8[3,3,3,2] -> cud8x3233 + I cud8Ix3233,
	Cquqd1[3,3,3,2] -> cquqd1x3233 + I cquqd1Ix3233,
	Cquqd8[3,3,3,2] -> cquqd8x3233 + I cquqd8Ix3233,

	Cquqd1[2,3,3,3] -> cquqd1x3323 + I cquqd1Ix3323,
	Cquqd8[2,3,3,3] -> cquqd8x3323 + I cquqd8Ix3323,

	Cquqd1[3,3,2,3] -> cquqd1x2333 + I cquqd1Ix2333,
	Cquqd8[3,3,2,3] -> cquqd8x2333 + I cquqd8Ix2333,
	
(*
for a in 1 2; do
for i in 1 2; do
for lab in qq11 qq13 qq81 qq83 uu1 uu8 ud1 ud8 qd1 qd8 qu1 qu8; do
	echo -e "\tC${lab}[${i},${i},3,${a}] -> c${lab}x3${a}ii + I c${lab}Ix3${a}ii,"
done
echo ""
for lab in qu1 qu8; do
	echo -e "\tC${lab}[3,${a},${i},${i}] -> c${lab}xii3${a} + I c${lab}Ixii3${a},"
done
echo ""
done
done
*)
	Cqq11[1,1,3,1] -> cqq11x31ii + I cqq11Ix31ii,
	Cqq13[1,1,3,1] -> cqq13x31ii + I cqq13Ix31ii,
	Cqq81[1,1,3,1] -> cqq81x31ii + I cqq81Ix31ii,
	Cqq83[1,1,3,1] -> cqq83x31ii + I cqq83Ix31ii,
	Cuu1[1,1,3,1] -> cuu1x31ii + I cuu1Ix31ii,
	Cuu8[1,1,3,1] -> cuu8x31ii + I cuu8Ix31ii,
	Cud1[1,1,3,1] -> cud1x31ii + I cud1Ix31ii,
	Cud8[1,1,3,1] -> cud8x31ii + I cud8Ix31ii,
	Cqd1[1,1,3,1] -> cqd1x31ii + I cqd1Ix31ii,
	Cqd8[1,1,3,1] -> cqd8x31ii + I cqd8Ix31ii,
	Cqu1[1,1,3,1] -> cqu1x31ii + I cqu1Ix31ii,
	Cqu8[1,1,3,1] -> cqu8x31ii + I cqu8Ix31ii,

	Cqu1[3,1,1,1] -> cqu1xii31 + I cqu1Ixii31,
	Cqu8[3,1,1,1] -> cqu8xii31 + I cqu8Ixii31,

	Cqq11[2,2,3,1] -> cqq11x31ii + I cqq11Ix31ii,
	Cqq13[2,2,3,1] -> cqq13x31ii + I cqq13Ix31ii,
	Cqq81[2,2,3,1] -> cqq81x31ii + I cqq81Ix31ii,
	Cqq83[2,2,3,1] -> cqq83x31ii + I cqq83Ix31ii,
	Cuu1[2,2,3,1] -> cuu1x31ii + I cuu1Ix31ii,
	Cuu8[2,2,3,1] -> cuu8x31ii + I cuu8Ix31ii,
	Cud1[2,2,3,1] -> cud1x31ii + I cud1Ix31ii,
	Cud8[2,2,3,1] -> cud8x31ii + I cud8Ix31ii,
	Cqd1[2,2,3,1] -> cqd1x31ii + I cqd1Ix31ii,
	Cqd8[2,2,3,1] -> cqd8x31ii + I cqd8Ix31ii,
	Cqu1[2,2,3,1] -> cqu1x31ii + I cqu1Ix31ii,
	Cqu8[2,2,3,1] -> cqu8x31ii + I cqu8Ix31ii,

	Cqu1[3,1,2,2] -> cqu1xii31 + I cqu1Ixii31,
	Cqu8[3,1,2,2] -> cqu8xii31 + I cqu8Ixii31,

	Cqq11[1,1,3,2] -> cqq11x32ii + I cqq11Ix32ii,
	Cqq13[1,1,3,2] -> cqq13x32ii + I cqq13Ix32ii,
	Cqq81[1,1,3,2] -> cqq81x32ii + I cqq81Ix32ii,
	Cqq83[1,1,3,2] -> cqq83x32ii + I cqq83Ix32ii,
	Cuu1[1,1,3,2] -> cuu1x32ii + I cuu1Ix32ii,
	Cuu8[1,1,3,2] -> cuu8x32ii + I cuu8Ix32ii,
	Cud1[1,1,3,2] -> cud1x32ii + I cud1Ix32ii,
	Cud8[1,1,3,2] -> cud8x32ii + I cud8Ix32ii,
	Cqd1[1,1,3,2] -> cqd1x32ii + I cqd1Ix32ii,
	Cqd8[1,1,3,2] -> cqd8x32ii + I cqd8Ix32ii,
	Cqu1[1,1,3,2] -> cqu1x32ii + I cqu1Ix32ii,
	Cqu8[1,1,3,2] -> cqu8x32ii + I cqu8Ix32ii,

	Cqu1[3,2,1,1] -> cqu1xii32 + I cqu1Ixii32,
	Cqu8[3,2,1,1] -> cqu8xii32 + I cqu8Ixii32,

	Cqq11[2,2,3,2] -> cqq11x32ii + I cqq11Ix32ii,
	Cqq13[2,2,3,2] -> cqq13x32ii + I cqq13Ix32ii,
	Cqq81[2,2,3,2] -> cqq81x32ii + I cqq81Ix32ii,
	Cqq83[2,2,3,2] -> cqq83x32ii + I cqq83Ix32ii,
	Cuu1[2,2,3,2] -> cuu1x32ii + I cuu1Ix32ii,
	Cuu8[2,2,3,2] -> cuu8x32ii + I cuu8Ix32ii,
	Cud1[2,2,3,2] -> cud1x32ii + I cud1Ix32ii,
	Cud8[2,2,3,2] -> cud8x32ii + I cud8Ix32ii,
	Cqd1[2,2,3,2] -> cqd1x32ii + I cqd1Ix32ii,
	Cqd8[2,2,3,2] -> cqd8x32ii + I cqd8Ix32ii,
	Cqu1[2,2,3,2] -> cqu1x32ii + I cqu1Ix32ii,
	Cqu8[2,2,3,2] -> cqu8x32ii + I cqu8Ix32ii,

	Cqu1[3,2,2,2] -> cqu1xii32 + I cqu1Ixii32,
	Cqu8[3,2,2,2] -> cqu8xii32 + I cqu8Ixii32
	};


Print["--- Save full Lagrangian"]
Ltot >> dim6top_Ltot ;


(*Ltot = << dim6top_Ltot_clean ;*)


Print["--- Compute vertices"]
verttot = FeynmanRules[
        Ltot(*,
        MaxParticles->4*)];


Print["--- Simplify vertices"]
verttot = verttot /.{cw^2->1-sw^2} // Simplify  ;


Print["--- Save vertices"]
verttot >> saved_verttot ;


(*verttot = << verttot_tmp ;*)


Print["--- Write UFO model with those"]
WriteUFO[
        {},Input->verttot,
        AddDecays->False,
        Output->"dim6top_LO_UFO"
        (*Output->"dim6top_LO_UFO_each_coupling_order"*)
];



