(*****************************************************************************************)

(* SM Lagrangian *)

(******************** Gauge F^2 Lagrangian terms*************************)
(*Sign convention from Lagrangian in between Eq. (A.9) and Eq. (A.10) of Peskin & Schroeder.*)
LGauge := Block[{mu,nu,ii,aa},
  ExpandIndices[-1/4 FS[B,mu,nu] FS[B,mu,nu] - 1/4 FS[Wi,mu,nu,ii] FS[Wi,mu,nu,ii] - 1/4 FS[G,mu,nu,aa] FS[G,mu,nu,aa], FlavorExpand->SU2W]];


(********************* Fermion Lagrangian terms*************************)
(*Sign convention from Lagrangian in between Eq. (A.9) and Eq. (A.10) of Peskin & Schroeder.*)
LFermions := Block[{mu},
  ExpandIndices[I*(
    QLbar.Ga[mu].DC[QL, mu] + LLbar.Ga[mu].DC[LL, mu] + uRbar.Ga[mu].DC[uR, mu] + dRbar.Ga[mu].DC[dR, mu] + lRbar.Ga[mu].DC[lR, mu]),
  FlavorExpand->{SU2W,SU2D}]/.{CKM[a_,b_] Conjugate[CKM[a_,c_]]->IndexDelta[b,c], CKM[b_,a_] Conjugate[CKM[c_,a_]]->IndexDelta[b,c]}];




(**************Ghost terms**************************)
(* Now we need the ghost terms which are of the form:             *)
(* - g * antighost * d_BRST G                                     *)
(* where d_BRST G is BRST transform of the gauge fixing function. *)

LGhost := Block[{LGh1,LGhw,LGhs,LGhphi,mu, generators,gh,ghbar,Vectorize,phix1,phix2,togoldstones,doublet,doublet0},
  (* Pure gauge piece *)
  LGh1 = -ghBbar.del[DC[ghB,mu],mu];
  LGhw = -ghWibar.del[DC[ghWi,mu],mu];
  LGhs = -ghGbar.del[DC[ghG,mu],mu];

  (* Scalar pieces: see Peskin pages 739-742 *)
  (* phix1 and phix2 are the real degrees of freedom of GP *)
  (* Vectorize transforms a doublet in a vector in the phi-basis, i.e. the basis of real degrees of freedom *)
  gh    = {ghB, ghWi[1], ghWi[2], ghWi[3]};
  ghbar = {ghBbar, ghWibar[1], ghWibar[2], ghWibar[3]};
  generators = {-I/2 g1 IdentityMatrix[2], -I/2 gw PauliSigma[1], -I/2 gw PauliSigma[2], -I/2 gw PauliSigma[3]};
  doublet = Expand[{(-I phix1 - phix2)/Sqrt[2], Phi1[2]} /. MR$Definitions /. vev -> 0];
  doublet0 = {0, vev/Sqrt[2]};
  Vectorize[{a_, b_}]:= Simplify[{Sqrt[2] Re[Expand[a]], Sqrt[2] Im[Expand[a]], Sqrt[2] Re[Expand[b]], Sqrt[2] Im[Expand[b]]}/.{Im[_]->0, Re[num_]->num}];
  togoldstones := {phix1 -> (GP + GPbar)/Sqrt[2], phix2 -> (-GP + GPbar)/(I Sqrt[2])};
  LGhphi=Plus@@Flatten[Table[-ghbar[[kkk]].gh[[lll]] Vectorize[generators[[kkk]].doublet0].Vectorize[generators[[lll]].(doublet+doublet0)],{kkk,4},{lll,4}]] /.togoldstones;

ExpandIndices[ LGhs + If[FeynmanGauge, LGh1 + LGhw + LGhphi,0], FlavorExpand->SU2W]];


(*Higgs sectors*)


LHiggs := Block[{ii,mu, feynmangaugerules,V1HDpSM,LHiggsKin},
  feynmangaugerules = If[Not[FeynmanGauge], {G0|GP|GPbar ->0}, {}];

  P11 := Module[{jj},Phi1bar[jj] Phi1[jj]];

  V1HDpSM := - mu1 P11 - mu2 SH0^2
           + l1 P11^2 + l2 SH0^4 + l3 P11 SH0^2;

  LHiggsKin := DC[Phi1bar[ii],mu] DC[Phi1[ii],mu] + del[SH0, mu]^2;

  ExpandIndices[LHiggsKin - V1HDpSM, FlavorExpand->{SU2D,SU2W}]/.feynmangaugerules
 ];


LYukawa := Block[{sp,ii,jj,cc,ff1,ff2,ff3,yuk,feynmangaugerules},
  feynmangaugerules = If[Not[FeynmanGauge], {G0|GP|GPbar ->0}, {}];

  yuk = ExpandIndices[
   -yd[ff2, ff3] CKM[ff1, ff2] QLbar[sp, ii, ff1, cc].dR [sp, ff3, cc] Phi1[ii] -
    yl[ff1, ff3] LLbar[sp, ii, ff1].lR [sp, ff3] Phi1[ii] -
    yu[ff1, ff2] QLbar[sp, ii, ff1, cc].uR [sp, ff2, cc] Phi1bar[jj] Eps[ii, jj], FlavorExpand -> SU2D];
  yuk = yuk /. { CKM[a_, b_] Conjugate[CKM[a_, c_]] -> IndexDelta[b, c], CKM[b_, a_] Conjugate[CKM[c_, a_]] -> IndexDelta[b, c]};

  yuk+HC[yuk]/.feynmangaugerules
 ];


L1HDpSM := LGauge + LFermions + LHiggs + LYukawa + LGhost + LCPEven ;