BeginPackage["SpinorTechniques`", "CrossSection`",
"RelativisticKinematics`", "Spinor`", "NonCommutativeMultiply`"]

Clear[Squared, Opposite]

Clear[SpinorS, SpinorT, ConvertToST, STToTraces,
LightlikeVectorCollinearWith, PolarizationCombinations]

LightlikeVectorDecayedFrom::usage = "LightlikeVectorDecayedFrom[p, n]
represents a lightlike vector used for the polarizations of massive
vector bosons.  It obeys p = LightlikeVectorDecayedFrom[p,1] +
LightlikeVectorDecayedFrom[p,2] where p is the momentum of the
external massive vector boson."

Opposite::usage = "Opposite[p] represents a four-momentum which has
the same energy but opposite three momentum to the four-momentum p."

Squared::usage = "Squared[me] indicates to
CEvaluatePhaseSpaceIntegrate that the matrix-element needs to be
squared."

SetAttributes[Squared, Listable]

PolarizationCombinations::usage =
"PolarizationCombinations[me] returns a list of the possible
polarization combinations of the initial and final particles for the
matrix-element me over the polarizations which are not explicitly
written in the matrix-element."

LightlikeVectorNotCollinearWith::usage = "LightlikeVectorNotCollinearWith[p]
represents a lightlike vector that is not collinear with p."

SpinorS::usage = "SpinorS[p, q] is the same as
SpinorUbar[p, Right]**SpinorU[q, Left]. Both p and q have to be
massless."

SpinorT::usage = "SpinorT[p, q] is the same as
SpinorUbar[p, Left]**SpinorU[q, Right]. Both p and q have to be
massless."

ConvertToST::usage = "ConvertToST[expr] attempts to convert DiracGamma5
and HelicityProjection occurrences in expr to SpinorS and SpinorT
objects."

Options[ConvertToST] = {ConvertToMassless -> {}}

ConvertToMassless::usage = "ConvertToMassless is an option for ConvertToST."

STToTraces::usage = "STToTraces[expr] attempts to convert products of
SpinorS's and SpinorT's in expr to traces using the rule
SpinorS[p1,p2]*SpinorT[p2, p3]...SpinorS[pn, p1] =
GammaTrace[Slash[p1]**Slash[p2]**...**Slash[pn]**HelicityProjection[Right]]"

SetAttributes[STToTraces, Listable]



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


Begin["SpinorTechniques`Private`"]

Clear[tr, SpinorST]

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

LeftRightToNumber[Left] = -1
LeftRightToNumber[Right] = 1

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

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

(*  PolarizationCombinations  *)

PolarizationCombinations[me_] := Block[
	{tmp = Level[me, {0,Infinity}], u, ubar, v, vbar, mless, rules1,
	rules2, rules},
	u = Union[Cases[tmp, SpinorU[_]]];
	ubar = Union[Cases[tmp, SpinorUbar[_]]];
	v = Union[Cases[tmp, SpinorV[_]]];
	vbar = Union[Cases[tmp, SpinorVbar[_]]];
	mless = Union[Cases[tmp, MasslessVectorPolarization[_,_]]];

	MasslessInAndOutList = Join[u, ubar, v, vbar, mless];
	rules1 = {Insert[#, Left, -1], Insert[#, Right, -1]}& /@ 
		MasslessInAndOutList;
	rules1 = Flatten[Outer[List, Sequence @@ rules1],
		Length[MasslessInAndOutList]-1];
	rules1 = Apply[Rule,
		Transpose[{MasslessInAndOutList, #}]& /@ rules1, {2}];
	rules2 = {Insert[#, Right, -1], Insert[#, Left, -1]}& /@ 
		MasslessInAndOutList;
	rules2 = Flatten[Outer[List, Sequence @@ rules2],
		Length[MasslessInAndOutList]-1];
	rules2 = Apply[Rule,
		Transpose[{#, Table[0, {Length[MasslessInAndOutList]}]}]& /@ 
		rules2, {2}];
	rules = Flatten /@ Transpose[{rules1, rules2}];
	(me /. #)& /@ rules]

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

(*  PhysicsForm1  *)

(*  these won't work because they go AFTER the /@ rules for PhysicsForm1 *)

Spinor`Private`PhysicsForm1[
	MasslessVectorPolarization[p_, mu_, pol___]] :=
	SequenceForm["eps",ColumnForm[{mu," ",
	Spinor`Private`PhysicsForm1[pol]}, Left, Center], "[", p, "]"]

Spinor`Private`PhysicsForm1[
	HeavyVectorPolarization[p_, mu_, pol___]] :=
	SequenceForm["eps",ColumnForm[{mu," ",
	Spinor`Private`PhysicsForm1[pol]}, Left, Center], "[", p, "]"]

Spinor`Private`PhysicsForm1[SpinorS[p_, q_]] := 
    	SequenceForm["s[", p, ",", q, "]"]

Spinor`Private`PhysicsForm1[SpinorT[p_, q_]] := 
    	SequenceForm["t[", p, ",", q, "]"]

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

(* Squared *)

Squared[0] = 0

Squared[] = 1

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

(* Complex rules *)

ComplexQ[_SpinorS] = True

ComplexQ[_SpinorT] = True

ComplexQ[_DotProduct] = False

ComplexQ[_Squared] = False


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

(* Conjugation rules *)


Unprotect[Abs]
Abs/: Conjugate[x_Abs] := x
Protect[Abs]

SpinorS/: Conjugate[Literal[SpinorS[p1_, p2_]]] := 
	  -SpinorT[Conjugate[p1], Conjugate[p2]]

SpinorT/: Conjugate[Literal[SpinorT[p1_, p2_]]] := 
	  -SpinorS[Conjugate[p1], Conjugate[p2]]

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

LightlikeVectorDecayedFrom/: Conjugate[
	LightlikeVectorDecayedFrom[p_, n_]] := 
	   LightlikeVectorDecayedFrom[p,n]


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

(*  SpinorS & SpinorT rules  *)

Literal[SpinorS[p_, p_]] = 0

Literal[SpinorT[p_, p_]] = 0

Literal[SpinorST[__, p_, p_, __]] = 0

Literal[SpinorST[pol1_, p1_, p___, pn_, pol2_]] := 0 /;
		(pol1 == pol2 && EvenQ[Length[{p}]]) ||
		(pol1 == -pol2 && OddQ[Length[{p}]])

Literal[SpinorS[p1_, p2_]] := -SpinorS[p2, p1] /; !OrderedQ[{p1, p2}]

Literal[SpinorT[p1_, p2_]] := -SpinorT[p2, p1] /; !OrderedQ[{p1, p2}]


SpinorS/: Literal[SpinorS[p1_, p2_]] * Literal[SpinorT[p1_, p2_]] := 
    -2 DotProduct[p1, p2]

SpinorS/: SpinorS[p1_, p2_] * Power[SpinorT[p1_, p2_], n_] :=
	-2 DotProduct[p1,p2] SpinorT[p1, p2]^(n-1) /; n > 1

SpinorT/: SpinorT[p1_, p2_] * Power[SpinorS[p1_, p2_], n_] :=
	-2 DotProduct[p1,p2] SpinorS[p1, p2]^(n-1) /; n > 1


SpinorT[p1:{_, _, _, _}, p2:{_, _, _, _}] := Conjugate[SpinorS[p2, p1]]

SpinorS[{px_, py_, pz_, pt_}, {kx_, ky_, kz_, kt_}] :=
	(py + I pz) Sqrt[(kt - kx)/(pt - px)] - 
	(ky + I kz) Sqrt[(pt - px)/(kt - kx)]

SpinorST/: Literal[SpinorST[pol1_, p___, pn_, pol2_] *
   SpinorST[pol2_, pn_, pp___, pol3_]] :=
   SpinorST[pol1, p, pn, pp, pol3]

SpinorST/: Literal[SpinorST[pol1_, p1_, p___, pol2_] *
   SpinorST[pol3_, p1_, pp___, pol4_]] :=
   LeftRightToNumber[pol1] LeftRightToNumber[pol2] *
   SpinorST[-pol2, ReverseInPlace[p], p1, pp, pol4] /; pol1 == -pol3

SpinorST/: Literal[SpinorST[pol1_, p___, p1_, pol2_] *
   SpinorST[pol3_, pp___, p1_, pol4_]] :=
   LeftRightToNumber[pol3] LeftRightToNumber[pol4] *
   SpinorST[pol1, p, p1, ReverseInPlace[pp], -pol3] /; pol4 == -pol2

SpinorST/: Literal[SpinorST[_, p1_, p___, pn_, pol_]] *
   Literal[Conjugate[SpinorST[_, p1_, p___, pn_, pol_]]] :=
   GammaTrace[HelicityProjection[(-1)^Length[{p}] pol] ** 
   (NonCommutativeMultiply @@ 
   (Slash /@ {p, pn, ReverseInPlace[p], p1}))]

SpinorST/: SpinorST[pol_, p1_, p___, p1_, pol_] :=
   GammaTrace[HelicityProjection[pol] ** (NonCommutativeMultiply @@
   (Slash /@ {p1, p}))]

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

(*  Converting S,T expr to Traces  *)

STToTraces[x_Plus] := STToTraces /@ x

STToTraces[a_ b_] := a STToTraces[b] /; 
            FreeQ[a, SpinorS] && FreeQ[a, SpinorT] && FreeQ[a, SpinorST]

STToTraces[expr_] := expr //.
{
SpinorS[p1_, p2_]^n1_. * SpinorT[p1_, p2_]^n2_. :>
	(-2 DotProduct[p1,p2])^n1 SpinorT[p1, p2]^(n2-n1) /; (n1>0 && n2>=n1),

SpinorS[p1_, p2_]^n1_. * SpinorT[p1_, p2_]^n2_. :>
	(-2 DotProduct[p1,p2])^n2 SpinorS[p1, p2]^(n1-n2) /; (n2>0 && n1>=n2),

SpinorS[p1_, p2_] :> SpinorST[Right, p1, p2, Left],

SpinorT[p1_, p2_] :> SpinorST[Left, p1, p2, Right],

SpinorST[pol1_, p1___, p2_, pol2_]^n1_. *
	SpinorST[pol2_, p2_, p3___, pol3_]^n2_. :>
	SpinorST[pol1, p1, p2, p3, pol3]^n1 * 
	SpinorST[pol2, p2, p3, pol3]^(n2-n1) /; (n1>0 && n2>=n1),

SpinorST[pol1_, p1___, p2_, pol2_]^n1_. *
	SpinorST[pol2_, p2_, p3___, pol3_]^n2_. :>
	SpinorST[pol1, p1, p2, p3, pol3]^n2 * 
	SpinorST[pol1, p1, p2, pol2]^(n1-n2) /; (n2>0 && n1>=n2),


SpinorST[pol1_, p1_, p___, pol2_]^n1_. *
	SpinorST[pol3_, p1_, pp___, pol4_]^n2_. :>
   	(LeftRightToNumber[pol1] LeftRightToNumber[pol2] *
	SpinorST[-pol2, ReverseInPlace[p], p1, pp, pol4])^n1 *
	SpinorST[pol3, p1, pp, pol4]^(n2-n1) /;
	(pol1 == -pol3 && n1>0 && n2>=n1),

SpinorST[pol1_, p1_, p___, pol2_]^n1_. *
	SpinorST[pol3_, p1_, pp___, pol4_]^n2_. :>
   	(LeftRightToNumber[pol1] LeftRightToNumber[pol2] *
	SpinorST[-pol2, ReverseInPlace[p], p1, pp, pol4])^n2 *
	SpinorST[pol1, p1, p, pol2]^(n1-n2) /;
	(pol1 == -pol3 && n2>0 && n1>=n2),

SpinorST[pol1_, p___, p1_, pol2_]^n1_. *
	SpinorST[pol3_, pp___, p1_, pol4_]^n2_. :>
   	(LeftRightToNumber[pol3] LeftRightToNumber[pol4] *
	SpinorST[pol1, p, p1, ReverseInPlace[pp], -pol3])^n1 *
	SpinorST[pol3, pp, p1, pol4]^(n2-n1) /;
	(pol4 == -pol2 && n1>0 && n2>=n1),

SpinorST[pol1_, p___, p1_, pol2_]^n1_. *
	SpinorST[pol3_, pp___, p1_, pol4_]^n2_. :>
   	(LeftRightToNumber[pol3] LeftRightToNumber[pol4] *
	SpinorST[pol1, p, p1, ReverseInPlace[pp], -pol3])^n2 *
	SpinorST[pol1, p, p1, pol2]^(n1-n2) /;
	(pol4 == -pol2 && n2>0 && n1>=n2),

SpinorST[pol1_, p1_, p___, pn_, pol2_]^n1_. *
	Conjugate[SpinorST[pol1_, p1_, p___, pn_, pol2_]]^n2_. :>
	(GammaTrace[HelicityProjection[(-1)^Length[{p}] pol2] ** 
        (NonCommutativeMultiply @@
	(Slash /@ {p, pn, ReverseInPlace[p], p1}))])^n1 *
	Conjugate[SpinorST[pol1, p1, p, pn, pol2]]^(n2-n1) /; (n1>0 && n2>=n1),

SpinorST[pol1_, p1_, p___, pn_, pol2_]^n1_. *
	Conjugate[SpinorST[pol1_, p1_, p___, pn_, pol2_]]^n2_. :>
	(GammaTrace[HelicityProjection[(-1)^Length[{p}] pol2] ** 
        (NonCommutativeMultiply @@
	(Slash /@ {p, pn, ReverseInPlace[p], p1}))])^n2 *
	SpinorST[pol1, p1, p, pn, pol2]^(n1-n2) /; (n2>0 && n1>=n2)

}


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

(* Converting to ST expression *)

LightlikeVectorNotCollinearWith/:  Mass[_LightlikeVectorNotCollinearWith] = 0

LightlikeVectorDecayedFrom/:  Mass[_LightlikeVectorDecayedFrom] = 0

LightlikeVectorDecayedFrom/:
	DotProduct[LightlikeVectorDecayedFrom[p_, 1], 
		LightlikeVectorDecayedFrom[p_, 2]] := Mass[p]^2/2

ConvertToST[expr_, opts___] := Block[
	{x = ConvertToMassless /. {opts} /. Options[ConvertToST]},
	x = Switch[x,
		{}, {},
		{_Rule..},
			Flatten[Block[{a = #[[1]], b = #[[2]]},
			{Slash[a] -> ExpandSlash[Slash[b]],
			a[mu_] :> Release[b[mu]]}]& /@ x]];
	ConvertToST1[Expand[ExpandSlash[expr]], x]]


ConvertToST1[a_Plus, rules_] := ConvertToST1[#, rules]& /@ a

ConvertToST1[a_, rules_] := Block[
	{tmp = a /. rules //. ConvertToSTRules},
	If[tmp === a,  tmp,  ConvertToST1[Expand[ExpandSlash[tmp]], rules]]]

ConvertToSTRules = {
	SpinorV[p_, pol_] :> SpinorU[p, pol] /; Mass[p] == 0,

	SpinorVbar[p_, pol_] :> SpinorUbar[p, pol] /; Mass[p] == 0,

	Literal[NonCommutativeMultiply[a___, b_Plus, c___]] :>
		(NonCommutativeMultiply[a, #, c]&) /@ b,

	Literal[NonCommutativeMultiply[
	___, SpinorUbar[_, pol_],SpinorU[_, pol_], ___]] :> 0,

	Literal[ncm:NonCommutativeMultiply[___, DiracGamma[mu_], ___]] *
	Literal[pp:Plus[(_[___, mu_,___])..]] :> Expand[ncm pp],

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

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

	Literal[NonCommutativeMultiply[a___, Slash[p_], b___]] :>
		NonCommutativeMultiply[
		a, SpinorU[p, Left], SpinorUbar[p, Left], b] +
		NonCommutativeMultiply[
		a, SpinorU[p, Right], SpinorUbar[p, Right], b] /;Mass[p] === 0,

	Literal[NonCommutativeMultiply[a___, Slash[p_], b___]] :>
		NonCommutativeMultiply[a,
		Slash[LightlikeVectorDecayedFrom[p, 1]] +
		Slash[LightlikeVectorDecayedFrom[p, 2]], b] /; Mass[p] =!= 0,

	MasslessVectorPolarization[p_, mu_, pol_] :> 
		Block[{tmp = LightlikeVectorNotCollinearWith[p]},
		SpinorUbar[p, pol] ** DiracGamma[mu] ** SpinorU[tmp, pol] / 
		(2 Sqrt[Abs[DotProduct[p, tmp]]])],

	HeavyVectorPolarization[p_, mu_] :>
		SpinorUbar[LightlikeVectorDecayedFrom[p, 1], Left] **
		DiracGamma[mu] **
		SpinorU[LightlikeVectorDecayedFrom[p, 2], Left] *
		Sqrt[3/(2 Pi)]/(2 Mass[p]),

(*  I think the next rule is only valid for massless p1 & p2  *)
	Literal[NonCommutativeMultiply[a___, DiracGamma[mu_], b___]] *
	Literal[NonCommutativeMultiply[
	SpinorUbar[p1_, pol_], DiracGamma[mu_], SpinorU[p2_, pol_]]] :>
		2 (NonCommutativeMultiply[
		a, SpinorU[p2, pol], SpinorUbar[p1, pol], b] +
		NonCommutativeMultiply[
		a, SpinorU[p1, -pol], SpinorUbar[p2, -pol], b]) /;
		Mass[p1] == 0 && Mass[p2] == 0,

	Literal[NonCommutativeMultiply[a___, SpinorUbar[p1_, pol1___],
	SpinorU[p2_, pol2___], b___]] :> SpinorT[p1, p2] *
		NonCommutativeMultiply[a] NonCommutativeMultiply[b] /; 
		Mass[p1] == 0 && Mass[p2] == 0 && 
		((pol1 == Left && pol2 == Right) || 
		(pol1 == Left && {pol2} === {}) ||
		({pol1} === {} && pol2 == Right)),

	Literal[NonCommutativeMultiply[a___, SpinorUbar[p1_, pol1___],
	SpinorU[p2_, pol2___], b___]] :> SpinorS[p1, p2] *
		NonCommutativeMultiply[a] NonCommutativeMultiply[b] /; 
		Mass[p1] == 0 && Mass[p2] == 0 && 
		((pol1 == Right && pol2 == Left) || 
		(pol1 == Right && {pol2} === {}) ||
		({pol1} === {} && pol2 == Left)),

(* the rule also appears in STToTraces *)
SpinorS[p1_, p2_]^n1_. * SpinorT[p1_, p2_]^n2_. :>
	(-2 DotProduct[p1,p2])^n1 SpinorT[p1, p2]^(n2-n1) /; (n1>0 && n2>=n1),


SpinorS[p1_, p2_]^n1_. * SpinorT[p1_, p2_]^n2_. :>
	(-2 DotProduct[p1,p2])^n2 SpinorS[p1, p2]^(n1-n2) /; (n1>=n2 && n2>0),

(* make sure these rules come last *)
Eps[a_, b_, c_, d_] :> 
	- I/4 tr[NonCommutativeExpand[NonCommutativeMultiply[
	DiracGamma[a], DiracGamma[b], DiracGamma[c], DiracGamma[d],
	HelicityProjection[Right] - HelicityProjection[Left]]]],

Literal[tr[NonCommutativeMultiply[a___, b_, c___]]] :>
	tr[NonCommutativeExpand[NonCommutativeMultiply[a, Slash[b], c]]] /;
	Mass[b] === 0

}

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

(* rules about trace, assumes SpaceTimeDimension is 4 *)

tr[] = 4

tr[a_Plus] := tr /@ a

tr[a_ b_] := a tr[b] /; 
	NumberQ[a] || Head[a] === SpinorS || Head[a] === SpinorT

tr[SpinorUbar[p1_, pol_] SpinorU[p2_, -pol_]] := 
	4 NonCommutativeMultiply[SpinorUbar[p2, pol], SpinorU[p1, -pol]]

tr[a_] := 4 a /;
	FreeQ[a, NonCommutativeMultiply] && FreeQ[a, SpinorU] &&
	FreeQ[a, SpinorUbar]

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

End[];

EndPackage[]
