(* Mathematica program for Problem 1.9

  date : 27/10/2009
  author: Jorge C. RomÃ£o
  email : jorge.romao@ist.utl.pt
*)

(* Initial Data *)
Clear[E1,E3,E4,m,a,teta]

E1=100.;  (* in MeV *)
m=0.5;    (* in MeV *)
a=1;      (* head on collision *)
teta=Pi;  (* back in the CM system *)

(* CM Kinematics *)
p=Function[E,Sqrt[E^2-m^2]];

beta=Function[E2,(p[E1]-a*E2)/(E1+E2)];
gama=Function[E2,1/Sqrt[1-beta[E2]^2]];
s=Function[E2,(E1+E2)^2-(p[E1]-a*E2)^2];

E3cm=Function[E2,(s[E2]+m^2)/2/Sqrt[s[E2]]];
E4cm=Function[E2,(s[E2]-m^2)/2/Sqrt[s[E2]]];
pcm=Function[E2,E4cm[E2]];

(* Lab Kinematics *)

E3=Function[E2,gama[E2]*(E3cm[E2]+beta[E2]*(pcm[E2]*Cos[teta]))];
E4=Function[E2,gama[E2]*(E4cm[E2]+beta[E2]*(-pcm[E2]*Cos[teta]))];

p3=Function[E2,gama[E2]*(pcm[E2]*Cos[teta]+beta[E2]*E3cm[E2])];
p4=Function[E2,gama[E2]*(-pcm[E2]*Cos[teta]+beta[E2]*E4cm[E2])];

(* Energy loss *)

DeltaE=Function[E2,E1-E3[E2]];
Print["Energy Loss = ",DeltaE[4.14*10^(-10)], " MeV"]

(* Plot of Energy Loss *)
LogLogPlot[E1-E3[E2],{E2,10^(-10),100},PlotRange->{{10^(-10),100},
           {10^(-6),200}},AxesLabel->{"Eγ (MeV)","ΔE (MeV)"}]