Exit[]

Neutral Pion Photoproduction

List of data:
	ANDERSON 1970 http://hepdata.cedar.ac.uk/view/ins67154
	BRAUNSCHWEIG 1973  http://hepdata.cedar.ac.uk/view/ins84302
	BRAUNSCHWEIG 1973  http://hepdata.cedar.ac.uk/view/ins84301
	BOOTH 1972 http://hepdata.cedar.ac.uk/view/ins75858
	BELLENGER 1969 http://hepdata.cedar.ac.uk/view/ins55099
	OSBORNE 1972 http://hepdata.cedar.ac.uk/view/ins84571
	BRAUNSCHWEIG 1968 http://hepdata.cedar.ac.uk/view/red7610
	BOLON 1967 http://hepdata.cedar.ac.uk/view/ins51286
	BARISH 1974 http://hepdata.cedar.ac.uk/view/ins94529
	GOMEZ 1971 http://hepdata.cedar.ac.uk/view/p5186
	DEUTSCH 1972 http://hepdata.cedar.ac.uk/view/ins84562

MP=0.938;MPI=0.138;
foldername = "/Users/vincent/Dropbox/mes articles/Regge/data";
SetDirectory[foldername];
Needs["PlotLegends`"]
Needs["ErrorBarPlots`"]
Print["Need to change the variable foldername to the folder where you have the data files"]
Print["Need to install the package ErrorBarLogPlots from http://library.wolfram.com/infocenter/MathSource/6747/"]
Needs["ErrorBarLogPlots`"]

January 2015

Execute the commands in the order

Prepare data for fit

dat2fit := Block[{sig2,sig3,sig4,sig5,sig58,
sig6g,sig18g,sig12g,sig6,sig9,sig12,sig15,dat,
datfit1,datfit},
(* create a file with the entries 
    Eγ  t  dσ/dt  Δdσ/dt *)
SetDirectory[foldername];
datfit = {}; (* initialization *)

(* BRAUNSCHWEIG 1968 http://hepdata.cedar.ac.uk/view/red7610 *)
sig2 = Import["Braun68_Sig_PhotPi0_2.txt","Table"];
sig3 = Import["Braun68_Sig_PhotPi0_3.txt","Table"];
sig4 = Import["Braun68_Sig_PhotPi0_4.txt","Table"];
sig5 = Import["Braun68_Sig_PhotPi0_5.txt","Table"];
sig58 = Import["Braun68_Sig_PhotPi0_58.txt","Table"];

dat = sig2;
datfit1 = Table[{2,-dat[[i,4]],dat[[i,10]],dat[[i,11]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];
dat = sig3;
datfit1 = Table[{3,-dat[[i,4]],dat[[i,10]],dat[[i,11]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];
dat = sig4;
datfit1 = Table[{4,-dat[[i,4]],dat[[i,10]],dat[[i,11]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];
dat = sig5;
datfit1 = Table[{5,-dat[[i,4]],dat[[i,10]],dat[[i,11]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];
dat = sig58;
datfit1 = Table[{5.8,-dat[[i,4]],dat[[i,10]],dat[[i,11]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];

(* GOMEZ 1971 http://hepdata.cedar.ac.uk/view/p5186 *)
sig6g   = Import["Gomez71_Sig_PhotPi0_6.txt","Table"];
sig12g = Import["Gomez71_Sig_PhotPi0_12.txt","Table"];
sig18g   = Import["Gomez71_Sig_PhotPi0_18.txt","Table"];

dat = sig6g;
datfit1 = Table[{6,-dat[[i,1]],10^3dat[[i,4]],10^3dat[[i,5]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];
dat = sig12g;
datfit1 = Table[{12,-dat[[i,1]],10^3dat[[i,4]],10^3dat[[i,5]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];
dat = sig18g;
datfit1 = Table[{18,-dat[[i,1]],10^3dat[[i,4]],10^3dat[[i,5]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];

(* ANDERSON 1970 http://hepdata.cedar.ac.uk/view/ins67154 *)
sig6 = Import["Anderson_SigPhotPi0_6.txt","Table"];
sig9 = Import["Anderson_SigPhotPi0_9.txt","Table"];
sig12 = Import["Anderson_SigPhotPi0_12.txt","Table"];
sig15 = Import["Anderson_SigPhotPi0_15.txt","Table"];

dat = sig6;
datfit1 = Table[{6,dat[[i,1]],dat[[i,4]],dat[[i,5]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];
dat = sig9;
datfit1 = Table[{9,dat[[i,1]],dat[[i,4]],dat[[i,5]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];
dat = sig12;
datfit1 = Table[{12,dat[[i,1]],dat[[i,4]],dat[[i,5]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];
dat = sig15;
datfit1 = Table[{15,dat[[i,1]],dat[[i,4]],dat[[i,5]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];

datfit (* return the table *)
];


dataBA:=Block[{BeamA4,BeamA6,BeamA10,dat,datfit,datfit1},
SetDirectory[foldername];

datfit = {}; (* initialization *)

(* ANDERSON 1970 http://hepdata.cedar.ac.uk/view/ins67154 *)
BeamA4 = Import["Anderson_BeamAPhotPi0_4.txt","Table"];
BeamA6 = Import["Anderson_BeamAPhotPi0_6.txt","Table"];
BeamA10 = Import["Anderson_BeamAPhotPi0_10.txt","Table"];

dat = BeamA4;
datfit1 = Table[{4,dat[[i,1]],dat[[i,4]],dat[[i,5]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];
dat = BeamA6;
datfit1 = Table[{6,dat[[i,1]],dat[[i,4]],dat[[i,5]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];
dat = BeamA10;
datfit1 = Table[{10,dat[[i,1]],dat[[i,4]],dat[[i,5]]},{i,Length[dat]}];
datfit = Join[datfit,datfit1];

datfit (* return the table *)
];

For the plots

SetDirectory[foldername];
Braun68:=Module[{sig2,sig3,sig4,sig5,sig58,dat,
ps2,ps3,ps4,ps5,ps6,ps},
(* load data *)
sig2 = Import["Braun68_Sig_PhotPi0_2.txt","Table"];
sig3 = Import["Braun68_Sig_PhotPi0_3.txt","Table"];
sig4 = Import["Braun68_Sig_PhotPi0_4.txt","Table"];
sig5 = Import["Braun68_Sig_PhotPi0_5.txt","Table"];
sig58 = Import["Braun68_Sig_PhotPi0_58.txt","Table"];
ScaleFacB68 = {10^4,10^3,10^2,10^1,10^0};

dat = sig2;
ps2 = ErrorListLogPlot[
Table[{{dat[[i,4]],dat[[i,10]]ScaleFacB68[[1]]}, ErrorBar[ ScaleFacB68[[1]]dat[[i,11]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Blue];
dat = sig3;
ps3 = ErrorListLogPlot[
Table[{{dat[[i,4]],dat[[i,10]]ScaleFacB68[[2]]}, ErrorBar[ScaleFacB68[[2]] dat[[i,11]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Darker[Green]];
dat = sig4;
ps4 = ErrorListLogPlot[
Table[{{dat[[i,4]],dat[[i,10]]ScaleFacB68[[3]]}, ErrorBar[ScaleFacB68[[3]]dat[[i,11]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Red];
dat = sig5;
ps5 = ErrorListLogPlot[
Table[{{dat[[i,4]],dat[[i,10]]ScaleFacB68[[4]]},ErrorBar[ScaleFacB68[[4]] dat[[i,11]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Black];
dat = sig58;
ps6 = ErrorListLogPlot[
Table[{{dat[[i,4]],dat[[i,10]]ScaleFacB68[[5]]}, ErrorBar[ScaleFacB68[[5]] dat[[i,11]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Orange];

ps = {ps2,ps3,ps4,ps5,ps6};

{ps}
];

Bolon67:=Module[{sig2,sig25,sig3,sig35,sig4,sig5,sig45,dat,
ps2,ps25,ps3,ps35,ps4,ps45,ps5,ps},
(* load data *)
sig2   = Import["Bolon67_Sig_PhotPi0_2.txt","Table"];
sig25 = Import["Bolon67_Sig_PhotPi0_2.5.txt","Table"];
sig3   = Import["Bolon67_Sig_PhotPi0_3.txt","Table"];
sig35 = Import["Bolon67_Sig_PhotPi0_3.5.txt","Table"];
sig4   = Import["Bolon67_Sig_PhotPi0_4.txt","Table"];
sig45 = Import["Bolon67_Sig_PhotPi0_4.5.txt","Table"];
sig5   = Import["Bolon67_Sig_PhotPi0_5.txt","Table"];
ScaleFacB67 =(1/10) {10^7,10^6,10^5,10^4,10^3,10^2,10^1};

dat = sig2;
ps2 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]10^3ScaleFacB67[[1]]}, ErrorBar[dat[[i,1]]-dat[[i,2]], 10^3ScaleFacB67[[1]]dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Blue];

dat = sig25;
ps25 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]10^3ScaleFacB67[[2]]}, ErrorBar[dat[[i,1]]-dat[[i,2]],10^3 ScaleFacB67[[2]]dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Darker[Green]];

dat = sig3;
ps3 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]10^3ScaleFacB67[[3]]}, ErrorBar[dat[[i,1]]-dat[[i,2]],10^3ScaleFacB67[[3]] dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Red];

dat = sig35;
ps35 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]10^3ScaleFacB67[[4]]}, ErrorBar[dat[[i,1]]-dat[[i,2]],10^3 ScaleFacB67[[4]]dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Black];

dat = sig4;
ps4 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]10^3ScaleFacB67[[5]]}, ErrorBar[dat[[i,1]]-dat[[i,2]],10^3ScaleFacB67[[5]]dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Orange];

dat = sig45;
ps45 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]10^3ScaleFacB67[[6]]}, ErrorBar[dat[[i,1]]-dat[[i,2]],10^3ScaleFacB67[[6]]dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Cyan];

dat = sig5;
ps5 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]10^3ScaleFacB67[[7]]},ErrorBar[dat[[i,1]]-dat[[i,2]],10^3ScaleFacB67[[7]] dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle-> Magenta];

ps = {ps2,ps25,ps3,ps35,ps4,ps45,ps5};

{ps}
];

(* Gomez 71 and Barish 74 *)
(* Barish 74 is included in Gomez 71*)
Gomez71:=Module[{sig6g,sig12g,sig18g,dat,ps6g,ps12g,ps18g,ps},
(* load data *)
sig6g   = Import["Gomez71_Sig_PhotPi0_6.txt","Table"];
sig12g = Import["Gomez71_Sig_PhotPi0_12.txt","Table"];
sig18g   = Import["Gomez71_Sig_PhotPi0_18.txt","Table"];
ScaleFacG71 ={10^0,10^0,10^0};

dat = sig6g;
ps6g = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]10^3ScaleFacG71[[1]]}, ErrorBar[ 10^3ScaleFacG71[[1]]dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Blue];

dat = sig12g;
ps12g = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]10^3ScaleFacG71[[2]]}, ErrorBar[ 10^3ScaleFacG71[[2]]dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Darker[Green]];

dat = sig18g;
ps18g = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]10^3ScaleFacG71[[3]]}, ErrorBar[10^3 ScaleFacG71[[3]]dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Red];

ps = {ps6g,ps12g,ps18g};

{ps}
];


Anderson:=
Module[{
Sig6,Sig9,Sig12, Sig15,BeamA4,BeamA6, BeamA10,
ps6,ps9,ps12,ps15,pb4,pb6,pb10,ps,pb,
dat,datasig, errsig,databa,errba},
(* load data *)
Sig6 = Import["Anderson_SigPhotPi0_6.txt","Table"];
Sig9 = Import["Anderson_SigPhotPi0_9.txt","Table"];
Sig12 = Import["Anderson_SigPhotPi0_12.txt","Table"];
Sig15 = Import["Anderson_SigPhotPi0_15.txt","Table"];
ScaleFac = {10^6,10^4,10^2,10^0};
ScaleFac = {1,1,1,1};
BeamA4 = Import["Anderson_BeamAPhotPi0_4.txt","Table"];
BeamA6 = Import["Anderson_BeamAPhotPi0_6.txt","Table"];
BeamA10 = Import["Anderson_BeamAPhotPi0_10.txt","Table"];


(* Differential cross section *)
(* Create Plots and the data object for the fit*)
(* There is no error on x *)
dat = Sig6;
ps6 = ErrorListLogPlot[
Table[{{-dat[[i,1]],dat[[i,4]]ScaleFac[[1]]}, ErrorBar[ScaleFac[[1]]dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Blue];

dat = Sig9;
ps9 = ErrorListLogPlot[
Table[{{-dat[[i,1]],dat[[i,4]]ScaleFac[[2]]},ErrorBar[ ScaleFac[[2]]dat[[i,5]] ]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Darker[Green]];

dat = Sig12;
ps12 = ErrorListLogPlot[
Table[{{-dat[[i,1]],dat[[i,4]]ScaleFac[[3]]}, ErrorBar[ScaleFac[[3]]dat[[i,5]] ]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Red];

dat = Sig15;
ps15 = ErrorListLogPlot[
Table[{{-dat[[i,1]],dat[[i,4]]ScaleFac[[4]]}, ErrorBar[ScaleFac[[4]]dat[[i,5]] ]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Black];

ps = {ps6,ps9,ps12,ps15};

(* Beam Asymmetry *)
(* Create Plots and the data object for the fit *)
dat = BeamA4;
pb4 = ErrorListPlot[
Table[{{-dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Blue];
dat = BeamA6;
pb6 = ErrorListPlot[
Table[{{-dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]] ]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Darker[Green]];
dat = BeamA10;
pb10 = ErrorListPlot[
Table[{{-dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]] ]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Red];

pb = {pb4,pb6,pb10};
{ps,pb}
];

Ratio := Module[{Rat4,Rat47,Rat82,pr4,pr47,pr82,dat,plotrat},
Rat4 = Import["Braun_RatioNP_PhotPi0_4.txt","Table"];
Rat47 = Import["Osborn2_RatioNP_PhotPi0_47.txt","Table"];
Rat82 = Import["Osborn2_RatioNP_PhotPi0_82.txt","Table"];

dat = Rat4;
pr4 = ErrorListPlot[
Table[{{-dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Blue];

dat = Rat47;
pr47 = ErrorListPlot[
Table[{{-dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Darker[Green]];

dat = Rat82;
pr82 = ErrorListPlot[
Table[{{-dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->Red];

plotrat ={pr4,pr47,pr82};

{Rat4,Rat47,Rat82,plotrat}]


(* Plot for Anderson data *)
pAfitS[name_] := LogPlot[{name[6,-t],name[9,-t],name[12,-t],name[15,-t]},{t,0.05,1.5},
PlotStyle-> {{Blue,Thick},{Thick,Darker[Green]},{Thick,Red},{Thick,Black}}, PlotPoints->2];
pAfitD[name_] := LogPlot[{name[6,-t],name[9,-t],name[12,-t],name[15,-t]},{t,0.05,1.5},
PlotStyle-> {{Blue,Thick,Dashed},{Thick,Darker[Green],Dashed},{Thick,Red,Dashed},{Thick,Black,Dashed}}, PlotPoints->2];

(* Plot for Braun68 data *)
pBfitS[name_]:= LogPlot[{
ScaleFacB68[[1]]name[2,-t],
ScaleFacB68[[2]]name[3,-t],
ScaleFacB68[[3]]name[4,-t],
ScaleFacB68[[4]]name[5,-t],
ScaleFacB68[[5]]name[5.8,-t]},{t,0.01,2},
PlotStyle-> {{Blue,Thick},{Thick,Darker[Green]},{Thick,Red},{Thick,Black},{Thick, Orange}}, PlotPoints->5];


SetDirectory[foldername];
ander = Anderson;
rat = Ratio;
gomez = Gomez71;
braun = Braun68;
bolon = Bolon67;

Amplitudes

(* Compute the differential cross section in μ barn/Gev^2 for pseudoscalar-photoproduction processes*)
DiffCrossSectPhot[plab_,t_,flag1_,flag2_,flag3_,a1ω_,a4ω_,a2B_,a1ρ_,a4ρ_,a1c_,a4c_,α0ω_,αpω_,α0c_,αpc_,α0B_,αpB_]:=
Module[{DCSp,DCSn,BA,res1,res2, res3,Fi,S,kt,pt,Σ,cost,sint,ϕ,u,αρ,αω,αB,αωc,ξω,ξB,ξρ,ξωc,f01p,f11p,f01m,f11m,F1,F2,F3,F4},
Σ = 2 MP^2+MPI^2;
S = MP^2+2MP Sqrt[plab^2];
kt = (t-MPI^2)/2/Sqrt[t];
pt = Sqrt[t-4MP^2]/2;
cost = (2S+t-Σ)/(4 kt pt);
sint = Sqrt[1-cost^2];

αρ = α0ω+αpω t;αω =α0ω+αpω t; αB  =α0B+αpB t; αωc =α0c+αpc t;

ξρ =Gamma[1-αρ](1-Exp[-I π αρ])/2 S^(αρ-1);
ξω =Gamma[1-αω](1-Exp[-I π αω])/2 S^(αω-1);
ξB =Gamma[1-αB] (1-Exp[-I π αB])/2 S^(αB-1);
ξωc =Gamma[1-αωc](1-Exp[-I π αωc])/2 S^(αωc-1)/Log[S];

F1 =-(t a1ω - 2MP a4ω)ξω -(t a1c- 2MP a4c )ξωc  - (t a1ρ  - 2MP a4ρ )ξρ ;
F2 = t a2B ξB;
F3 =t(2MP  a1ω - a4ω  )ξω + t(2MP  a1c  - a4c )ξωc + t(2MP  a1ρ  - a4ρ )ξρ;
F4 = 0;

f01p = kt Sqrt[t]/2. F1;
f11p = kt/2 F3;
f01m =  pt kt F2;
f11m = -pt kt Sqrt[t]F4;

res1 = 2(Abs[f01p]^2  + Abs[f01m]^2) Abs[sint]^2;
res1 = res1  + Abs[f11p-f11m]^2 (1-cost)^2 +  Abs[f11p+f11m]^2 (1+cost)^2;
res2 = 2(Abs[f01p]^2  - Abs[f01m]^2) Abs[sint]^2;
res2 = res2  + Abs[f11p]^2 (1-cost)^2 +  Abs[f11p]^2 (1+cost)^2;
res2 = res2  - Abs[f11m]^2 (1-cost)^2 -  Abs[f11m]^2 (1+cost)^2;

DCSp=Re[res1]/(64π MP^2plab^2)389.3;

(* for neutron target ; I=1 changes sign *)
F1 =-(t a1ω  - 2MP a4ω )ξω -(t a1c - 2MP a4c )ξωc  + (t a1ρ  - 2MP a4ρ )ξρ ;
F2 =- t a2B ξB;
F3 =t(2MP  a1ω- a4ω  )ξω + t(2MP  a1c  - a4c  )ξωc - t(2MP  a1ρ  - a4ρ )ξρ;
f01p = kt Sqrt[t]/2. F1;
f11p = kt/2 F3;
f01m =  pt kt F2;
f11m = -pt kt Sqrt[t]F4;

res3 = 2(Abs[f01p]^2  + Abs[f01m]^2) Abs[sint]^2;
res3 = res3  + Abs[f11p-f11m]^2 (1-cost)^2 +  Abs[f11p+f11m]^2 (1+cost)^2;

BA=Re[res2/res1];
DCSn=Re[res3/res1];
flag1 DCSp + flag2 BA + DCSn flag3];


(* Compute the differential cross section in μ barn/Gev^2 for pseudoscalar-photoproduction processes 
IN THE LEADING S APPROXIMATION *)
HELObserv[plab_,t_,flag1_,flag2_,flag3_,a1ω_,a4ω_,a2B_,a1ρ_,a4ρ_,a1c_,a4c_,α0ω_,αpω_,α0c_,αpc_,α0B_,αpB_]:=
Module[{DCSp,DCSn,BA,res1,res2, res3,S,
αρ,αω,αB,αωc,
ξω,ξB,ξρ,ξωc,
F1,F2,F3,F4},
S = MP^2+2MP Sqrt[plab^2];

αρ = α0ω+αpω t;αω =α0ω+αpω t; αB  =α0B+αpB t; αωc =α0c+αpc t;

ξρ =Gamma[1-αρ](1-Exp[-I π αρ])/2.0 S^(αρ-1);
ξω =Gamma[1-αω](1-Exp[-I π αω])/2.0 S^(αω-1);
ξB =Gamma[1-αB] (1-Exp[-I π αB])/2.0 S^(αB-1);
ξωc =Gamma[1-αωc](1-Exp[-I π αωc])/2.0 S^(αωc-1)/Log[S];

F1 =-(t a1ω - 2MP a4ω)ξω -(t a1c- 2MP a4c )ξωc  - (t a1ρ  - 2MP a4ρ )ξρ ;
F2 = t a2B ξB;
F3 =t(2MP  a1ω - a4ω  )ξω + t(2MP  a1c  - a4c )ξωc + t(2MP  a1ρ  - a4ρ )ξρ;
F4 = 0;

res1 = (Abs[F3]^2  -t Abs[F1]^2)/(4MP^2-t) + (Abs[F2]^2  -t Abs[F4]^2);
res2 = (Abs[F3]^2  -t Abs[F1]^2)/(4MP^2-t) - (Abs[F2]^2  -t Abs[F4]^2);

(* for neutron target ; I=1 changes sign *)
F1 =-(t a1ω  - 2MP a4ω )ξω -(t a1c - 2MP a4c )ξωc  + (t a1ρ  - 2MP a4ρ )ξρ ;
F2 =- t a2B ξB;
F3 =t(2MP  a1ω- a4ω  )ξω + t(2MP  a1c  - a4c  )ξωc - t(2MP  a1ρ  - a4ρ )ξρ;
res3 = (Abs[F3]^2  -t Abs[F1]^2)/(4MP^2-t) + (Abs[F2]^2  -t Abs[F4]^2);

DCSp=Re[res1]/(32π)389.3;
BA=Re[res2/res1];
DCSn=Re[res3/res1];
flag1 DCSp + flag2 BA + DCSn flag3];

Mass = {0.,MP,MPI,MP};
Cos2T[plab_,cos_,mass_]:=Block[
{ma,mb,mc,md,qi,qf,s,t,coef,tmin,tmax,som},
ma = mass[[1]];mb = mass[[2]]; mc = mass[[3]]; md = mass[[4]];
s=ma^2+mb^2+2mb Sqrt[plab^2+ma^2];
qi = Sqrt[(s-(ma+mb)^2)(s-(ma-mb)^2)/(4.s)]; 
qf = Sqrt[(s-(mc+md)^2)(s-(mc-md)^2)/(4.s)]; 
som = ma^2+mb^2+mc^2+md^2;
t = 2 qi qf cos+ som/2 - s/2. - (ma^2-mb^2)(mc^2-md^2)/(2. s);
coef = π/(qi qf);
tmin = (ma^2-mb^2-mc^2+md^2)^2/4./s-(qi-qf)^2;
tmax = (ma^2-mb^2-mc^2+md^2)^2/4./s-(qi+qf)^2; 
{t,coef, tmin, tmax}
]

T2Cos[plab_,t_,mass_]:=Block[
{ma,mb,mc,md,qi,qf,s,cos,coef,som},
ma = mass[[1]];mb = mass[[2]]; mc = mass[[3]]; md = mass[[4]];
s=ma^2+mb^2+2mb Sqrt[plab^2+ma^2];
qi = Sqrt[(s-(ma+mb)^2)(s-(ma-mb)^2)/(4.s)]; 
qf = Sqrt[(s-(mc+md)^2)(s-(mc-md)^2)/(4.s)]; 
som = ma^2+mb^2+mc^2+md^2;
cos = (t + (ma^2-mb^2)(mc^2-md^2)/(2. s)-som/2 + s/2 ) /(2 qi qf ); 
coef = π/(qi qf);
{cos,coef}
]

Fit with Eγ ≥ 4 GeV and t ≤ 3.0

(* select data *)
dat1 = dat2fit ;
pdat05 =  Select[dat1,0.01<Abs[Part[#,2]]≤0.5&];(* consider 0.01 < |t| < 0.5 *)
pdat10 =  Select[dat1,0.01<Abs[Part[#,2]]≤1.0&];(* consider 0.01 < |t| < 1.0 *)
pdat15 =  Select[dat1,0.01<Abs[Part[#,2]]≤1.5&];(* consider 0.01 < |t| < 1.5 *)
pdat20 =  Select[dat1,0.01<Abs[Part[#,2]]≤2.0&];(* consider 0.01 < |t| < 2.0 *)
pdat25 =  Select[dat1,0.01<Abs[Part[#,2]]≤2.5&];(* consider 0.01 < |t| < 2.5 *)
pdat30 =  Select[dat1,0.01<Abs[Part[#,2]]≤3.0&];(* consider 0.01 < |t| < 3.0 *)

(* pole, cut and trajectories : 8 parameters *)
FitFullCutα[datsig_,errsig_]:=NonlinearModelFit[datsig,
{DiffCrossSectPhot[pl,t,1,0,0,0,a4ω,0,a1ρ,0,a1ρc,a4ωc,α0v,αpv,α0c,αpc,0,0]},
{{a4ω,-20 },{a1ρ,-1},{a1ρc,0},{a4ωc,0},{α0v,0.5},{αpv,0.9},{α0c,0.5},{αpc,0.18}},
{pl,t}, Weights->1/errsig^2,VarianceEstimatorFunction->(1&),AccuracyGoal->5];

(* consider Eγ ≥ 4 and t ≤ 3.0 *)
xpdat = pdat30;
pdat =  Select[xpdat,Part[#,1]≥4.0&]; 
datsig = pdat[[All,1;;3]]; errsig = pdat[[All,4]];
(* with Full expression *)
nlmγπ=FitFullCutα[datsig,errsig]; χ2 =nlmγπ["EstimatedVariance",VarianceEstimatorFunction->Automatic];
Print["Eγ ≥ 4 GeV ; χ^2 = ",χ2]
plotSigA = Show[{ander[[1]],pAfitS[nlmγπ]},PlotRange-> {{0,1.5},{-5,0.7}}, AxesOrigin-> {0,-5}];
plotSigB= Show[{pBfitS[nlmγπ],braun[[1]]}, PlotRange->  {{0,2},{-4,15}}, AxesOrigin-> {0,-4},Frame->True,FrameStyle->FontSize->15];
GraphicsRow[{nlmγπ["ParameterTable"],plotSigA,plotSigB}]
(* The dimensions of the parameters are 
a4 GeV-3
a2 GeV-4
a1 GeV-2
αp GeV-2 
*)

Plot differential cross section with data

(* Anderson *)
leg = Graphics[Text[Column[{
Row[{Style["•",Blue,20],Style["  6 GeV",20]},Spacer[10]] ,
Row[{Style["•",Darker[Green],20],Style["  9 GeV",20]},Spacer[10]],
Row[{Style["•",Red,20],Style["12 GeV",20]},Spacer[10]],
Row[{Style["•",Black,20],Style["15 GeV",20]},Spacer[10]]
}],{1,-0.6}]];

psig1 = Show[{ander[[1]],plotSigA,leg},PlotRange-> {{0,1.5},{-5,0.7}}, AxesOrigin-> {0,-5},FrameLabel->{Style["-t ((GeV^2))",20],Style[ "dσ/dt (μb.(GeV^(-2)))",20]},Frame->True,FrameStyle->FontSize->15];

(* Braun *)
SetDirectory[foldername];
braun = Braun68;
leg = Graphics[Text[Column[{
Row[{Style["•",Blue,15],Style["   2 GeV x (10^4)",15]},Spacer[10]] ,
Row[{Style["•",Darker[Green],15],Style["   3 GeV x (10^3)",15]},Spacer[10]] ,
Row[{Style["•",Red,15],Style["   4 GeV x (10^2)",15]},Spacer[10]],
Row[{Style["•",Black,15],Style["   5 GeV x (10^1)",15]},Spacer[10]],
Row[{Style["•",Orange,15],Style["5.8 GeV x (10^0)",15]},Spacer[10]]
}],{1.5,12.5}]];

psig2 = Show[{plotSigB ,leg}, PlotRange->  {{0,2},{-4,15}}, AxesOrigin-> {0,-4}, FrameLabel->{Style["-t ((GeV^2))",20],Style[ "dσ/dt (μb.(GeV^(-2)))",20]},Frame->True,FrameStyle->FontSize->15];


(* Bolon67 *)
SetDirectory[foldername];
bolon = Bolon67;
pfit6 = LogPlot[{
ScaleFacB67[[1]]nlmγπ[2.0,-t],
ScaleFacB67[[2]]nlmγπ[2.5,-t],
ScaleFacB67[[3]]nlmγπ[3.0,-t],
ScaleFacB67[[4]]nlmγπ[3.5,-t],
ScaleFacB67[[5]]nlmγπ[4.0,-t],
ScaleFacB67[[6]]nlmγπ[4.5,-t],
ScaleFacB67[[7]]nlmγπ[5.0,-t]},{t,0.01,5},
PlotStyle-> {{Blue,Thick},{Thick,Darker[Green]},{Thick,Red},{Thick,Black},{Thick, Orange},{Thick, Cyan},{Thick, Magenta}}, PlotPoints->10];

leg = Graphics[Text[Column[{
Row[{Style["•",Blue,15],Style["2.0 GeV x (10^6)",15]},Spacer[10]] ,
Row[{Style["•",Darker[Green],15],Style["2.5 GeV x (10^5)",15]},Spacer[10]] ,
Row[{Style["•",Red,15],Style["3.0 GeV x (10^4)",15]},Spacer[10]],
Row[{Style["•",Black,15],Style["3.5 GeV x (10^3)",15]},Spacer[10]],
Row[{Style["•",Orange,15],Style["4.0 GeV x (10^02)",15]},Spacer[10]],
Row[{Style["•",Cyan,15],Style["4.5 GeV x (10^1)",15]},Spacer[10]],
Row[{Style["•",Magenta,15],Style["5.0 GeV x (10^0)",15]},Spacer[10]]
}],{0.5,-9}]];

psig3 = Show[{pfit6,bolon[[1]],leg}, PlotRange->  {{0,4.5},{-15,18}},AxesOrigin-> {0,-15}, FrameLabel->{Style["-t ((GeV^2))",20],Style[ "dσ/dt (μb.(GeV^(-2)))",20]},Frame->True,FrameStyle->FontSize->15];

(* Gomez 71 *)
SetDirectory[foldername];
gom71 = Gomez71;
pfit7 = LogPlot[{
ScaleFacG71[[1]]nlmγπ[6,-t],
ScaleFacG71[[2]]nlmγπ[12,-t],
ScaleFacG71[[3]]nlmγπ[18,-t]},{t,0.01,3.1},
PlotStyle-> {{Blue,Thick},{Thick,Darker[Green]},{Thick,Red}}];

leg = Graphics[Text[Column[{
Row[{Style["•",Blue,20],Style["  6 GeV",20]},Spacer[10]] ,
Row[{Style["•",Darker[Green],20],Style["12 GeV",20]},Spacer[10]] ,
Row[{Style["•",Red,20],Style["18 GeV",20]},Spacer[10]]
}],{.4,-6}]];

psig4=Show[pfit7,gom71[[1]],leg, FrameLabel->{Style["-t ((GeV^2))",20],Style[ "dσ/dt (μb.(GeV^(-2)))",20]}, Frame->True,FrameStyle->FontSize->15];

plotall=GraphicsGrid[{{psig1,psig2},{psig3,psig4}}]

Beam Asymmetry

fit

fitModelBA[plab_,t_,a2B_,α0B_,αpB_]:=DiffCrossSectPhot[plab,t,0,1,0,0,a4ω,a2B,a1ρ,0,a1ρc,a4ωc,α0v,αpv,α0c,αpc,α0B,αpB]/.nlmγπ["BestFitParameters"]

datBA = dataBA[[All,1;;3]];
errBA = dataBA[[All,4]];

(* axial pole and trajectories : 3 parameters *)
nlmBA:=NonlinearModelFit[datBA,{fitModelBA[pl,t,a2B,α0B,αpB]},
{{a2B,-1 },{α0B,0},{αpB,0.7}}, {pl,t}, Weights->1/errBA^2,VarianceEstimatorFunction->(1&),AccuracyGoal->5];

χ2 =nlmBA["EstimatedVariance",VarianceEstimatorFunction->Automatic];
Print["χ^2 = ",χ2]
nlmBA["ParameterTable"]
(* The dimensions of the parameters are 
a4 GeV-3
a2 GeV-4
a1 GeV-2
αp GeV-2 
*)

plot

modelBA[plab_,t_]:=DiffCrossSectPhot[plab,t,0,1,0,0,a4ω,-8.887,a1ρ,0,a1ρc,a4ωc,α0v,αpv,α0c,αpc,-0.193,1.021]/.nlmγπ["BestFitParameters"]
modelRat[plab_,t_]:=DiffCrossSectPhot[plab,t,0,0,1,0,a4ω,-8.887,a1ρ,0,a1ρc,a4ωc,α0v,αpv,α0c,αpc,-0.193,1.021]/.nlmγπ["BestFitParameters"]

pfit2 = Plot[{modelBA[4.,-t],modelBA[6.,-t],modelBA[10.,-t]},{t,0.1,1.5},
PlotStyle-> {{Blue,Thick},{Thick,Darker[Green]},{Thick,Red}}, PlotPoints->2];

leg = Graphics[Text[Column[{
Row[{Style["•",Blue,20],Style["  4 GeV",20]},Spacer[10]] ,
Row[{Style["•",Darker[Green],20],Style["  6 GeV",20]},Spacer[10]] ,
Row[{Style["•",Red,20],Style["10 GeV",20]},Spacer[10]]
}],{1.,0.5}]];

plotBA = Show[{ander[[2]],pfit2,leg},PlotRange-> {{0,1.5},{0,1.2}}, AxesOrigin-> {0,0},FrameLabel->{Style["-t ((GeV^2))",20],Style[ "Σ",20]},Frame->True,FrameStyle->FontSize->15];

(* Ratio neutron/proton *)
leg = Graphics[Text[Column[{
Row[{Style["•",Blue,20],Style["4.0 GeV",20]},Spacer[10]] ,
Row[{Style["•",Darker[Green],20],Style["4.7 GeV",20]},Spacer[10]] ,
Row[{Style["•",Red,20],Style["8.2 GeV",20]},Spacer[10]]
}],{0.3,1.5}]];

pfit3 = Plot[{modelRat[4,-t],modelRat[4.7,-t],modelRat[8.2,-t]},{t,0,1.5},
PlotStyle-> {{Blue,Thick},{Thick,Darker[Green]},{Thick,Red}}, PlotPoints->2];

plotRat = Show[{rat[[4]],pfit3,leg}, PlotRange->  {{0,1.5},{0,2}}, AxesOrigin-> {0,0},
FrameLabel->{Style["-t ((GeV^2))",20],Style[ "dσ(n)/dσ(p)",20]},Frame->True,FrameStyle->FontSize->15];

GraphicsRow[{plotBA,plotRat}]

Target and Recoil Asymmetries

SetDirectory[foldername];

tarA1 = Import["Beinlein_TargetA_PhotPi0_4.txt","Table"];
tarA2 = Import["Booth_TargetA_PhotPi0_4.txt","Table"];
recA = Import["Deutsch72_RecoilA4-6.txt","Table"];

dat = tarA1;
ps1 = ErrorListPlot[
Table[{{-dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->{●,15},PlotStyle->Blue];
dat = tarA2;
ps2 = ErrorListPlot[
Table[{{-dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->{■,15},PlotStyle->Red];
dat = recA;
ps3 = ErrorListPlot[
Table[{{-dat[[i,5]],dat[[i,7]]}, ErrorBar[dat[[i,8]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->{●,15},PlotStyle->Black];


(* Target and recoil Asymmetries are the same since F4=0 *)
targetA[plab_,t_,a1ω_,a4ω_,a2B_,a1ρ_,a4ρ_,a1c_,a4c_,α0ω_,αpω_,α0c_,αpc_,α0B_,αpB_]:=
Module[{res1,S,
αρ,αω,αB,αωc,
ξω,ξB,ξρ,ξωc,
F1,F2,F3,F4},
S = MP^2+2MP Sqrt[plab^2];

(*αρ = 0.5+0.8t;αω = 0.5 +0.8t; αB  = 0 + 0.7t; αωc = 0.5+0.15t;*)
αρ = α0ω+αpω t;αω =α0ω+αpω t; αB  =α0B+αpB t; αωc =α0c+αpc t;

ξρ =-Gamma[1-αρ](1-Exp[-I π αρ])/2 S^(αρ-1);
ξω =-Gamma[1-αω](1-Exp[-I π αω])/2 S^(αω-1);
ξB =-Gamma[1-αB] (1-Exp[-I π αB])/2 S^(αB-1);
ξωc =Gamma[1-αωc](1-Exp[-I π αωc])/2 S^(αωc-1)/Log[S];

F1 =-(t a1ω - 2MP a4ω)ξω -(t a1c- 2MP a4c )ξωc  - (t a1ρ  - 2MP a4ρ )ξρ ;
F2 = t a2B ξB;
F3 =t(2MP  a1ω - a4ω  )ξω + t(2MP  a1c  - a4c )ξωc + t(2MP  a1ρ  - a4ρ )ξρ;
F4 = 0;

res1 = F1 Conjugate[F3]/(4MP^2-t) - F4 Conjugate[F2];
Im[res1]Sqrt[-t]/(16π)389.3];

tarA[pl_,t_]:=targetA[pl,t,0,a4ω,-8.887,a1ρ,0,a1ρc,a4ωc,α0v,αpv,α0c,αpc,-0.193,1.021]/.nlmγπ["BestFitParameters"];

pth1 =Plot[{tarA[4.,-t],tarA[5.,-t],tarA[6.,-t]},{t,0,2},PlotStyle->{{Thick, Blue},{Thick, Darker[Green]},{Thick, Red}}, PlotRange->{{0,2},{-1.1,0.3}},FrameLabel->{Style["-t ((GeV^2))",20],Style[ "R",20]},Frame->True,FrameStyle->FontSize->15];

pth2 =Plot[tarA[4.,-t],{t,0,2},PlotStyle->{Thick, Black}, PlotRange->{{0,2},{-1.1,0.3}},FrameLabel->{Style["-t ((GeV^2))",20],Style[ "T",20]},Frame->True,FrameStyle->FontSize->15];

leg = Graphics[Text[Column[{
Row[{Style["-",Red,40],Style["6 GeV",20]},Spacer[10]] ,
Row[{Style["-",Darker[Green],40],Style["5 GeV",20]},Spacer[10]],
Row[{Style["-",Blue,40],Style["4 GeV",20]},Spacer[10]]
}],{1.3,-0.65}]];
ptar = Show[pth1,ps3,leg];
prec = Show[pth2,ps1,ps2];

GraphicsGrid[{{Show[pth2,ps1,ps2],Show[pth1,ps3,leg]},{Show[pth2,pth1,ps1,ps2,ps3]}}]

HIgh Energy Limit vs Full Expression

FullExp[plab_,t_,fl1_,fl2_,fl3_]:=DiffCrossSectPhot[plab,t,fl1,fl2,fl3,0,a4ω,-8.887,a1ρ,0,a1ρc,a4ωc,α0v,αpv,α0c,αpc,-0.193,1.021]/.nlmγπ["BestFitParameters"]
Approx[plab_,t_,fl1_,fl2_,fl3_]:=HELObserv[plab,t,fl1,fl2,fl3,0,a4ω,-8.887,a1ρ,0,a1ρc,a4ωc,α0v,αpv,α0c,αpc,-0.193,1.021]/.nlmγπ["BestFitParameters"]


papprox = LogPlot[{Approx[2,-t,1,0,0],Approx[3,-t,1,0,0],Approx[4,-t,1,0,0],Approx[5,-t,1,0,0],Approx[6,-t,1,0,0],Approx[9,-t,1,0,0]},{t,0.05,1.5}, PlotStyle->{{Blue,Thick},{Thick,Darker[Green]},{Thick,Red},{Thick,Black},{Thick, Orange},{Thick, Magenta}}, PlotPoints->2];

pfull = LogPlot[{FullExp[2,-t,1,0,0],FullExp[3,-t,1,0,0],FullExp[4,-t,1,0,0],FullExp[5,-t,1,0,0],FullExp[6,-t,1,0,0],FullExp[9,-t,1,0,0]},{t,0.05,1.5}, PlotStyle->{{Blue,Thick,Dashed},{Thick,Darker[Green],Dashed},{Thick,Red,Dashed},{Thick,Black,Dashed},{Thick, Orange,Dashed},{Thick, Magenta,Dashed}}, PlotPoints->2];

psig = Show[papprox,pfull, PlotRange->{{0,1.5},All}, AxesOrigin->{0,-4}];

pfit2b = Plot[{FullExp[4,-t,0,1,0],FullExp[6,-t,0,1,0],FullExp[10,-t,0,1,0]},{t,0,1.5},
PlotStyle-> {{Blue,Thick,Dashed},{Thick,Darker[Green],Dashed},{Thick,Red,Dashed}}, PlotPoints->2];
pfit3b = Plot[{FullExp[4,-t,0,0,1],FullExp[4.7,-t,0,0,1],FullExp[8.2,-t,0,0,1]},{t,0,1.5},
PlotStyle-> {{Blue,Thick,Dashed},{Thick,Darker[Green],Dashed},{Thick,Red},Dashed}, PlotPoints->2];

pfit2c = Plot[{Approx[4,-t,0,1,0],Approx[6,-t,0,1,0],Approx[10,-t,0,1,0]},{t,0,1.5},
PlotStyle-> {{Blue,Thick},{Thick,Darker[Green]},{Thick,Red}}, PlotPoints->2];
pfit3c = Plot[{Approx[4,-t,0,0,1],Approx[4.7,-t,0,0,1],Approx[8.2,-t,0,0,1]},{t,0,1.5},
PlotStyle-> {{Blue,Thick},{Thick,Darker[Green]},{Thick,Red}}, PlotPoints->2];

pba = Show[pfit2c,pfit2b];
prat = Show[pfit3c,pfit3b];

GraphicsGrid[{{psig,pba,prat}}]

leg = Graphics[Text[Column[{
Row[{Style["-",Blue,20],Style["2 GeV",20]},Spacer[10]] ,
Row[{Style["-",Darker[Green],20],Style["3 GeV",20]},Spacer[10]],
Row[{Style["-",Red,20],Style["4 GeV",20]},Spacer[10]],
Row[{Style["-",Black,20],Style["5 GeV",20]},Spacer[10]],
Row[{Style["-",Orange,20],Style["6 GeV",20]},Spacer[10]],
Row[{Style["-",Magenta,20],Style["9 GeV",20]},Spacer[10]]
}],{.1,-2}]];
psig6 = Show[papprox,leg,pfull, PlotRange->{{0,1.5},All}, AxesOrigin->{0,-4},FrameLabel->{Style["-t ((GeV^2))",20],Style[ "dσ/dt (μb.(GeV^(-2)))",20]},Frame->True,FrameStyle->FontSize->15]

FullExp[5,-1,1,0,0]/Approx[5,-1,1,0,0]

papprox = Plot[{Approx[2,-t,1,0,0]/FullExp[2,-t,1,0,0],
Approx[3,-t,1,0,0]/FullExp[3,-t,1,0,0],
Approx[4,-t,1,0,0]/FullExp[4,-t,1,0,0],
Approx[5,-t,1,0,0]/FullExp[5,-t,1,0,0],
Approx[6,-t,1,0,0]/FullExp[6,-t,1,0,0],
Approx[9,-t,1,0,0]/FullExp[9,-t,1,0,0]},{t,0.05,1.5}, PlotStyle->{{Blue,Thick},{Thick,Darker[Green]},{Thick,Red},{Thick,Black},{Thick, Orange},{Thick, Magenta}}, PlotPoints->2,FrameLabel->{Style["-t ((GeV^2))",20],Style[ "dσ/dt (μb.(GeV^(-2)))",20]}];
leg = Graphics[Text[Column[{
Row[{Style["-",Blue,20],Style["2 GeV",20]},Spacer[10]] ,
Row[{Style["-",Darker[Green],20],Style["3 GeV",20]},Spacer[10]],
Row[{Style["-",Red,20],Style["4 GeV",20]},Spacer[10]],
Row[{Style["-",Black,20],Style["5 GeV",20]},Spacer[10]],
Row[{Style["-",Orange,20],Style["6 GeV",20]},Spacer[10]],
Row[{Style["-",Magenta,20],Style["9 GeV",20]},Spacer[10]]
}],{.2,1.5}]];
psig7 = Show[leg,papprox,Frame->True,FrameStyle->FontSize->15, PlotRange->{{0,1.5},All}, AxesOrigin->{0,1},FrameLabel->{Style["-t ((GeV^2))",20]},Frame->True,FrameStyle->FontSize->15]

Cut vs Daugther

(* Compute the differential cross section in μ barn/Gev^2 for pseudoscalar-photoproduction processes*)
SimpleModel[plab_,t_,gv_,gc_,gd_,α0v_,αpv_,α0c_,αpc_,α0d_,αpd_]:=
Module[{DCSp,res1,S,kt,pt,Σ,cost,sint,ϕ,u,αv,αc,αd,ξv,ξd,ξc,f01p,f11p,f01m,f11m,F1,F2,F3,F4},
Σ = 2 MP^2+MPI^2;
S = MP^2+2MP Sqrt[plab^2];
kt = (t-MPI^2)/2/Sqrt[t];
pt = Sqrt[t-4MP^2]/2;
cost = (2S+t-Σ)/(4 kt pt);
sint = Sqrt[1-cost^2];

αv = α0v+αpv t; αc =α0c+αpc t;αd =α0d+αpd t;

ξv =Gamma[1-αv](1-Exp[-I π αv])/2 S^(αv-1);
ξd =Gamma[1-αd](1-Exp[-I π αd])/2 S^(αd-1);
ξc =Gamma[1-αc](1-Exp[-I π αc])/2 S^(αc-1)/Log[S];

F1 =-( 0- 2MP gv )ξv -(0- 2MP gc )ξc  - ( 0 - 2MP gd )ξd ;
F2 = 0;
F3 =t( 0- gv  )ξv + t(0  - gc )ξc + t(0  - gd )ξd;
F4 = 0;

f01p = kt Sqrt[t]/2. F1;
f11p = kt/2 F3;
f01m =  pt kt F2;
f11m = -pt kt Sqrt[t]F4;

res1 = 2(Abs[f01p]^2  + Abs[f01m]^2) Abs[sint]^2;
res1 = res1  + Abs[f11p-f11m]^2 (1-cost)^2 +  Abs[f11p+f11m]^2 (1+cost)^2;

DCSp=Re[res1]/(64π MP^2plab^2)389.3;

 DCSp ];


Clear[fig]
pole[plab_,t_] := SimpleModel[plab,t,1,0,0,α0v,αpv,0,0,0,0]/.nlmγπ["BestFitParameters"];
polecut[plab_,t_] := SimpleModel[plab,t,1,.1,0,α0v,αpv,α0c,αpc,0,0]/.nlmγπ["BestFitParameters"];
poledau[plab_,t_] := SimpleModel[plab,t,1,0,0.5,α0v,αpv,0,0,α0v-1,αpv]/.nlmγπ["BestFitParameters"];
polecutdau[plab_,t_] := SimpleModel[plab,t,1,0.1,0.5,α0v,αpv,α0c,αpc,α0v-1,αpv]/.nlmγπ["BestFitParameters"];

leg = Graphics[Text[Column[{
Row[{Style["-",Blue,30],Style["Pole",20]},Spacer[10]] ,
Row[{Style["-",Darker[Green],30],Style["Pole+Cut",20]},Spacer[10]],
Row[{Style["-",Orange,30],Style["Pole+Daugther",20]},Spacer[10]],
Row[{Style["-",Red,30],Style["Pole+Cut+Daugther",20]},Spacer[10]]
}],{-1,-17}]];

En1 = 9; En2 = 3;
fig[1]=LogPlot[{pole[En1,t],polecut[En1,t],poledau[En1,t],polecutdau[En1,t]},{t,-3,-.01}, PlotPoints->3,
PlotStyle->{{Thick,Blue},{Thick,Darker[Green]},{Thick,Orange},{Thick,Red}}];
fig[2]=LogPlot[{100pole[En2,t],100polecut[En2,t],100poledau[En2,t],100polecutdau[En2,t]},{t,-3,-.01}, PlotPoints->3,PlotStyle->{{Thick,Blue},{Thick,Darker[Green]},{Thick,Orange},{Thick,Red}}];
Show[fig[1],fig[2],leg,PlotRange->All,Frame->True,FrameStyle->FontSize->15]

leg = Graphics[Text[Column[{
Row[{Style["-",Blue,40],Style["Pole+Cut",20]},Spacer[10]],
Row[{Style["--",Darker[Green],20],Style["Pole+Daugther",20]},Spacer[10]]
}],{1,-17}]];

fig[3]=LogPlot[{polecut[En1,-t],poledau[En1,-t]},{t,.01,3}, PlotPoints->3,
PlotStyle->{{Thick,Blue},{Thick,Darker[Green],Dashed},{Thick,Orange},{Thick,Red}}];
fig[4]=LogPlot[{polecut[En2,-t],poledau[En2,-t]},{t,.01,3}, PlotPoints->3,PlotStyle->{{Thick,Blue},{Thick,Darker[Green],Dashed},{Thick,Orange},{Thick,Red}}];

pcutdau=Show[fig[3],fig[4],leg,PlotRange->All,Frame->True,FrameStyle->FontSize->15,FrameLabel->{Style["-t ((GeV^2))",20],Style[ "dσ/dt (μb.(GeV^(-2)))",20]},
Epilog->{Text[Style["Eγ = 9 GeV",Black,20],{1.3,-8}],Text[Style["Eγ = 3 GeV",Black,20],{1.3,-5}]}]

Prediction for CLAS

tvec = Range[-3.0,-0.01,0.02];
Elab1 = 3;  (* Beam energy in GeV *)
Elab2 = 4;  (* Beam energy in GeV *)
Elab3 = 5;  (* Beam energy in GeV *)
Elab4 = 6;  (* Beam energy in GeV *)
Elab5 = 9;  (* Beam energy in GeV *)

fullmodel[plab_,t_]:=DiffCrossSectPhot[plab,t,1,0,0,0,a4ω,-8.887,a1ρ,0,a1ρc,a4ωc,α0v,αpv,α0c,αpc,-0.193,1.021]/.nlmγπ["BestFitParameters"];

ord1 = Table[fullmodel[Elab1,tvec[[it]]]/T2Cos[Elab1,tvec,Mass][[2]],{it,Length[tvec]}];
ord2 = Table[fullmodel[Elab2,tvec[[it]]]/T2Cos[Elab2,tvec,Mass][[2]],{it,Length[tvec]}];
ord3 = Table[fullmodel[Elab3,tvec[[it]]]/T2Cos[Elab3,tvec,Mass][[2]],{it,Length[tvec]}];
ord4 = Table[fullmodel[Elab4,tvec[[it]]]/T2Cos[Elab4,tvec,Mass][[2]],{it,Length[tvec]}];
ord5 = Table[fullmodel[Elab5,tvec[[it]]]/T2Cos[Elab4,tvec,Mass][[2]],{it,Length[tvec]}];

pcos1 = ListLogPlot[Transpose[{T2Cos[Elab1,tvec,Mass][[1]],ord1}],
PlotRange->{{-1,1},{10^-4,1}}, AxesOrigin->{-1,0}, Joined->True, PlotStyle->{Blue, Thick}, Frame->True];
pcos2 = ListLogPlot[Transpose[{T2Cos[Elab2,tvec,Mass][[1]],ord2}],
PlotRange->{{-1,1},{10^-4,1}}, AxesOrigin->{-1,0}, Joined->True, PlotStyle->{Darker[Green], Thick}, Frame->True];
pcos3 = ListLogPlot[Transpose[{T2Cos[Elab3,tvec,Mass][[1]],ord3}],
PlotRange->{{-1,1},{10^-4,1}}, AxesOrigin->{-1,0}, Joined->True, PlotStyle->{Red, Thick}, Frame->True];
pcos4 = ListLogPlot[Transpose[{T2Cos[Elab4,tvec,Mass][[1]],ord4}],
PlotRange->{{-1,1},{10^-4,1}}, AxesOrigin->{-1,0}, Joined->True, PlotStyle->{Black, Thick}, Frame->True];
pcos5 = ListLogPlot[Transpose[{T2Cos[Elab5,tvec,Mass][[1]],ord4}],
PlotRange->{{-1,1},{10^-4,1}}, AxesOrigin->{-1,0}, Joined->True, PlotStyle->{Orange, Thick}, Frame->True];

leg = Graphics[Text[Column[{
Row[{Style["-",Blue,40],Style["3 GeV",20]},Spacer[10]] ,
Row[{Style["-",Darker[Green],40],Style["4 GeV",20]},Spacer[10]],
Row[{Style["-",Red,40],Style["5 GeV",20]},Spacer[10]],
Row[{Style["-",Black,40],Style["6 GeV",20]},Spacer[10]],
Row[{Style["-",Orange,40],Style["9 GeV",20]},Spacer[10]]
}],{-0.7,-3}]];
pclas1=Show[pcos1,pcos2,pcos3,pcos4,pcos5,leg,FrameLabel->{Style["Cos θ",20],Style[ "dσ/dt (μb.(GeV^(-2)))",20]},FrameStyle->FontSize->15];


pfull = LogPlot[{fullmodel[Elab1,t],fullmodel[Elab2,t],fullmodel[Elab3,t],fullmodel[Elab4,t],fullmodel[Elab5,t]},{t,-3,-0.03}, PlotStyle->{{Blue,Thick},{Thick,Darker[Green]},{Thick,Red},{Thick,Black},{Thick, Orange}}, PlotPoints->5];

leg = Graphics[Text[Column[{
Row[{Style["-",Blue,40],Style["3 GeV",20]},Spacer[10]] ,
Row[{Style["-",Darker[Green],40],Style["4 GeV",20]},Spacer[10]],
Row[{Style["-",Red,40],Style["5 GeV",20]},Spacer[10]],
Row[{Style["-",Black,40],Style["6 GeV",20]},Spacer[10]],
Row[{Style["-",Orange,40],Style["9 GeV",20]},Spacer[10]]
}],{-3.5,-3}]];

pclas2 = Show[pfull,leg,PlotRange->{{-4,-0.08},All},AxesOrigin->{0,-12},Frame->True,FrameLabel->{Style["t ((GeV^2))",20],Style[ "dσ/dt (μb.(GeV^(-2)))",20]},FrameStyle->FontSize->15];

GraphicsRow[{pclas1,pclas2}]

(*Export["SigCLAS1.pdf",%,ImageSize->Scaled[1.5]];*)

Export Figures

(* SetDirectory["/Users/vincent/Dropbox/mes articles/Regge/Pi0Phot/latex"];*)
Export["BeamAsym.pdf",plotBA,ImageSize->Scaled[1.5]];
Export["Ratio.pdf",plotRat,ImageSize->Scaled[1.5]];
Export["Sig1.pdf",psig1,ImageSize->Scaled[1.5]];
Export["Sig2.pdf",psig2,ImageSize->Scaled[1.5]];
Export["Sig3.pdf",psig3,ImageSize->Scaled[1.5]];
Export["Sig4.pdf",psig4,ImageSize->Scaled[1.5]];
Export["Sig6.pdf",psig6,ImageSize->Scaled[1.5]];
Export["CLAS1.pdf",pclas1,ImageSize->Scaled[1.5]];
Export["CLAS2.pdf",pclas2,ImageSize->Scaled[1.5]];
Export["SigAll.pdf",plotall,ImageSize->Scaled[2.0]];
Export["Target.pdf",prec,ImageSize->Scaled[1.5]];
Export["Recoil.pdf",ptar,ImageSize->Scaled[1.5]];
Export["CutDau.pdf",pcutdau,ImageSize->Scaled[1.5]];