BeginPackage["StandardModel`", "Class`",
"RelativisticKinematics`", "Spinor`"]

(* General *)

NumberOfFamilies::usage = "NumberOfFamilies is the number of complete
fermion families.  It is set to 3."

Vertex::usage = "Vertex[line1, line2, ...] is the Feynman rule for the
vertex.  Each line is of the form {particle} or {particle, p} or
{particle, p, indexes} where particle is the particle type, p is the
particle's momentum, and indexes are the appropriate Lorentz or gauge
group indexes."

KobayashiMaskawa::usage = "KobayashiMaskawa[p1, p2] is the
Kobayashi-Moskawa matrix element connecting particles p1 and p2."

(* SU3 *)

SU3F::usage = "SU3F[a, b, c] stands for the completely anti-symmetric
symbol f(a,b,c).  Commutator[T[a], T[b]] -> I SU3F[a, b, c] T[c] where
T's are the generators of SU3."

SU3D::usage = "SU3D[a, b, c] stands for the completely symmetric
symbol d(a,b,c).  T[a] T[b] + T[b] T[a] -> Delta[a, b]
IdentityMatrix/3 + SU3D[a, b, c] T[c] where T's are the generators of
SU3."

SU3T::usage = "SU3T[a] is the SU3 generator T[a]. SU3T[a, b, c] is the
b-c matrix element of SU3T[a]."

SetAttributes[SU3D, Orderless]

Delta::usage = "Delta[a, b] is 0 unless a == b."

(* Properties *)

ElectricCharge::usage = "ElectricCharge is the property giving the
electric charge of a particle in units of the proton charge."

WeakIsospin3::usage = "WeakIsospin3 is the property giving the third
component of the weak isospin of a particle."

Generation::usage = "Generation is the property giving the generation
number of a fermion."

Spin::usage = "Spin is the property giving the spin of a particle."

Mass::usage = StringJoin[
		If[Head[Mass::usage] === String, Mass::usage, ""],
		"Mass used as a property gives the mass of a particle."]

Propagator::usage = "Propagator[line] is the factor associated with
the particle line, where line can be {type, p} or {type, p, mu, nu}
where type is the type of particle, p is the momentum, and mu and nu
are Lorentz indexes.  Propagator[{type, p__}] is the same as
Propagator[type][p]"

(* Particle Classes *)

Particle::usage = "Particle is the super-class of all particle
classes."

Fermion::usage = "Fermion is the super-class of all the fermions."

Lepton::usage = "Lepton is the super-class of all the leptons."

Quark::usage = "Quark is the super-class of all the quarks."

ChargedLepton::usage = "ChargedLepton is the super-class of all the
charged leptons."

Neutrino::usage = "Neutrino is the super-class of all Neutrinos."

UpTypeQuark::usage = "UpTypeQuark is the super-class of all up-type
quarks."

DownTypeQuark::usage = "DownTypeQuark is the super-class of all
down-type quarks."

Boson::usage = "Boson is the super-class of all bosons."

Vector::usage = "Vector is the super-class of all the vector bosons."

Scalar::usage = "Scalar is the super-class of all the scalar bosons."

MasslessVector::usage = "MasslessVector is the super-class of the
Photon and Gluon."

HeavyVector::usage = "HeavyVector is the super-class of the WBoson and
ZBoson."

FirstGeneration::usage = "FirstGeneration is the super-class of all
the first-generation fermions."

SecondGeneration::usage = "SecondGeneration is the super-class of all
the second-generation fermions."

ThirdGeneration::usage = "ThirdGeneration is the super-class of all
the third-generation fermions."

Electron::usage = "Electron is the lightest charged particle."

Photon::usage = "Photon is the mediator of electromagnetic
interactions."

ElectronNeutrino::usage = "ElectronNeutrino is the neutral isospin partner of
the Electron."

DownQuark::usage = "DownQuark is the lightest d-type quark."

UpQuark::usage = "UpQuark is the lightest u-type quark."

Gluon::usage = "Gluon is the mediator of QCD interactions."

WBoson::usage = "WBoson is the charged mediator of weak interactions."

ZBoson::usage = "ZBoson is the neutral mediator of weak interactions."

Muon::usage = "Muon is the second generation charged lepton."

MuonNeutrino::usage = "MuonNeutrino is the second generation
neutrino."

StrangeQuark::usage = "StrangeQuark is the second generation down-type
quark."

CharmQuark::usage = "CharmQuark is the second generation up-type
quark."

Tau::usage = "Tau is the third-generation charged lepton."

TauNeutrino::usage = "TauNeutrino is the third generation neutrino."

BottomQuark::usage = "BottomQuark is the third generation down-type
quark."

TopQuark::usage = "TopQuark is the (speculated) third generation
up-type quark."

HiggsBoson::usage = "HiggsBoson is the (speculated) minimal standard
Higgs."

(* Units *)

MeV::usage = "MeV is GeV/1000."

GeV::usage = "GeV is 10^9 eV and is used as the basic energy/mass unit."

TeV::usage = "TeV is 1000 GeV."

MilliBarn::usage = "MilliBarn is 10^-31 Meter^2, also used as a cross-section
unit."

GeVToMilliBarn::usage = "GeVToMilliBarn[expr] converts GeV^-2 to
MilliBarns in expr."

(* Coupling Constants *)

Alpha::usage = "Alpha is the electromagnetic coupling constant alpha."

FermiGF::usage = "FermiGF is the fermi G_F constant."

ThetaW::usage = "ThetaW is the weak mixing angle."

StrongAlpha::usage = "StrongAlpha is the QCD coupling constant."

WeakG::usage = "WeakG is the weak interaction coupling constant g."

ElectricE::usage = "ElectricE is the electro-magnetic coupling constant e."

StrongG::usage = "StrongG is the QCD coupling constant g_s."

SetReal[ChargedLepton, Neutrino, UpTypeQuark, DownTypeQuark,
MasslessVector, HeavyVector, Particle, Fermion, Lepton, Quark, Boson,
Vector, Scalar, FirstGeneration, SecondGeneration, ThirdGeneration,
Electron, Photon, ElectronNeutrino, DownQuark, UpQuark, Gluon, WBoson,
ZBoson, Muon, MuonNeutrino, StrangeQuark, CharmQuark, Tau,
TauNeutrino, BottomQuark, TopQuark, HiggsBoson, MeV, GeV, TeV,
MilliBarn, GeVToMilliBarn]

Begin["`Private`"]

(* General *)

NumberOfFamilies = 3

(* Particle Masses *)

N[Mass[Electron]] ^= 0.510999 MeV
N[Mass[Muon]] ^= 105.65839 MeV
N[Mass[Tau]] ^= 1784.1 MeV
N[Mass[UpQuark]] ^= 5.6 MeV
N[Mass[DownQuark]] ^= 9.9 MeV
N[Mass[StrangeQuark]] ^= 199. MeV
N[Mass[CharmQuark]] ^= 1.35 GeV
N[Mass[BottomQuark]] ^= 5. GeV
N[Mass[WBoson]] ^= 80.9 GeV
N[Mass[ZBoson]] ^= 91.0 GeV

(* Kobayashi-Maskawa matrix *)

KMValues = {{0.9755, 0.220, 0.006}, {0.220, 0.9744, 0.046}, {0.011,
0.046, 0.9988}}

KobayashiMaskawa[ChargedLepton, Neutrino] = 1
KobayashiMaskawa[Neutrino, ChargedLepton] = 1

KobayashiMaskawa[a_, b_] := If[Generation[a] === Generation[b], 1, 0] /;
	(DescendentQ[a, ChargedLepton] && DescendentQ[b, Neutrino]) ||
	DescendentQ[b, ChargedLepton] && DescendentQ[a, Neutrino]

N[KobayashiMaskawa[a_?(DescendentQ[UpTypeQuark]), 
	b_?(DescendentQ[DownTypeQuark])]] :=
    	Block[
	{n1 = Generation[a], n2 = Generation[b]}, 
	KMValues[[n1, n2]] /; IntegerQ[n1] && IntegerQ[n2]
    ] /; DescendentQ[b, DownTypeQuark] && DescendentQ[a, UpTypeQuark]

N[KobayashiMaskawa[a_, b_]] := N[KobayashiMaskawa[b, a]] /;
	DescendentQ[a, DownTypeQuark] && DescendentQ[b, UpTypeQuark]

KobayashiMaskawa[a_, b_] := 0 /; 
	(DescendentQ[a, Lepton] && DescendentQ[b, Quark]) ||
	(DescendentQ[b, Lepton] && DescendentQ[a, Quark])

KobayashiMaskawa[a_, b_] := 0 /; ElectricCharge[a] == ElectricCharge[b]

(* Class Structure *)

SetClass[Particle]

SetClass[Fermion, Particle, {Spin->1/2}]

SetClass[Lepton, Fermion]
SetClass[Quark, Fermion]

SetClass[ChargedLepton, Lepton, {ElectricCharge->-1, WeakIsospin3->-1/2}]
SetClass[Neutrino, Lepton, {ElectricCharge->0, Mass->0, WeakIsospin3->1/2}]

SetClass[UpTypeQuark, Quark, {ElectricCharge->2/3, WeakIsospin3->1/2}]
SetClass[DownTypeQuark, Quark, {ElectricCharge->-1/3, WeakIsospin3->-1/2}]

SetClass[Boson, Particle]

SetClass[Vector, Boson, {Spin->1}]

SetClass[MasslessVector, Vector, {Mass->0}]
SetClass[HeavyVector, Vector]

SetClass[Scalar, Boson, {Spin->0}]

SetClass[FirstGeneration, Fermion, {Generation->1}]
SetClass[SecondGeneration, Fermion, {Generation->2}]
SetClass[ThirdGeneration, Fermion, {Generation->3}]

SetClass[Photon, MasslessVector, {ElectricCharge->0}]
SetClass[Gluon, MasslessVector, {ElectricCharge->0}]
SetClass[WBoson, HeavyVector, {ElectricCharge->1}]
SetClass[ZBoson, HeavyVector, {ElectricCharge->0}]
SetClass[HiggsBoson, Scalar, {ElectricCharge->0}]
SetClass[Electron, {FirstGeneration, ChargedLepton}]
SetClass[ElectronNeutrino, {FirstGeneration, Neutrino}]
SetClass[UpQuark, {FirstGeneration, UpTypeQuark}]
SetClass[DownQuark, {FirstGeneration, DownTypeQuark}]
SetClass[Muon, {SecondGeneration, ChargedLepton}]
SetClass[MuonNeutrino, {SecondGeneration, Neutrino}]
SetClass[CharmQuark, {SecondGeneration, UpTypeQuark}]
SetClass[StrangeQuark, {SecondGeneration, DownTypeQuark}]
SetClass[Tau, {ThirdGeneration, ChargedLepton}]
SetClass[TauNeutrino, {ThirdGeneration, Neutrino}]
SetClass[TopQuark, {ThirdGeneration, UpTypeQuark}]
SetClass[BottomQuark, {ThirdGeneration, DownTypeQuark}]

(* Units *)

MeV = GeV/1000
TeV = 1000 GeV
GeVToMilliBarn[expr_] := expr /. GeV-> 1/Sqrt[0.38937966 MilliBarn]

(* Other constants *)

N[Alpha] = 1/137.036
N[FermiGF] = 1.16637 10^-5 / GeV^2
N[ThetaW] = ArcSin[Sqrt[0.228]]
SetReal[Alpha, FermiGF, ThetaW, StrongAlpha]
Unprotect[Sin, Cos]
SetReal[Sin[ThetaW], Cos[ThetaW]]
Protect[Sin, Cos]

ElectricE = Sqrt[4 Pi Alpha]
StrongG = Sqrt[4 Pi StrongAlpha]
WeakG = ElectricE / Sin[ThetaW]

(* Utilitiy Functions *)

FermionQ = DescendentQ[Fermion]

(* Propagators *)

Propagator[{type_?ClassQ, p__}] := Propagator[type][{p}]

Propagator[a_?FermionQ][{p_}] := Propagator[a][p]

Propagator[a_?FermionQ][p_] :=
    I (Slash[p] + Mass[a]) / (DotProduct[p, p] - Mass[a]^2)

Propagator[a_][{p_, mu_, nu_}] := 
    -I (G[mu, nu] - p[mu] p[nu] / Mass[a]^2) / 
			(DotProduct[p, p] - Mass[a]^2) /; 
					DescendentQ[a, HeavyVector]

Propagator[a_][{p_, mu_, nu_}] := -I G[mu, nu] / DotProduct[p, p] /;
					DescendentQ[a, MasslessVector]

Propagator[a_][pl_] := 
     Block[{p = First[Flatten[{pl}]]}, -I  / (DotProduct[p, p] - Mass[a]^2)] /;
					DescendentQ[a, Vector]

Propagator[a_][pl_] :=
     Block[{p = First[Flatten[{pl}]]}, I / (DotProduct[p, p] - Mass[a]^2)] /;
					DescendentQ[a, Scalar]

(* Vertices *)

Vertex[a___, b_, c___] := Vertex[a, {b}, c] /; Head[b] =!= List

Vertex[{a_?FermionQ, ___}, {a_, ___}, {Photon, _, mu_, ___}] := 
	-I ElectricCharge[a] ElectricE DiracGamma[mu]

Vertex[{Photon, _, mu_, ___}, {a_?FermionQ, ___}, {a_, ___}] :=
	-I ElectricCharge[a] ElectricE DiracGamma[mu]

Vertex[{a_?FermionQ, ___}, {Photon, _, mu_, ___}, {a_, ___}] :=
	-I ElectricCharge[a] ElectricE DiracGamma[mu]

Vertex[{a_?FermionQ, ___}, {b_?FermionQ, ___}, {WBoson, _, mu_}] :=
    -I (WeakG / Sqrt[2]) DiracGamma[mu] ** HelicityProjection[Left] *
	KobayashiMaskawa[a, b]

Vertex[{a_?FermionQ, ___}, {WBoson, _, mu_}, {b_?FermionQ, ___}] :=
    -I (WeakG / Sqrt[2]) DiracGamma[mu] ** HelicityProjection[Left] *
	KobayashiMaskawa[a, b]

Vertex[{WBoson, _, mu_}, {a_?FermionQ, ___}, {b_?FermionQ, ___}] :=
    -I (WeakG / Sqrt[2]) DiracGamma[mu] ** HelicityProjection[Left] *
	KobayashiMaskawa[a, b]

Vertex[{a_?FermionQ, ___}, {a_, ___}, {ZBoson, _, mu_}] :=
    -I (WeakG / Cos[ThetaW]) DiracGamma[mu] **
 (WeakIsospin3[a] HelicityProjection[Left] - Sin[ThetaW]^2 ElectricCharge[a])

Vertex[{a_?FermionQ, ___}, {ZBoson, _, mu_}, {a_, ___}] :=
    -I (WeakG / Cos[ThetaW]) DiracGamma[mu] **
 (WeakIsospin3[a] HelicityProjection[Left] - Sin[ThetaW]^2 ElectricCharge[a])

Vertex[{ZBoson, _, mu_}, {a_?FermionQ, ___}, {a_, ___}] :=
    -I (WeakG / Cos[ThetaW]) DiracGamma[mu] **
 (WeakIsospin3[a] HelicityProjection[Left] - Sin[ThetaW]^2 ElectricCharge[a])

threeVectorVertex[{p1_, mu1_}, {p2_, mu2_}, {p3_, mu3_}] :=
   (p1-p2)[mu3]G[mu1, mu2] + (p2-p3)[mu1]G[mu2, mu3] + (p3-p1)[mu2]G[mu3, mu1]

(* all momenta incoming and order of particles is w+, w- , photon *)
Vertex[{WBoson, p1_, mu1_, ___}, {WBoson, p2_, mu2_, ___}, 
		{Photon, p3_, mu3_, ___}] := 
     -I ElectricE threeVectorVertex[{p1, mu1}, {p2, mu2}, {p3, mu3}] 

Vertex[l1:{Photon, ___}, l2:{WBoson, ___}, l3:{WBoson, ___}] :=
		   Vertex[l2, l3, l1]

Vertex[l1:{WBoson, ___}, l2:{Photon, ___}, l3:{WBoson, ___}] :=
		   Vertex[l3, l1, l2]

(* all momenta incoming and order of particles is w+, w- , z *)
Vertex[{WBoson, p1_, mu1_, ___}, {WBoson, p2_, mu2_, ___}, 
		{ZBoson, p3_, mu3_, ___}] := 
     -I WeakG Cos[ThetaW] threeVectorVertex[{p1, mu1}, {p2, mu2}, {p3, mu3}] 

Vertex[l1:{ZBoson, ___}, l2:{WBoson, ___}, l3:{WBoson, ___}] :=
		   Vertex[l2, l3, l1]

Vertex[l1:{WBoson, ___}, l2:{ZBoson, ___}, l3:{WBoson, ___}] :=
		   Vertex[l3, l1, l2]

fourVectorVertex[mu1_, mu2_, mu3_, mu4_] :=
       2 G[mu1, mu2] G[mu3, mu4] - G[mu1, mu3] G[mu2, mu4] - 
       G[mu1, mu4] G[mu2, mu3]

Vertex[{WBoson, _, mu1_, ___}, {WBoson, _, mu2_, ___}, 
		{WBoson, _, mu3_, ___}, {WBoson, _, mu4_, ___}] := 
	I WeakG^2 fourVectorVertex[mu1, mu2, mu3, mu4]       

Vertex[{WBoson, _, mu1_, ___}, {WBoson, _, mu2_, ___}, 
		{ZBoson, _, mu3_, ___}, {ZBoson, _, mu4_, ___}] := 
	-I WeakG^2 Cos[ThetaW]^2 fourVectorVertex[mu1, mu2, mu3, mu4]       

Vertex[{WBoson, _, mu1_, ___}, {WBoson, _, mu2_, ___}, 
		{ZBoson, _, mu3_, ___}, {Photon, _, mu4_, ___}] := 
	-I ElectricE WeakG Cos[ThetaW] fourVectorVertex[mu1, mu2, mu3, mu4]

Vertex[{WBoson, _, mu1_, ___}, {WBoson, _, mu2_, ___}, 
		{Photon, _, mu3_, ___}, {Photon, _, mu4_, ___}] := 
	-I ElectricE^2 fourVectorVertex[mu1, mu2, mu3, mu4]       

Vertex[{a_?FermionQ, ___}, {a_, ___}, {HiggsBoson, ___}] :=
		     -I WeakG Mass[a] / (2 Mass[WBoson])

Vertex[{WBoson, _, mu1_, ___}, {WBoson, _, mu2_, ___}, {HiggsBoson, ___}] :=
		-I WeakG Mass[WBoson] G[mu1, mu2]

Vertex[{ZBoson, _, mu1_, ___}, {ZBoson, _, mu2_, ___}, {HiggsBoson, ___}] :=
		-I WeakG Mass[ZBoson] G[mu1, mu2]/Cos[ThetaW]

Vertex[{q_?(DescendentQ[Quark]), _, a_, ___}, {q_, _, b_, ___},
       {Gluon, _, mu_, c_, ___}] :=
    	-I StrongG SU3T[c, a, b] DiracGamma[mu]

Vertex[{q_?(DescendentQ[Quark]), _, a_, ___}, {Gluon, _, mu_, c_, ___},
	{q_, _, b_, ___}] := 
	-I StrongG SU3T[c, a, b] DiracGamma[mu]

Vertex[{Gluon, _, mu_, c_, ___}, {q_?(DescendentQ[Quark]), _, a_, ___}, 
	{q_, _, b_, ___}] := 
	-I StrongG SU3T[c, a, b] DiracGamma[mu]

Vertex[{Gluon, p1_, mu1_, a1_, ___}, {Gluon, p2_, mu2_, a2_, ___},
       {Gluon, p3_, mu3_, a3_, ___}] :=
 -I StrongG SU3F[a1, a2, a3] threeVectorVertex[{p1, mu1}, {p2, mu2}, {p3, mu3}]

oneElement[{a_, b_, c_, d_, e_}, {m1_, m2_, m3_, m4_}] :=
   SU3F[a, b, e] SU3F[c, d, e] (G[m1, m3]G[m2, m4] - G[m1, m4]G[m2, m3])

Vertex[{Gluon, p1_, mu1_, a1_, ___}, {Gluon, p2_, mu2_, a2_, ___},
       {Gluon, p3_, mu3_, a3_, ___}, {Gluon, p4_, mu4_, a4_, ___}] :=
    -I StrongG^2 Block[{a5 = Unique["a"]},
	       oneElement[{a1, a2, a3, a4, a5}, {mu1, mu2, mu3, mu4}] +
	       oneElement[{a1, a3, a4, a2, a5}, {mu1, mu3, mu4, mu2}] +
	       oneElement[{a1, a4, a2, a3, a5}, {mu1, mu4, mu2, mu3}]]

SpaceTimeDimension = 4;
DiracGammaSize = 4;

(* SU3 stuff *)

Literal[SU3F[a___, b_, c_, d___]] := -SU3F[a, c, b, d] /; !OrderedQ[{b, c}]
Literal[SU3F[___, a_, ___, a_, ___]] := 0
SU3F[1,2,3] = 1
SU3F[1,4,7] = SU3F[2,4,6] = SU3F[2,5,7] = SU3F[3,4,5] = 1/2
SU3F[1,5,6] = SU3F[3,6,7] = -1/2
SU3F[4,5,8] = SU3F[6,7,8] = Sqrt[3] / 2

SU3D[1,1,8] = SU3D[2,2,8] = SU3D[3,3,8] = 1/Sqrt[3]
SU3D[8,8,8] = -1/Sqrt[3]
SU3D[1,4,6] = SU3D[1,5,7] = SU3D[2,5,6] = SU3D[3,4,4] = SU3D[3,5,5] = 1/2
SU3D[2,4,7] = SU3D[3,6,6] = SU3D[3,7,7] = -1/2
SU3D[4,4,8] = SU3D[5,5,8] = SU3D[6,6,8] = SU3D[7,7,8] = -1/2/Sqrt[3]
SU3D[_Integer, _Integer, _Integer] = 0

End[]

EndPackage[]

