
(*         Copyright (C) 1990  Eran Yehudai, Alexander C. K. Hsieh          *)

(* Changes after 3/19/91

5/29/91 - Change all ...PhaseSpaceIntegrate to ...PhaseSpaceIntegral

        - Set the output form of PhaseSpaceIntegral to
	"-PhaseSpaceIntegral-"
	
*)


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

LightlikeVectorDecayedFrom::usage = "Take this away for any releases."

Spherical::usage = "Spherical[p1, p2] indicates a spherically
symmetric two body decay into p1 and p2."

Cylindrical::usage = "Cylindrical[p1, p2] indicates a cylindrically
symmetric two body decay into p1 and p2 about the direction of motion
of the decaying particle."

DecayWidth::usage = "DecayWidth[me2, p -> {px, py, pz, e}, outGoing]
returns a PhaseSpaceIntegral. me2 is the matrix element squared, {px,
py, pz, e} is explicit four-vector of the decaying particle, and
outGoing specifies the order of phase-space evaluation."

Options[DecayWidth] = {Mandelstam -> {}}

CrossSection::usage = "CrossSection[me2, p1 -> {p1x, p1y, p1z, e1}, p2
-> {p2x, p2y, p2z, e2}, outGoing] returns a PhaseSpaceIntegral. me2
is the matrix element squared, {p1x, p1y, p1z, e1} and {p2x, p2y, p2z,
e2} are the explicit four-vectors of the incoming particles, and outGoing
specifies the order of phase-space evaluation."

Options[CrossSection] = {Mandelstam -> {}}

PhaseSpaceIntegral::usage = "PhaseSpaceIntegral[integrand, {p1,
(p2)}, decays, integrals, extmom] specifies a phase-space integral.
Use EvaluatePhaseSpaceIntegral to evaluate the integral."

SetAttributes[PhaseSpaceIntegral, HoldRest]

EvaluatePhaseSpaceIntegral::usage =
"EvaluatePhaseSpaceIntegral[PhaseSpaceIntegral[...]] will evaluate
the Phase-space-integral."

NEvaluatePhaseSpaceIntegral::usage =
"NEvaluatePhaseSpaceIntegral[PhaseSpaceIntegral[...]] will evaluate
the Phase-space-integral numerically."

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

Begin["`Private`"]

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

combinations[l_List, k_Integer] := Flatten[
             Block[{i, ll=Length[l],
             m = Table[Unique["t"], {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]

headtest[x_] := MemberQ[{List,Cylindrical,Spherical},x]

createexpr[x_,p_,y__] := Hold[x = Decay[p,y]]

SetAttributes[holding,HoldAll]

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

formParentsChildrenSymmetriesEtc[list_] := Block[
		{p0,variables = {},extmom,list1 = list},
	list1 = list1 //. {a_,b__,c_}->{a,{b,c}};
	list1 = list1 //.
	(h_?(headtest))[a_?(!headtest[Head[#]]&),b_?(!headtest[Head[#]]&)] :>
		Block[{tmp = Unique["ms"]},
		PrependTo[variables,tmp];
		children[tmp] ^= {a,b};
		Mass[tmp] ^= Sqrt[tmp];
		mom[tmp] ^= Unique["q"];
		parent[a] ^= tmp;
		parent[b] ^= tmp;
		Symmetry[tmp] ^= h /. List->None;
		p0 = tmp];
	extmom = Complement[Union[Flatten[children /@ variables]],
                 variables];
	Scan[(sumMassChildren[#] ^= Mass[#])&,extmom];
        Scan[setSiblings,Rest[variables]];
	Scan[setSiblings,extmom];
        Scan[(children[#] ^= #;mom[#] ^= #)&, extmom];
	{variables,extmom}
]


sumMassChildren[x_] := Plus @@ sumMassChildren /@ children[x]

angularIntegrals[z_] := Block[
		{factor, m = Mass /@ children[z]},
		factor = 1/(8*(2 Pi)^6) * 
		    Sqrt[(1-(Plus@@m/Mass[z])^2) (1-(Subtract@@m/Mass[z])^2)];
		factor *= Switch[Symmetry[z],
			None, angles[z] ^= {Unique["x"], Unique["phi"]}; 1,
			Cylindrical, angles[z] ^= {Unique["x"], 0}; 2 Pi,
			Spherical, angles[z] ^= {1,0}; 4 Pi ]
		]

setSiblings[x_] :=  Sibling[x] ^= First[Select[children[parent[x]],(#=!=x)&]]

limits[x_Symbol] := (LowerLimit[x] ^= sumMassChildren[x]^2;
		     UpperLimit[x] ^= (Mass[parent[x]] - Mass[Sibling[x]])^2)

limits[{x_,y_}]  := (LowerLimit[x] ^= sumMassChildren[x]^2;
                     UpperLimit[x] ^= (Mass[parent[x]]-sumMassChildren[y])^2;
                     LowerLimit[y] ^= sumMassChildren[y]^2;
                     UpperLimit[y] ^= (Mass[parent[y]]-Mass[x])^2)

massIntegrals[midmom_] := Scan[limits,
		Rest[Prepend[midmom,0] //.
			{0,b___,a_Symbol,c___,e_Symbol,d___} :> 
			{0,b,{a,e},c,d} /; a == Sibling[e]]]

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

Format[PhaseSpaceIntegral[___]] = "-PhaseSpaceIntegral-"

DecayWidth[me2_, mom1_Rule, list_, opts___] := Block[
		{momentum = mom1[[2]]},
		setUp[(2 Pi)^4/(2 Mass[momentum]) me2, {mom1}, list, opts]]

CrossSection[me2_, mom1_Rule, mom2_Rule, list_, opts___] := Block[
		{momentum1 = mom1[[2]], momentum2 = mom2[[2]]},
		setUp[(2 Pi)^4 /
		(4 Simplify[Sqrt[DotProduct[momentum1,momentum2]^2 - 
		Mass[momentum1]^2 Mass[momentum2]^2]]) me2,
		{mom1,mom2}, list, opts]]
		
(* ************************************************************* *)

setUp[me2_, initmomrules_List, list_, opts___] := Block[
		{psumexplicit = Plus @@ (#[[2]]& /@ initmomrules),
		psumsymbolic = Plus @@ (#[[1]]& /@ initmomrules),
		p0symlist = #[[1]]& /@ initmomrules,
		factor, midmom, extmom,
		decayrules, int, tmp, tmp1, tmp2, newxphilist = {}},
		{midmom,extmom} = formParentsChildrenSymmetriesEtc[list];
		p0 = First[midmom];
		Mass[p0] ^= Mass[psumexplicit];
		mom[p0] ^= psumsymbolic;
		factor = (2 Pi)^(3 (Length[midmom]-1));
		massIntegrals[midmom];
		factor *= Times @@ angularIntegrals /@ midmom;
		decayrules = (Function[{a,ms,b,c},
			Block[{t = (Symmetry[ms] === Cylindrical)},
			If[t,
			createexpr[a,mom[ms],b,c,
				ZAxis -> Momentum[mom[ms]]/.{0,0,0}->{0,0,1}],
			createexpr[a,mom[ms],b,c]]]
			][mom /@ children[#],#,
			angles[#], Mass /@ children[#]])& /@ midmom;
		tmp = Union[Cases[Level[me2, {0,Infinity}],
			LightlikeVectorDecayedFrom[_, 1]]];
		tmp = #[[1]]& /@ tmp;
		If[(tmp1 = Intersection[tmp, p0symlist]) === {}, Null,
			Scan[Block[{temp1, temp2},
			PrependTo[decayrules,
			createexpr[{LightlikeVectorDecayedFrom[#, 1],
			LightlikeVectorDecayedFrom[#, 2]}, #, 
			{temp1 = Unique["x"], temp2 = Unique["phi"]},
			{0, 0}]];
			AppendTo[newxphilist, {temp1, temp2}]]&, tmp1]];
		If[(tmp1 = Complement[tmp, p0symlist]) === {}, Null,
			Scan[Block[{temp1, temp2},
			AppendTo[decayrules,
			createexpr[{LightlikeVectorDecayedFrom[#, 1],
			LightlikeVectorDecayedFrom[#, 2]}, #, 
			{temp1 = Unique["x"], temp2 = Unique["phi"]},
			{0, 0}, ZAxis -> Momentum[#]]];
			AppendTo[newxphilist, {temp1, temp2}]]&, tmp1]];
		decayrules = Apply[CompoundExpression,
			tmp2[Append[decayrules,Null]] /.tmp2 -> holding, 1] /. 
			Hold[c_] :> c;
		int = {factor me2, initmomrules, decayrules,
			Join[{#,LowerLimit[#],UpperLimit[#]}& /@ Rest[midmom],
			{#,-1,1}& /@ Complement[First /@ angles /@midmom,{1}],
			{#,0,2 Pi}& /@ Complement[Last /@angles /@midmom,{0}],
			Flatten[{{#[[1]], -1, 1}, {#[[2]],0, 2 Pi}}& /@ 
			newxphilist, 1]],
			Sequence[Rest[Prepend[extmom,0]
				//.{0,b___,a_Symbol,c___,e_Symbol,d___} :> 
				{0,b,{a,e},c,d} /; a == Sibling[e]], opts]};
		(PhaseSpaceIntegral @@ int /. 
			holding[c_] :> {c}) /.
			{mom[p0] -> psumsymbolic, Mass[p0]->Mass[psumexplicit]}
		]

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

EvaluatePhaseSpaceIntegral[PhaseSpaceIntegral[
	integrand_, initmomrules_List, decayrules_List, integrals_List,
	extmom_List, opts___], opts1___] :=
		Block[{man = Mandelstam/.{opts}/.Options[CrossSection], 
		numerical = First[{opts1, False}],
		manrules, dplist, dprules, tmp, integrand1, dupextmom, mtable,
		rules, integrals1, decayrules1 = Hold[decayrules],
		flatextmom = Flatten[extmom], dpshtcutrules, M,
		psumexplicit = Plus @@ (#[[2]]& /@ initmomrules),
		p0symlist = #[[1]]& /@ initmomrules},
		Scan[(Mass[#[[1]]] ^= Mass[#[[2]]])&, initmomrules];
		mtable = Table[Unique["m"], {Length[flatextmom]}];
		dupextmom = Table[Unique["k"], {Length[flatextmom]}];
		rules = Join[
			Apply[Rule,Cases[
			Transpose[{Mass /@ flatextmom,mtable}],
			{_?(!NumberQ[#]&),_}],1],
			Apply[Rule,Transpose[{flatextmom,dupextmom}],1]];
		rules = If[NumberQ[Mass[psumexplicit]],
			rules,
			Join[rules, {Mass[psumexplicit]^2 -> M^2,
						Mass[psumexplicit] -> M}]];
		dplist = (Append[#,Unique["dp"]])& /@ 
				combinations[Join[flatextmom, p0symlist], 2];
		If[man === {} ||
		(manrules = Select[{DotProduct[#[[1]],#[[2]]],#[[3]]}& /@ 
		dplist, (Head[First[#]]=!=DotProduct)&]) === {}, manrules = {},
		manrules = Flatten[{Reduce[Apply[Equal, manrules, 1], man] /.
				And->List /. Or->List}];
		manrules = ToRules[Union[
				Select[manrules,MemberQ[man, First[#]]&]] /.
				List->And]];
		dpshtcutrules = (DotProduct[#[[1]],#[[2]]] ->
		(Mass[parent[#[[1]]]]^2 -Mass[#[[1]]]^2-Mass[#[[2]]]^2)/2)& /@
			Select[extmom, (MatchQ[#,{_,_}] && 
				(Symmetry[parent[#[[1]]]] === Spherical ||
				Symmetry[parent[#[[1]]]] === Cylindrical))&];
		integrand1  = integrand /. manrules /. dpshtcutrules /. rules;
		decayrules1 = decayrules1 /. rules;
		integrals1  = integrals /. rules;
		dplist = dplist /. rules;
		dprules = (DotProduct[#[[1]],#[[2]]] -> #[[3]])& /@ dplist;
		integrand1 = integrand1 /. dprules;
		Scan[(SetDelayed @@ {#[[3]], DotProduct @@ Drop[#,-1]})&,
			dplist /. initmomrules];
		tmp = Release[decayrules1 /. initmomrules];
		If[integrals1=!={},
			If[numerical === True,
			ClearAttributes[NIntegrate, HoldAll];
			Off[NIntegrate::precwarn];
			tmp = NIntegrate[integrand1, Sequence @@ integrals1];
			SetAttributes[NIntegrate, HoldAll];
			On[NIntegrate::precwarn]; tmp,
			Integrate[integrand1, Sequence @@ integrals1]],
		integrand1] /.
			Apply[Rule,Transpose[{mtable, Mass /@ flatextmom}],1]
			/. M -> Mass[psumexplicit]
		]

NEvaluatePhaseSpaceIntegral[PhaseSpaceIntegral[
	integrand_, initmomrules_List, decayrules_List, integrals_List,
	extmom_List, opts___]] :=
	EvaluatePhaseSpaceIntegral[PhaseSpaceIntegral[
	integrand, initmomrules, decayrules, integrals,
	extmom, opts], True]

End[]

EndPackage[]
