(* Changes after 3/17/91:

5/29/91:

Conditionally change DownValue to DownValues in version 2.0. The
packages are now 2.0 compatible.

Remove MixedSlash.

3/19/91: Change 

HelicityProjection/:
  Literal[NonCommutativeMultiply[a___, s_[p_, pol___], 
  HelicityProjection[pol2_], b___]] := If[
   SameQ[pol, -pol2], NonCommutativeMultiply[a, s[p, -pol2], b], 0] /;
	      MemberQ[{SpinorUbar, SpinorVbar}, s]

to

HelicityProjection/:
  Literal[NonCommutativeMultiply[a___, s_[p_, pol_],
   HelicityProjection[pol2_], b___]] := 
   (1 + pol pol2) / 2 NonCommutativeMultiply[a, s[p, -pol2], b] /;
    MemberQ[{SpinorUbar, SpinorVbar}, s]

HelicityProjection/:
  Literal[NonCommutativeMultiply[a___, s_[p_],
   HelicityProjection[pol2_], b___]] := 
   NonCommutativeMultiply[a, s[p, -pol2], b] /;
    MemberQ[{SpinorUbar, SpinorVbar}, s]

*)

(* ------------------ Spinor.m ------------------ *)
BeginPackage["Spinor`", "NonCommutativeMultiply`", "RelativisticKinematics`"];


Clear[Contract, DiracGamma, G, Slash, GammaTrace, AbsSquared, DiracGamma5,
SpinorU, SpinorUbar, HelicityProjection, PrepareIndex,
PhysicsForm, PhysicsFormPrint, ConvertToGamma5,
ConvertToChiralFermions, ExpandSlash, TraceAllTheWay]

MasslessVectorPolarization::usage = "MasslessVectorPolarization[
p, mu, pol] represents the polarization vector of a massless vector particle
with momentum p, polarization pol and index mu."

HeavyVectorPolarization::usage = "HeavyVectorPolarization[p, mu,
pol] represents the polarization vector of a massive vector particle
with momentum p, polarization pol and index mu."

G::usage = "G[mu, nu] is the ordinary spacetime metric. mu and nu range
from 1 to SpaceTimeDimension."

Contract::usage = "Contract[expr, index] contracts index in expr.
All subexpressions of form G[index, _] or p_[index, _] can potentially
be effected. Contract[expr, {index1, index2, ...}] does the contraction
sequentially."

DiracGamma::usage = "DiracGamma[mu] is the Dirac Gamma matrix.
It is set to be non-commutative. Note that DotProduct[DiracGamma, p]
is simplified to Slash[p]."

Slash::usage = "Slash[p] is the same as Contract[p[mu]DiracGamma[mu], mu].
Slash is set to be non-commutative."

GammaTrace::usage = "GammaTrace[expr] is the trace (in spinor space) of expression.
Note that no division by DiracGammaSize is carried out, ie
GammaTrace[1] gives DiracGammaSize."

Options[GammaTrace] = {DoEverything->True}

TraceAllTheWay = False

DiracGammaSize::usage = "DiracGammaSize is the dimension of the Dirac
spinors. Its default value is 4."

AbsSquared::usage = "AbsSquared[expr] is the absolute value squared of
expr."

SetAttributes[AbsSquared, Listable]

DiracGamma5::usage = "DiracGamma5 is mainly meaningfull for
SpaceTimeDimension of 4. Some of its features will work in any
even spacetime dimension."

SpinorU::usage = "SpinorU[p] is a spinor object of a four-vector p.
SpinorU[p, pol] is a spinor object of particular polarization."

SpinorUbar::usage = "SpinorUbar[p] is a spinor object of a four-vector p.
SpinorUbar[p, pol] is a spinor object of particular polarization."

SpinorV::usage = "SpinorV[p] is a spinor object of an anti-particle
of four-vector p. SpinorV[p, pol] is a spinor object of particular
polarization."

SpinorVbar::usage = "SpinorVbar[p] is a spinor object of an anti-particle
of four-vector p. SpinorVbar[p, pol] is a spinor object of particular
polarization."

HelicityProjection::usage = "HelicityProjection[pol] is the helicity
projection operator (1 +/- Gamma5)/2."

PrepareIndex::usage = "PrepareIndex[a, b, c, ...] makes future
contractions involving a, b, c, ... work faster."

ConvertToGamma5::usage = "ConvertToGamma5[expr] attempts to convert all
occurrences of HelicityProjection and spinors of specific chirality
into products involving DiracGamma5."

SetAttributes[ConvertToGamma5, Listable]

ConvertToChiralFermions::usage = "ConvertToChiralFermions[expr] attempts
to convert all occurrences of DiracGamma5 and HelicityProjection in expr
into products of chiral spinors."

SetAttributes[ConvertToChiralFermions, Listable]

PhysicsForm::usage = "PhysicsForm[expr_] outputs expr in a form
resembling ordinary Physics notation."

PhysicsFormPrint::usage = "PhysicsFormPrint[expr_] prints expr in a form
resembling ordinary Physics notation."

ExpandSlash::usage = "ExpandSlash[expr] expands slash of sums to sums
of slashes in expr."

SetAttributes[ExpandSlash, Listable]

DebugMetric :=
  On[AbsSquared, NonCommutativeMultiply, GammaTrace, SpinorU, SpinorUbar,
  SpinorV, SpinorVbar]

Begin["`Private`"];


Clear[PhysicsForm1, AbsSquared1, Trace1, ContractRules, ConvertToST1]

(*  Some small functions  *)

ExpandSlash[expr_] := NonCommutativeExpand[
                   expr //. {Slash[a_+b_] :> Slash[a]+Slash[b],
				Slash[c_?NumberQ p_] :> c Slash[p]}]

PrepareIndex[a__] :=
  (Function[
    x,
    (ContractRules[x] = Dispatch[ContractRules[x]];
    	DotProduct[p_, Conjugate[x]] := p[Conjugate[x]];
     	DotProduct[Conjugate[x], p_] := p[Conjugate[x]];
     	prot = Unprotect[Conjugate];
     	Conjugate[x][y_] := G[Conjugate[x], y];
     	Conjugate[x][mu_] := G[Conjugate[x], mu];
        PhysicsForm1[p_[x]] := Subscripted[p[x]];
     	Protect[Release[prot]];
     Block[
     	{h = x},
     	While[Head[h]=!=Symbol, h = Head[h]];
     	Release[h]/: DotProduct[x, p_] := p[x];
     	Release[h]/: x[y_] := G[x, y];
     	Release[h]/: DotProduct[p_, x] := p[x];
     ];
    )
  ] /@ {a};)

Combinations[l_List, k_Integer] := Flatten[
             Block[{i, ll=Length[l], m = Table[Unique["a"], {k}]},
             Table[l[[Release[m]]], Release[{m[[1]], ll-k+1}],
             Release[Sequence @@ Table[{m[[i]], m[[i-1]]+1, ll-k+i}, {i,
             2, k}]]]], k-1]

SignedCombinations[l_List, k_Integer] := Flatten[
             Block[{c = (-1)^(k (k+1)/2), i, ll=Length[l],
             m = Table[Unique["a"], {k}]},
             Table[{c (-1)^(Plus @@ m), l[[Release[m]]]},
             Release[{m[[1]], ll-k+1}],
             Release[Sequence @@ Table[{m[[i]], m[[i-1]]+1, ll-k+i}, {i,
             2, k}]]]], k-1]

ReverseInPlace[p___] := Sequence @@ Reverse[{p}]

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

(*  Set NonCommutative rules  *)

SetNonCommutative[HelicityProjection, DiracGamma, DiracGamma5, SpinorU,
SpinorUbar, SpinorV, SpinorVbar, Slash, maybeHP]


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

(*  Set Real rules  *)

SetReal[G]

Unprotect[Plus, Times]
SetReal[Plus, Times]
Protect[Plus, Times]



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

(* Conjugation rules *)

Unprotect[Pi]
Pi/: Conjugate[Pi] = Pi
Protect[Pi]

NonCommutativeMultiply/: 
	Conjugate[Literal[NonCommutativeMultiply[a__]]] := 
		NonCommutativeMultiply @@ (Conjugate /@ Reverse[{a}])

DiracGamma/: Conjugate[DiracGamma[nu_]] := DiracGamma[Conjugate[nu]]

Slash/: Conjugate[Slash[p_]] := Slash[Conjugate[p]]

DiracGamma5/: Conjugate[DiracGamma5] = -DiracGamma5

(*
SetReal[DiracGamma5]
*)

Eps/: Conjugate[Eps[a__]] := Eps[Sequence @@ Conjugate /@ {a}]

SpinorU/: Conjugate[SpinorU[p__]] := SpinorUbar[p]

SpinorUbar/: Conjugate[SpinorUbar[p__]] := SpinorU[p]

SpinorV/: Conjugate[SpinorV[p__]] := SpinorVbar[p]

SpinorVbar/: Conjugate[SpinorVbar[p__]] := SpinorV[p]

Unprotect[Conjugate]
Conjugate[Power[a_, b_?NumberQ]] := Power[Conjugate[a], b]
Protect[Conjugate]


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

(*  Commutator rules  *)

Slash/: Commutator[Slash[p1_], Slash[p2_]] :=
        2DotProduct[p1,p2] - Slash[p2]**Slash[p1]

Slash/: Commutator[Slash[p_], DiracGamma[mu_]] :=
        2 p[mu] - DiracGamma[mu]**Slash[p]
Slash/: Commutator[DiracGamma[mu_], Slash[p_]] :=
        2 p[mu] - Slash[p]**DiracGamma[mu]

Slash[-p_] := -Slash[p]

DiracGamma/: Commutator[DiracGamma[mu_], DiracGamma[nu_]] :=
             2 G[mu, nu] - DiracGamma[nu]**DiracGamma[mu]


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

(*  Dot Product Rules  *)

Literal[DotProduct[Slash[p_], q_]] := DotProduct[p, q]

DotProduct[DiracGamma[mu_], p_] := p[mu]

DotProduct[DiracGamma[mu_], DiracGamma[nu_]] := G[mu, nu]


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

(*  Gammalike rules  *)

GammaLike[_DiracGamma] = True
GammaLike[_Slash] = True
GammaLike[_]  = False
GammaLike[a_, b__] := And @@ (GammaLike /@ {a, b})


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

(*  SpinorQ and SpinorbarQ rules  *)

SpinorQ[_SpinorU] = True
SpinorQ[_SpinorV] = True
SpinorQ[_] = False
SpinorbarQ[_SpinorUbar] = True
SpinorbarQ[_SpinorVbar] = True
SpinorbarQ[_] = False


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

(*  CommutativeQ rules  *)

CommutativeQ[_NonCommutativeMultiply] := False


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

(*  Left & Right rules  (as opposed to wrong and right rules!!)  *)
Unprotect[Left, Right]

Left/: -Left = Right
Left/: Left==Right = False
Left/: Left^2 = 1
Left/: Left Right = -1
Right/: Right^2 = 1
Right/: Right==Left = False
Right/: -Right = Left


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

(*  Slash rules  *)

Slash/: Literal[NonCommutativeMultiply[a___, Slash[p_], Slash[p_], b___]] :=
  DotProduct[p, p] NonCommutativeMultiply[a, b]
  
Slash/:	Literal[NonCommutativeMultiply[
	a___, SpinorUbar[p_, pol___], Slash[p_],b___]
 	] := Mass[p] NonCommutativeMultiply[a, SpinorUbar[p, -pol], b]

Slash/:	Literal[NonCommutativeMultiply[
	a___, SpinorVbar[p_, pol___], Slash[p_],b___]
 	] := -Mass[p] NonCommutativeMultiply[a, SpinorVbar[p, -pol], b]

Slash/:	Literal[NonCommutativeMultiply[
	a___, Slash[p_], SpinorU[p_, pol___], b___]
	] := Mass[p] NonCommutativeMultiply[a, SpinorU[p, -pol], b]
		
Slash/:	Literal[NonCommutativeMultiply[
	a___, Slash[p_], SpinorV[p_, pol___], b___]
	] := -Mass[p] NonCommutativeMultiply[a, SpinorV[p, -pol], b]


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

(* Collecting NonCommutativeMultiply terms *)

GammaTrace/:
  GammaTrace[Literal[NonCommutativeMultiply[a__, SpinorU[p_, pol___]]], opts1___] *
  GammaTrace[Literal[NonCommutativeMultiply[SpinorUbar[p_, pol___], b__]], opts2___] :=
  	GammaTrace[
  		NonCommutativeMultiply[
  			a, 
  			Slash[p] + Mass[p], 
  			HelicityProjection[pol], 
  			b
  		], opts1, opts2
  	]

GammaTrace/:
  GammaTrace[Literal[NonCommutativeMultiply[a__, SpinorV[p_, pol___]]], opts1___] *
  GammaTrace[Literal[NonCommutativeMultiply[SpinorVbar[p_, pol___], b__]], opts2___] :=
  	GammaTrace[
  		NonCommutativeMultiply[
  			a, 
  			Slash[p] - Mass[p], 
  			HelicityProjection[pol], 
  			b
  		], opts1, opts2
  	]


NonCommutativeMultiply/:
	Literal[NonCommutativeMultiply[SpinorUbar[_,pol1_],
		b___, SpinorU[_,pol2_]]] := 0 /;
		GammaLike[b] &&
		((pol1 == pol2 && EvenQ[Length[{b}]]) ||
		(pol1 == -pol2 && OddQ[Length[{b}]]))

NonCommutativeMultiply/:
	Literal[NonCommutativeMultiply[SpinorVbar[_,pol1_],
		b___, SpinorV[_,pol2_]]] := 0 /;
		GammaLike[b] &&
		((pol1 == pol2 && EvenQ[Length[{b}]]) ||
		(pol1 == -pol2 && OddQ[Length[{b}]]))

NonCommutativeMultiply/:
	Literal[NonCommutativeMultiply[SpinorVbar[_, pol1_],
		b___, SpinorU[_, pol2_]]] := 0 /;
		GammaLike[b] &&
		((pol1 == pol2 && EvenQ[Length[{b}]]) ||
		(pol1 == -pol2 && OddQ[Length[{b}]]))

NonCommutativeMultiply/:
	Literal[NonCommutativeMultiply[SpinorUbar[_, pol1_],
		b___, SpinorV[_, pol2_]]] := 0 /;
		GammaLike[b] &&
		((pol1 == pol2 && EvenQ[Length[{b}]]) ||
		(pol1 == -pol2 && OddQ[Length[{b}]]))


(*
NonCommutativeMultiply/:
  Literal[
  	NonCommutativeMultiply[a__, SpinorU[p_, pol___]]] *
    Literal[NonCommutativeMultiply[SpinorUbar[p_, pol___], b__]
   ] :=
      NonCommutativeMultiply[
      	a, 
      	Slash[p] + Mass[p], 
      	HelicityProjection[pol], 
      	b
      ]

NonCommutativeMultiply/:
  Literal[NonCommutativeMultiply[a__, SpinorV[p_, pol___]]] *
  Literal[NonCommutativeMultiply[SpinorVbar[p_, pol___], b__]] :=
  	NonCommutativeMultiply[
  		a, 
  		Slash[p] - Mass[p], 
  		HelicityProjection[pol], 
  		b
  	]
*)


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


(*  other junk  *)

G[a_, b_] := G[b, a] /; !OrderedQ[{a, b}]


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

(* ConvertToGamma5 *)

ConvertToGamma5[expr_] := expr //. ConvertToGamma5Rules

ConvertToGamma5Rules = {
    HelicityProjection[Left] ->	(1 - DiracGamma5) / 2,
    HelicityProjection[Right] -> (1 + DiracGamma5) / 2,
    HelicityProjection[l_] -> (1 + l DiracGamma5) / 2,
    HelicityProjection[] -> 1,
    SpinorUbar[p_, pol_] :> SpinorUbar[p]**HelicityProjection[-pol],
    SpinorVbar[p_, pol_] :> SpinorVbar[p]**HelicityProjection[-pol],
    SpinorU[p_, pol_] :> HelicityProjection[pol]**SpinorU[p],
    SpinorV[p_, pol_] :> HelicityProjection[pol]**SpinorV[p]
    }


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

(* DiracGamma5 and HelicityProjection products: general rules *)

DiracGamma5/: 
	Literal[
		NonCommutativeMultiply[
			a___, 
			DiracGamma5, 
			DiracGamma5, 
			b___
		]
	] := NonCommutativeMultiply[a, b]

DiracGamma5/: 
	Literal[
		NonCommutativeMultiply[
			a___, 
			b_, 
			DiracGamma5, 
			c___
		]
	] :=	-NonCommutativeMultiply[a, DiracGamma5, b, c] /; GammaLike[b]

DiracGamma5/: 
	Literal[
		NonCommutativeMultiply[
			a___, 
			HelicityProjection[Left], 
			DiracGamma5, 
			b___
		]
	] :=	- NonCommutativeMultiply[a, HelicityProjection[Left], b]

DiracGamma5/:
  Literal[
  	NonCommutativeMultiply[
  		a___, 
  		HelicityProjection[Right], 
  		DiracGamma5, 
  		b___
  	]
  ] :=	NonCommutativeMultiply[a, HelicityProjection[Right], b]



HelicityProjection/: 
    NonCommutativeMultiply[a___, HelicityProjection[], b___] :=
    NonCommutativeMultiply[a, b]

HelicityProjection/:
  Literal[
  	NonCommutativeMultiply[
  		a___, 
  		b_, 
  		HelicityProjection[pol_], 
  		c___
  	]
  ] :=	NonCommutativeMultiply[a, HelicityProjection[-pol], b, c] /; 
  				GammaLike[b]

HelicityProjection/:
  Literal[NonCommutativeMultiply[a___, s_[p_, pol_],
   HelicityProjection[pol2_], b___]] := 
   (1 - pol pol2) / 2 NonCommutativeMultiply[a, s[p, -pol2], b] /;
    MemberQ[{SpinorUbar, SpinorVbar}, s]

HelicityProjection/:
  Literal[NonCommutativeMultiply[a___, s_[p_],
   HelicityProjection[pol2_], b___]] := 
   NonCommutativeMultiply[a, s[p, -pol2], b] /;
    MemberQ[{SpinorUbar, SpinorVbar}, s]

HelicityProjection/:
  Literal[
  	NonCommutativeMultiply[
  		a___, 
  		HelicityProjection[l1_], 
  		HelicityProjection[l2_], 
  		b___
  	]
  ] :=	(1 + l1 l2) / 2 NonCommutativeMultiply[a, HelicityProjection[l1], b]

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

(*  Vector Polarization rules  *)

MasslessVectorPolarization/: AbsSquared1[
	MasslessVectorPolarization[p_, mu_]] := - G[mu, Conjugate[mu]]

HeavyVectorPolarization/: AbsSquared1[
	HeavyVectorPolarization[p_, mu_]] := 
	- G[mu, Conjugate[mu]] + p[mu] p[Conjugate[mu]] / Mass[p]^2

MasslessVectorPolarization/:
	MasslessVectorPolarization[p_, mu_] *
	Conjugate[MasslessVectorPolarization[p_, nu_]] :=
	- G[mu, Conjugate[nu]]

HeavyVectorPolarization/:
	HeavyVectorPolarization[p_, mu_] *
	Conjugate[HeavyVectorPolarization[p_, nu_]] :=
	- G[mu, Conjugate[nu]] + p[mu] p[Conjugate[nu]] / Mass[p]^2

MasslessVectorPolarization/:
	MasslessVectorPolarization[p_, mu_, ___] * p_[mu_] = 0

HeavyVectorPolarization/:
	HeavyVectorPolarization[p_, mu_, ___] * p_[mu_] = 0



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

(*  SpinorU, Ubar, V and Vbar rules  *)

(*  temporarily leave these out ....
SpinorU/: Literal[NonCommutativeMultiply[
	a___, SpinorU[p_], SpinorUbar[p_], b___]] :=
	NonCommutativeMultiply[a, Slash[p] + Mass[p], b]

SpinorU/:
  Literal[
  	NonCommutativeMultiply[
  		a___, 
  		SpinorU[p_, pol1___], 
  		SpinorUbar[p_, pol2___], 
  		b___
  	]
  ] :=
    NonCommutativeMultiply[
    	a, 
	HelicityProjection[pol1],
    	Slash[p] + Mass[p], 
    	HelicityProjection[-pol2], 
    	b
    ]

SpinorV/: 
	Literal[
		NonCommutativeMultiply[a___, SpinorV[p_], SpinorVbar[p_], b___]
	] :=	NonCommutativeMultiply[a, Slash[p] - Mass[p], b]

SpinorV/:
  Literal[
  	NonCommutativeMultiply[
  		a___, 
  		SpinorV[p_, pol1___], 
  		SpinorVbar[p_, pol2___], 
  		b___
  	]
  ] :=
    NonCommutativeMultiply[
    	a,
	HelicityProjection[pol1], 
    	Slash[p] - Mass[p], 
    	HelicityProjection[-pol2], 
    	b
    ]
*)
(*  try leaving these out in place of the identical rules in AbsSq  
Literal[NonCommutativeMultiply[SpinorUbar[p_], a___, SpinorU[p_]]] :=
    GammaTrace[NonCommutativeMultiply[Slash[p] + Mass[p], a], DoEverything->False]

Literal[NonCommutativeMultiply[SpinorVbar[p_], a___, SpinorV[p_]]] :=
    GammaTrace[NonCommutativeMultiply[Slash[p] - Mass[p], a],DoEverything->False]

Literal[NonCommutativeMultiply[
		SpinorUbar[p_, pol1___],
		a___,
		SpinorU[p_, pol2___]]] :=
    GammaTrace[NonCommutativeMultiply[
		HelicityProjection[pol2],
		Slash[p] + Mass[p], 
		HelicityProjection[-pol1], 
		a],DoEverything->False]

Literal[NonCommutativeMultiply[
		SpinorVbar[p_, pol1___], 
		a___, 
		SpinorV[p_, pol2___]]] :=
    GammaTrace[NonCommutativeMultiply[
		HelicityProjection[pol2],
    		Slash[p] - Mass[p], 
    		HelicityProjection[-pol1], 
    		a],DoEverything->False]
*)

Literal[NonCommutativeMultiply[
  		a_?SpinorbarQ, 
  		b___,
  		c_?SpinorQ, 
  		e_?SpinorbarQ,
  		f___, 
  		g_?SpinorQ]] :=
    NonCommutativeMultiply[a, b, c] * NonCommutativeMultiply[e, f, g]

Literal[NonCommutativeMultiply[
  		a_?SpinorbarQ, 
  		b___, 
  		c_?SpinorbarQ, 
  		d___, 
  		e_?SpinorQ, 
  		f___, 
  		g_?SpinorQ]] :=
    NonCommutativeMultiply[a, b, f, g] * NonCommutativeMultiply[c, d, e]
    	

GammaTrace[Literal[NonCommutativeMultiply[
			a__, 
			g_?SpinorbarQ, 
			b___, 
			h_?SpinorQ, 
			c___]], opts___] :=
  GammaTrace[NonCommutativeMultiply[a, c] * NonCommutativeMultiply[g, b, h], opts]


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

(*  Contract  *)

Contract[a_List, x_] := Contract[#, x]& /@ a

Contract[expr_, x_] := Contract1[Expand[expr], x]

Contract1[a_Plus, x_] := Contract1[#, x]& /@ a

Contract1[expr_, {}] := expr

Contract1[expr_, l_List] :=
  Contract1[Contract1[expr, First[l]], Rest[l]]

Contract1[a_ b_, mu_] := a Contract1[b, mu] /; FreeQ[a, mu]

Contract1[expr_, mu_] := expr //. ContractRules[mu]

ContractRules[mu_] :=
  {
    Literal[G[mu, nu_]]^2 -> G[nu, nu],

    Literal[G[nu_, mu]]^2 -> G[nu, nu],

    G[mu, mu] -> SpaceTimeDimension,

    Literal[G[mu, nu_]]p_[a___, mu, b___] -> p[a, nu, b],

    Literal[G[nu_, mu]]p_[a___, mu, b___] -> p[a, nu, b],

    Literal[Eps[a___, mu, b___]] p_[mu] :> Eps[a, p, b],

    Literal[Eps[a1___, mu, b1___] Eps[a2___, mu, b2___]] :>
	- Block[
  	  {c = {mu}, m, l, ac = {a1, b1}, bc = {a2, b2},
		a = {a1, mu, b1}, b = {a2, mu, b2}},
  	  m = Length[a];
  	  l = Join[Table[1, {m}], Table[2, {m}]];
  	  Signature[Transpose[{Join[a, b], l}]] *
  	    Signature[
  	      Transpose[{Join[c, ac, c, bc], l}]
   	   ] *
   	   (SpaceTimeDimension -m + 1)! / (SpaceTimeDimension - m)! *
   	   Det[Outer[DotProduct, ac, bc]]
 	 ] /; Length[{a1, b1}] == Length[{a2, b2}],

    Literal[Power[Eps[a1___, mu, b1___], 2]] :>
	- Block[
	{c = {mu}, m, l, ac = {a1, b1}, a = {a1, mu, b1}},
	m = Length[a];
  	l = Join[Table[1, {m}], Table[2, {m}]];
  	Signature[Transpose[{Join[a, a], l}]] *
  	  Signature[
  	   Transpose[{Join[c, ac, c, ac], l}]
   	   ] *
   	   (SpaceTimeDimension -m + 1)! / (SpaceTimeDimension - m)! *
   	   Det[Outer[DotProduct, ac, ac]]],

    p_[mu] DiracGamma[mu] -> Slash[p],
    p_[mu]q_[mu] -> DotProduct[p, q],

    p_[mu]^2 -> DotProduct[p, p],

    Literal[NonCommutativeMultiply[a___, DiracGamma[mu], b___]] p_[mu] :>
    	NonCommutativeMultiply[a, Slash[p], b],

    Literal[NonCommutativeMultiply[a___, DiracGamma[mu], b___]] *
    	Literal[G[mu, nu_]] :>
    		NonCommutativeMultiply[a, DiracGamma[nu], b],

    Literal[NonCommutativeMultiply[a___, DiracGamma[mu], b___]] *
    	Literal[G[nu_, mu]] :>
    		NonCommutativeMultiply[a, DiracGamma[nu], b],


    Literal[
    	NonCommutativeMultiply[
    		a___, 
    		DiracGamma[mu], 
    		DiracGamma[mu], 
    		b___
    	]
    ] :>
      SpaceTimeDimension NonCommutativeMultiply[a, b],

    Literal[
    	NonCommutativeMultiply[
    		a___, 
    		DiracGamma[mu], 
    		b_, 
    		DiracGamma[mu], 
    		c___
    	]
    ] :>
      (2 - SpaceTimeDimension)NonCommutativeMultiply[a, b, c] /; 
      	GammaLike[b],

    Literal[
    	NonCommutativeMultiply[
    		a___, 
    		DiracGamma[mu], 
    		b_, 
    		c_, 
    		DiracGamma[mu], 
    		d___
    	]
    ] :>
      (SpaceTimeDimension - 4) NonCommutativeMultiply[a, b, c, d] + 
      	4 DotProduct[b, c] NonCommutativeMultiply[a, d] /;
        GammaLike[b, c],

    Literal[
    	NonCommutativeMultiply[
    		a___, 
    		DiracGamma[mu], 
    		b_, 
    		c_, 
    		d_, 
    		DiracGamma[mu], 
    		e___
    	]
    ] :>
      -2 NonCommutativeMultiply[a, d, c, b, e] - 
      (SpaceTimeDimension - 4) NonCommutativeMultiply[a, b, c, d, e] /;
        GammaLike[b, c, d],

    Literal[
    	NonCommutativeMultiply[
    		a___, 
    		DiracGamma[mu], 
    		b1_, b___, 
    		DiracGamma[mu], 
    		c___
    	]
    ] :>
      Block[
        {e1 = OddQ[Length[{b}]], e2 = EvenQ[Length[{a}]]}, 
        If[
          e1,
         2 (NonCommutativeMultiply[a, b, b1, c] + 
	 NonCommutativeMultiply @@ Join[{a}, Reverse[{b, b1}], {c}]),
	-2 NonCommutativeMultiply @@ Join[{a}, Reverse[{b1, b}], {c}]]
        ] /;
      SpaceTimeDimension === 4 && GammaLike[a, b] && Length[{b}] > 3,

    Literal[
      NonCommutativeMultiply[
      	a___, 
      	DiracGamma[mu], 
      	b1_, 
      	b2_, 
      	b3_, 
      	b__, 
      	DiracGamma[mu], 
      	c___
      ]
    ] :>
      Block[
        {m = Length[{b}] + 3, i},
        (-1)^m * (
          (SpaceTimeDimension - 4) *
          	NonCommutativeMultiply[a, b1, b2, b3, b, c] +
          2 NonCommutativeMultiply[a, b3, b2, b1, b, c] -
          2 Sum[
            (-1)^i * (NonCommutativeMultiply @@
              Flatten[{{a}, {b}[[i]], Drop[{b1, b2, b3, b}, {i+3}], {c}}]),
            {i, m-3}
          ]
        )
      ] /;
        GammaLike[b1, b2, b3, b],

(* hey, this rules is useless, it never gets to here does it?  DOn't
the TRaces get evaluated first by GammaTrace??? *)
    Literal[
      GammaTrace[NonCommutativeMultiply[a5___DiracGamma5, a___,
      DiracGamma[mu], b___], opts1___] *
      GammaTrace[NonCommutativeMultiply[c5___DiracGamma5, c___,
      DiracGamma[mu], d___], opts2___]] :>
        (DiracGammaSize/2)(GammaTrace[NonCommutativeMultiply[b, a5, a, d,
	c5, c]] +
	      GammaTrace[NonCommutativeMultiply @@ Join[{b, a5, a}, Reverse[{d,
	      c5, c}]], opts1, opts2]) /; SpaceTimeDimension===4 &&
	      GammaLike[a, b, c, d]
  }



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

(*  AbsSquared  *)

AbsSquared[x_] := Expand[Expand[NonCommutativeExpand[
		AbsSquared1[NonCommutativeFactor[x]]]] //.
		AbsSquaredNonCommutativeMultiplyRules]

AbsSquared1Rules = {
AbsSquared1[x_?(FreeQ[#, NonCommutativeMultiply]&)] :> x Conjugate[x],

AbsSquared1[a_Times] :> AbsSquared1 /@ a,

AbsSquared1[a_Plus] :> 
	Block[
		{d, dd, ddd, i}, 
		dd = Sum[d[i], {i, Length[a]}];
		ddd = Expand[dd Conjugate[dd]];
		ddd /. d[j_] :> a[[j]]
	] /; !FreeQ[a, NonCommutativeMultiply],

AbsSquared1[expr_] :> expr Conjugate[expr]

}

If[$VersionNumber > 1.9,
DownValues[AbsSquared1] = System`Private`ValueList @@ AbsSquared1Rules]

maybeHP[a_] := 1
maybeHP[pol_, a_] := HelicityProjection[a pol]

AbsSquaredNonCommutativeMultiplyRules = {

  Literal[NonCommutativeMultiply[a__, SpinorU[p_, pol1___]]] *
    Literal[NonCommutativeMultiply[SpinorUbar[p_, pol2___], b__]
   ] :>
      NonCommutativeMultiply[
      	a, 
      	maybeHP[pol1, 1], 
      	Slash[p] + Mass[p], 
      	maybeHP[pol2, -1], 
      	b
      ],

  Literal[NonCommutativeMultiply[a__, SpinorV[p_, pol1___]]] *
  Literal[NonCommutativeMultiply[SpinorVbar[p_, pol2___], b__]] :>
  	NonCommutativeMultiply[
  		a, 
		maybeHP[pol1, 1], 
  		Slash[p] - Mass[p], 
  		maybeHP[pol2, -1], 
  		b
  	]

(* Note: these rules are now implemented as permanent rules.  Beware!!!
*)
,
Literal[NonCommutativeMultiply[SpinorUbar[p_], a___, SpinorU[p_]]] :>
    GammaTrace[NonCommutativeMultiply[Slash[p] + Mass[p], a]],

Literal[NonCommutativeMultiply[SpinorVbar[p_], a___, SpinorV[p_]]] :>
    GammaTrace[NonCommutativeMultiply[Slash[p] - Mass[p], a]],

Literal[
	NonCommutativeMultiply[
		SpinorUbar[p_, pol1___], 
		a___, 
		SpinorU[p_, pol2___]
	]
] :>
    GammaTrace[
    	NonCommutativeMultiply[
		maybeHP[pol2, 1],
		Slash[p] + Mass[p], 
		maybeHP[pol1, -1], 
		a
	]],

Literal[
	NonCommutativeMultiply[
		SpinorVbar[p_, pol1___], 
		a___, 
		SpinorV[p_, pol2___]
		]
	] :>
    GammaTrace[
    	NonCommutativeMultiply[
		maybeHP[pol2, 1],
    		Slash[p] - Mass[p], 
    		maybeHP[pol1, -1], 
    		a
    	]
    ]

}


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

(*  GammaTrace  *)

(*
GammaTrace[expr_, opts___] := 
	     Block[{do = (Or @@ (DoEverything /. (List /@ {opts}))) /.
	                         Options[GammaTrace],
	inter = 
	       Trace1[NonCommutativeExpand[ConvertToGamma5[expr]]]},
	       (
	       If[do, inter //. TraceComplicatingRules, inter]
		    /; (do || FreeQ[inter, Trace1] || 
		    inter[[1]] =!= expr))]
*)

GammaTrace[expr_, opts___] := 
	     Block[{od = TraceAllTheWay, 
	inter = (TraceAllTheWay = TraceAllTheWay || 
	                 (Or @@ (DoEverything /. (List /@ {opts}))) /.
	                         Options[GammaTrace];
	       Trace1[NonCommutativeExpand[ConvertToGamma5[expr]]]),
	     do = TraceAllTheWay},
	       (
	       If[do, inter //. TraceComplicatingRules, inter /. Trace1->GammaTrace]
		    /; (TraceAllTheWay = od; (do || FreeQ[inter, Trace1] || 
		    inter[[1]] =!= expr)))]


Trace1Rules = {

Trace1[DiracGamma5] -> 0,

Trace1[_?GammaLike] -> 0,

Trace1[] -> DiracGammaSize,

Trace1[a_Plus] :> (GammaTrace[#, DoEverything->False]& /@ a),

Trace1[a_ b_Plus] :> Trace1[Expand[a b]],

(*  the cutout rules were originally put in for spinor techniques
Trace1[Times[_?SpinorQ, _?SpinorbarQ, DiracGamma5]] -> 0,

Trace1[Times[_?SpinorQ, _?SpinorbarQ, _DiracGamma]] -> 0,

Trace1[Times[_?SpinorQ, _?SpinorbarQ, _Slash]] -> 0,

Trace1[Times[a_?SpinorbarQ, b_?SpinorQ, c___NonCommutativeMultiply]] :>
	NonCommutativeMultiply[a, b] Trace1[c],

*)

Trace1[a_ b_] :> 
	a Trace1[b] /; FreeQ[a, NonCommutativeMultiply] && CommutativeAllQ[a],
 
Trace1[a_ b_] :> 
	a Trace1[b] /; FreeQ[a, NonCommutativeMultiply] && CommutativeQ[a],

Trace1[a_] :> DiracGammaSize a /; FreeQ[a, NonCommutativeMultiply],

Trace1[Literal[NonCommutativeMultiply[_HelicityProjection, b__]]] :>
	GammaTrace[b]/2 /; (Length[{b}] < 4 && GammaLike[b]),

Trace1[Literal[NonCommutativeMultiply[a__]]] :> 
	0 /; OddQ[Length[{a}]] && GammaLike[a],

Trace1[Literal[NonCommutativeMultiply[DiracGamma5, a___]]] :> 0 /;
  (OddQ[Length[{a}]] || Length[{a}] < 4) && GammaLike[a],

Trace1[
        Literal[
                n:NonCommutativeMultiply[
                        a___,
                        DiracGamma[m_],
                        b___,
                        DiracGamma[m_],
                        c___
                ]
        ]
] :>
  If[
        Length[{b}] > Length[{a,c}],
        Trace1[Contract[
                NonCommutativeMultiply[
                        DiracGamma[m],
                        c,
                        a,
                        DiracGamma[m],
                        b
                ], m]
        ],
        Trace1[Contract[n, m]]
      ],

Trace1[ Literal[NonCommutativeMultiply[Slash[p_], b__, Slash[p_]]]] :>
	DotProduct[p,p] Trace1[NonCommutativeMultiply[b]],

(* Almost correct, but might flip sign in traces involving
   Gamma5. Check.

Trace1[
	Literal[NonCommutativeMultiply[a___, Slash[p_], b__, Slash[p_], c___]]
] :>
  Block[
    {br = Reverse[{b}]},
    - DotProduct[p, p] Trace1[NonCommutativeMultiply @@ Join[{a}, br, {c}]] +
      Trace1[NonCommutativeMultiply[a, Slash[p], c]] Trace1[NonCommutativeMultiply @@ Join[{Slash[p]}, br]] / 2 +
      Trace1[NonCommutativeMultiply[a, Slash[p], DiracGamma5, c]] *
        Trace1[NonCommutativeMultiply @@ Join[{DiracGamma5, Slash[p]}, br]] /2
  ] /; SpaceTimeDimension === 4 && OddQ[Length[{b}]] && Length[{b}] > 2,
*)

(*  this has been moved upwards
Trace1[Literal[NonCommutativeMultiply[DiracGamma5, a___]]] :> 0 /;
  (OddQ[Length[{a}]] || Length[{a}] < 4) && GammaLike[a],
*)

Trace1[Literal[NonCommutativeMultiply[DiracGamma5, a_, b_, c_, d_]]] :>
  DiracGammaSize I Eps[a[[1]], b[[1]], c[[1]], d[[1]]] /; GammaLike[a, b, c, d],
(*  
Trace1[
	Literal[NonCommutativeMultiply[DiracGamma5, a_, b_, c_, d_, e_, f_]]
] :>
DiracGammaSize * (Plus @@
	Block[
		{
			l = {a,b,c,d,e,f},
			ii = 
		Union[Select[Sort /@ Flatten[Outer[List, Range[6], Range[6]],
		 			1], UnsameQ @@ #&]], 
		 	g,
		 	eps,
		 	iic
		},
		g = l[[#]]& /@ ii;
		iic = Complement[Range[6], #]& /@ ii;
		eps = l[[#]]& /@ iic;
		(
			Signature[Flatten[#]] * Signature[l] *
			DotProduct[#[[1,1,1]], #[[1,2,1]]] *
			Eps[#[[2,1,1]], #[[2,2,1]], #[[2,3,1]], #[[2,4,1]]]
		)& /@	Transpose[{g, eps}]
	]),
*)


(*  this has been moved 
Trace1[Literal[NonCommutativeMultiply[a__]]] :> 
	0 /; OddQ[Length[{a}]] && GammaLike[a],
*)

Trace1[Literal[NonCommutativeMultiply[a_, b_]]] :>
  DiracGammaSize DotProduct[a, b] /; GammaLike[a, b],

Trace1[Literal[NonCommutativeMultiply[a_, b_, c_, d_]]] :> 
	DiracGammaSize (
  DotProduct[a, b]DotProduct[c, d] + DotProduct[a, d]DotProduct[b, c] -
    DotProduct[a, c]DotProduct[b, d]
    ) /; GammaLike[a, b, c, d]

}

If[$VersionNumber > 1.9,
DownValues[Trace1] = System`Private`ValueList  @@ Trace1Rules]

TraceComplicatingRules = {

(* Check again, for traces involving Gamma5

Trace1[
	Literal[
		NonCommutativeMultiply[a___, Slash[p_], b__, Slash[p_], c___]
	]
] :>
  Block[
    {br = Reverse[{b}], ii = Unique["i"]},
    - DotProduct[p, p] *
    Trace1[NonCommutativeMultiply @@ Join[{a}, br, {c}]] +
      Contract[
        GammaTrace[NonCommutativeMultiply[a, DiracGamma[ii], Slash[p], c]] *
        GammaTrace[NonCommutativeMultiply @@ Join[{Slash[p], DiracGamma[ii]}, br]],
        ii
      ] /2
  ] /; SpaceTimeDimension === 4 && EvenQ[Length[{b}]] && Length[{b}] > 3,
*)

   Literal[Trace1[NonCommutativeMultiply[DiracGamma5, a__]]] :>
   Block[{l = Range[Length[{a}]]},
      I Plus @@ (
      First[#] (Eps @@ First /@ Part[{a}, Last[#]]) *
        GammaTrace[NonCommutativeMultiply @@ Part[{a}, Complement[l, Last[#]]],
	             DoEverything->True]& /@
       SignedCombinations[l, 4])],

Literal[Trace1[NonCommutativeMultiply[a_, b__]]] :>
  Block[
    {j},
    Sum[
      -(-1)^j DotProduct[a, {b}[[j]]] Trace1[Drop[NonCommutativeMultiply[b], {j}]],
      {j, Length[{b}]}
    ]
  ] /; GammaLike[a, b] && Length[{a, b}] > 5
}


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

(*  PhysicsForm1  *)


PhysicsForm1[G[x_, y_]] := SequenceForm["g",Subscript[SequenceForm[x,",",y]]]
PhysicsForm1[DotProduct[p_, q_]] := SequenceForm["(",p,".",q,")"]
PhysicsForm1[DiracGamma[x_]] := SequenceForm["gm", Subscript[x]]
PhysicsForm1[Slash[x_?((Length[#]<2)&)]] := 
	SequenceForm[Subscript["/"], x, Superscript["/"]]
PhysicsForm1[Slash[x_]] := 
	SequenceForm[Subscript["/"], "(", x, ")", Superscript["/"]]
PhysicsForm1[AbsSquared[x_]] := PrecedenceForm[
	SequenceForm["|",x,"|", Superscript[2]], 100]
PhysicsForm1[DiracGamma5] := SequenceForm["gm", Subscript[5]]
PhysicsForm1[SpinorU[x_]] := SequenceForm["u[",x,"]"]
PhysicsForm1[SpinorU[x_, pol_]] := 
	SequenceForm["u",Subscript[PhysicsForm1[pol]],"[",x,"]"]
PhysicsForm1[SpinorUbar[x_]] := ColumnForm[{"_",
	SequenceForm["u[",x,"]"]}, Left, Above]
PhysicsForm1[SpinorUbar[x_, pol_]] := 
ColumnForm[{"_", SequenceForm["u",Subscript[PhysicsForm1[pol]],
	"[",x,"]"]}, Left, Above]
PhysicsForm1[SpinorV[x_]] := SequenceForm["v[",x,"]"]
PhysicsForm1[SpinorV[x_, pol_]] := 
	SequenceForm["v",Subscript[PhysicsForm1[pol]],"[",x,"]"]
PhysicsForm1[SpinorVbar[x_]] := ColumnForm[{"_",
	SequenceForm["v[",x,"]"]}, Left, Above]
PhysicsForm1[SpinorVbar[x_, pol_]] := 
ColumnForm[{"_", SequenceForm["v",Subscript[PhysicsForm1[pol]],
	"[",x,"]"]}, Left, Above]
PhysicsForm1[HelicityProjection[pol_]] := 
    SequenceForm["P",Subscript[PhysicsForm1[pol]]]
PhysicsForm1[Left] := "L"
PhysicsForm1[Right] := "R"
PhysicsForm1[Eps[x__]] := Subscripted[Eps[x]]
PhysicsForm1[Literal[NonCommutativeMultiply[x__]]] :=
	SequenceForm["{", 
	SequenceForm[PhysicsForm1[#], " "]& /@ 
		SequenceForm[x], "}"]
PhysicsForm1[f_[x___]] := PhysicsForm1 /@ f[x]
PhysicsForm1[s_Symbol] := s
PhysicsForm1[c_?NumberQ] := c
PhysicsForm1[] := " "

PhysicsFormPrint[x_] := (Print[PhysicsForm1[x]]; x;)

PhysicsForm[x_] :=
	(
		If[
			ValueQ[$PrePrint],
			(
				temp = $PrePrint; 
				$PrePrint := ($PrePrint = temp; PhysicsForm1[#])&
			),
			$PrePrint := ($PrePrint =.; PhysicsForm1[#])&
		];
		x
	)

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

End[];

EndPackage[]

 
