(* Program to calculate the degree of Polarization in Coulomb 
   Scattering as described in Problem 5.18 of QFT. 
   Substitutions are made to give a result as close as possible
   to Eq. (5.217)

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



Remove["Global`*"]

(* Load FeynCalc *)

(* << FeynCalc` *)

(* Definitions *)

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

(* Program *)
Line1:= P[sf]. prop[pf,m] . GS[a] . P[si] . prop[pi,m] . GS[a]  
ans1= TR[Line1];

kin0={SP[a,a]->1, SP[a,pi]->Ei,SP[a,pf]->Ei,SP[pi,pf]->Ei^2 -  \[Beta]^2 Ei^2 Cos[\[Theta]]}//FCI

kin1={SP[a,si]->Ei/m \[Beta] li,SP[a,sf]->Ei/m \[Beta] lf,SP[pi,si]->0,SP[pf,sf]->0,SP[si,sf]->li lf(Ei^2/m^2 \[Beta]^2 -  Ei^2/m^2 Cos[\[Theta]]),SP[pi,sf]->lf Ei^2/m \[Beta] (1-Cos[\[Theta]]),SP[pf,si]->li Ei^2/m \[Beta] (1-Cos[\[Theta]]) }//FCI

kin2={SP[a,si]->0 ,SP[a,sf]->0,SP[pi,si]->0,SP[pf,sf]->0,SP[si,sf]->-li lf,SP[pi,sf]->0,SP[pf,si]-> 0 }//FCI


ans = (ans1 /. kin0) /. kin1
ans2 = (ans1 /. kin0) /. kin2

Amp=Function[{si,sf}, ans /.{li->si,lf->sf,\[Beta]^2->(Ei^2-m^2)/Ei^2}//Simplify]
Amp2=Function[{si,sf}, ans2 /.{li->si,lf->sf,m^2-> Ei^2(1-\[Beta]^2)}//Simplify]


PR = (Amp[1,1]-Amp[1,-1])/(Amp[1,1]+Amp[1,-1])//Simplify 
PR = PR /. Cos[\[Theta]] -> ct // Simplify
PR = PR /. {ct+1-> 2 Cos[\[Theta]/2]^2,ct-1->- 2 Sin[\[Theta]/2]^2}

PR = Simplify[PR, Trig -> False];

Print["1) Result for the degree of Polarization: "]
Print["PR = ",PR,"\n"]

Print["2) Verification of the sum over Polarizations: "]

Print["Amp(RR) = ", Amp[1,1]]
Print["Amp(RL) = ", Amp[1,-1]]
Print["Amp(LR) = ", Amp[-1,1]]
Print["Amp(LL) = ", Amp[-1,-1]]

Aux=Amp[1,1] + Amp[1,-1] + Amp[-1,1] + Amp[-1,-1]//Simplify
Aux=Aux /. Cos[\[Theta]] -> ct // Simplify
Aux=Aux /. {ct+1-> 2 - 2 Sin[\[Theta]/2]^2,ct-1->- 2 Sin[\[Theta]/2]^2}
Aux=Simplify[Aux,Trig->False]
Aux=Simplify[Aux /. m^2-> Ei^2(1-\[Beta]^2),Trig->False]
Aux2=Coefficient[Aux, Ei^2]
Res = Ei^2 Aux2
Print["Eq. (5.11) = \!\(\*UnderscriptBox[\(\[Sum]\), \(spins\)]\)|M\!\(\*SuperscriptBox[\(|\), \(2\)]\) = ",Res,"\n"]


Print["3) Take the polarization axis perpendicular to the scattering plane:  "]

Print["Amp(+,+) = ", Amp2[1,1]]
Print["Amp(+,-) = ", Amp2[1,-1]]
Print["Amp(-,+) = ", Amp2[-1,1]]
Print["Amp(-1,-1) = ", Amp2[-1,-1]]

Aux3=Amp2[1,1] + Amp2[1,-1] + Amp2[-1,1] + Amp2[-1,-1]//Simplify
Aux3=Aux3 /. Cos[\[Theta]] -> 1 - 2 Sin[\[Theta]/2]^2 
Aux4=Simplify[Aux3,Trig->False]
Aux4=Simplify[Aux4 /. m^2-> Ei^2(1-\[Beta]^2),Trig->False]
Aux5=Coefficient[Aux4, Ei^2]
Res2 = Ei^2 Aux5
Print["Eq. (5.11) = \!\(\*UnderscriptBox[\(\[Sum]\), \(spins\)]\)|M\!\(\*SuperscriptBox[\(|\), \(2\)]\) = ",Res2,"\n"]




(********************* End of Program **********************)
