(* Program to Calculate the Traces of the decay of positronium
   as in Problem 5.19 of ITC

 Last Version: 06/12/2024
 Compatible with FeynCalc 10.1.0
 Author: Jorge C. Romao
 email: jorge.romao@tecnico.ulisboa.pt
*)

Remove["Global`*"]

(* Definitions *)

prop[p_,m_]:=GS[p]+m
Ps[s_]:=(1 + GA[5] . GS[s])/2


  
Line1:= prop[p2,-m] .  GA[mu] . prop[p1-k1,m] . GA[nu] . Ps[s1] . prop[p1,m] . GA[nu] . prop[p1-k1,m] . GA[mu] . Ps[s2]

Line2:= prop[p2,-m] . GA[nu] . prop[p1-k2,m] . GA[mu] . Ps[s1] . prop[p1,m] . GA[mu] . prop[p1-k2,m] . GA[nu] . Ps[s2]

Line12:= prop[p2,-m] . GA[mu] . prop[p1-k1,m] . GA[nu] . Ps[s1] . prop[p1,m] . GA[mu] . prop[p1-k2,m] . GA[nu] . Ps[s2]

Line21:= prop[p2,-m] .  GA[nu] . prop[p1-k2,m] . GA[mu] . Ps[s1] . prop[p1,m] .  GA[nu] . prop[p1-k1,m] . GA[mu] . Ps[s2]


ans1= Simplify[Contract[TR[Line1]]/(t-m^2)^2]
ans2= Simplify[Contract[TR[Line2]]/(u-m^2)^2]
ans12= Simplify[Contract[TR[Line12]]/((t-m^2) (u-m^2))]
ans21= Simplify[Contract[TR[Line21]]/((t-m^2) (u-m^2))]


ans=  (ans1 + ans2 + ans12 + ans21)//Calc //FCE


onshell={SP[p1,p1]->m^2,SP[p2,p2]->m^2,SP[k1,k1]->0,SP[k2,k2]->0,SP[s1,s1]->-1,SP[s2,s2]->-1,SP[p1,s1]->0,SP[p2,s2]->0}


res = ans /. onshell  (* 1/2 for identical particles *)


kin={SP[p1,p2]->(s-2 m^2)/2 ,SP[p1,k1]->(m^2-t)/2 ,SP[p1,k2]->(m^2-u)/2,SP[p2,k1]->(m^2-u)/2,SP[p2,k2]->(m^2-t)/2,SP[k1,k2]->s/2}

kin2={SP[s1,s2]->l1 l2 s/4/m^2 (beta^2+1), SP[p1,s2]->s/2/m beta l2, SP[p2,s1]->s/2/m l1 beta, SP[k1,s1]-> s/4/m l1 (beta - x),SP[k1,s2]-> s/4/m l2 (beta + x), SP[k2,s1]-> s/4/m l1 (beta + x), SP[k2,s2]-> s/4/m l2 (beta - x)}


aux=Expand[(res /. kin)/. kin2]


final = Expand[aux /. {t-> Dt+m^2,u->Du+m^2}]


final2 = final /. {Du -> -s/2 (1 + beta x), Dt -> -s/2 (1 - beta x)}


xs = Integrate[final2, {x, -1, 1}, Assumptions -> {m > 0, s > 0, beta > 0, beta<=1}]

(* 1/2 for identical particles *)
xs = xs/32/Pi/s/(s-4 m^2)^(1/2)*(s)^(1/2) * (4 Pi \[Alpha])^2 * (1/2) 
xs = xs /. beta-> (s-4 m^2)^(1/2)/s^(1/2) // Expand // Simplify



XS=Function[{s1,s2}, xs /. {l1->s1,l2->s2}];

xsNonPol=(XS[1, 1] + XS[1, -1] + XS[-1, 1] + XS[-1, -1])/4 // Simplify


Print["1) Non-Polarized Cross Section "]
Print["\[Sigma] = ",xsNonPol]

Print["2) Singlet Cross Section "]
Print["\*SubscriptBox[\(\[Sigma]\), singlet] = ",4 xsNonPol]

Print["3) Non-Relativistic Limit of Singlet Cross Section:      ", s \[TildeEqual] 4 m^2 + m^2 \[Beta]^2]

xsNonPolNRL=Simplify[Normal[Series[(4 xsNonPol /. s-> 4 m^2 + m^2 \[Beta]^2),{\[Beta],0,0}]],Assumptions->{m>0,\[Beta]>0}];

Print["\*SubscriptBox[\(\[Sigma]\), singlet] \[TildeEqual] ", xsNonPolNRL]

(********************* Fim do Programa **********************)

