Exit[]

Neutral Pion Photoproduction

ANDERSON 1970 http://hepdata.cedar.ac.uk/view/ins67154

MP= 0.93827;MPI=0.1349766; Mass = {0.,MP,MPI,MP}; 

SetDirectory["/Users/vincent/Dropbox/mes articles/Regge/data"];
Needs["PlotLegends`"]
Needs["ErrorBarPlots`"]
Print["Need to install the package ErrorBarLogPlots from http://library.wolfram.com/infocenter/MathSource/6747/"]
Needs["ErrorBarLogPlots`"]

SetDirectory["/Users/vincent/Dropbox/mes articles/Regge/data"];
Print["Change the directory !"]
(* load data from: 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"];

Plot Data

Effective Trajectory

set = {1,2,3,4,5,6,7,8,10,11,12};
ratio6 =(1/2) Log[6^2Sig6[[All,4]]/12^2/Sig12[[set,4]]]/Log[(MP^2+2MP 6)/(MP^2+2MP 12)];
ratio9 =(1/2) Log[9^2Sig9[[All,4]]/12^2/Sig12[[All,4]]]/Log[(MP^2+2MP 9)/(MP^2+2MP 12)];
ratio15 =(1/2) Log[15^2Sig15[[All,4]]/12^2/Sig12[[2;;12,4]]]/Log[(MP^2+2MP 15)/(MP^2+2MP 12)];

err6 = (1/2)(Sig6[[All,5]]/Sig6[[All,4]] + Sig12[[set,5]]/Sig12[[set,4]]);
err9 = (1/2)(Sig9[[All,5]]/Sig9[[All,4]] + Sig12[[All,5]]/Sig12[[All,4]]);
err15 = (1/2)(Sig15[[All,5]]/Sig15[[All,4]] + Sig12[[2;;12,5]]/Sig12[[2;;12,4]]);

framestyle = {PlotMarkers->Automatic,Frame->True,FrameStyle->FontSize->20};
plot[6]=ErrorListPlot[Table[{{Sig6[[i,1]]+0.01,ratio6[[i]]},ErrorBar[err6[[i]]]},{i,11}],framestyle, PlotStyle->Blue];
plot[9]=ErrorListPlot[Table[{{Sig9[[i,1]]+0.00,ratio9[[i]]},ErrorBar[err9[[i]]]},{i,12}],framestyle, PlotStyle->Red];
plot[15]=ErrorListPlot[Table[{{Sig15[[i,1]]-0.01,ratio15[[i]]},ErrorBar[err15[[i]]]},{i,11}],framestyle, PlotStyle->Darker[Green]];

leg = Graphics[Text[Column[{
Row[{Style["•",Blue,20],Style["Ratio   6 GeV to 12 GeV",15]},Spacer[10]] ,
Row[{Style["•",Red,20],Style["Ratio   9 GeV to 12 GeV",15]},Spacer[10]],
Row[{Style["•",Darker[Green],20],Style["Ratio 15 GeV to 12 GeV",15]},Spacer[10]]
}],{-1.1,0.3}]];

plotAlpEff=Show[plot[6],plot[9],plot[15],leg,PlotRange->All,FrameLabel->{"t ((GeV^2))", "α(t)"}]

Plot Differential Cross Section

AndersonPlot:=Block[{ps6,ps9,ps12,ps15,dat,style},
style = {PlotMarkers->Automatic,Joined->True,FrameTicks->{{{0.01,0.05,0.1,0.5,1},None},{Range[-1.4,0,0.2],None}}};
dat = Sig6;
ps6 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]]]},{i,Length[dat]}],style,PlotStyle->{Thick,Blue}];
dat = Sig9;
ps9 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]},ErrorBar[ dat[[i,5]] ]},{i,Length[dat]}],style,PlotStyle->{Thick,Darker[Green]}];
dat = Sig12;
ps12 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]] ]},{i,Length[dat]}],style,PlotStyle->{Thick,Red}];
dat = Sig15;
ps15 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]] ]},{i,Length[dat]}],style,PlotStyle->{Thick,Black}];
(* return *)
{ps6,ps9,ps12,ps15}
];

ander = AndersonPlot;
framestyle = {PlotRange->{{-1.4,0},{-5,1}},Frame->True, FrameStyle->FontSize->25, AxesOrigin->{0,2},GridLines->Automatic};
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}]];
plotCS=Show[ander,leg,framestyle ,FrameLabel->{"t ((GeV^2))", "dσ/dt μb/(GeV^2)"}]

Fit Data

Prepare data for fit

dat2fit := Block[{dat,datfit1,datfit},
(* create a file with the entries 
    Eγ  t  dσ/dt  Δdσ/dt 
*)
datfit = {}; (* initialization *)
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 *)
];

(* without joined the points *)
AndersonPlot2:=Block[{ps6,ps9,ps12,ps15,dat,style},
style = {PlotMarkers->Automatic,FrameTicks->{{{0.01,0.05,0.1,0.5,1},None},{Range[-1.4,0,0.2],None}}};
dat = Sig6;
ps6 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]]]},{i,Length[dat]}],style,PlotStyle->{Thick,Blue}];
dat = Sig9;
ps9 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]},ErrorBar[ dat[[i,5]] ]},{i,Length[dat]}],style,PlotStyle->{Thick,Darker[Green]}];
dat = Sig12;
ps12 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]] ]},{i,Length[dat]}],style,PlotStyle->{Thick,Red}];
dat = Sig15;
ps15 = ErrorListLogPlot[
Table[{{dat[[i,1]],dat[[i,4]]}, ErrorBar[dat[[i,5]] ]},{i,Length[dat]}],style,PlotStyle->{Thick,Black}];
(* return *)
{ps6,ps9,ps12,ps15}
];

Model 1

(* Compute the differential cross section in μ barn/Gev^2 for pseudoscalar-photoproduction processes*)
DiffCrossSect1[plab_,t_,g1_,g4_,α0_,αp_]:=
Module[{S,kt,pt,Σ,sint, cost,App, Amm, Apm, Amp,Ai},
Σ = 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];

Ai =Model1[plab,t,g1,g4,α0,αp];

App = kt/Sqrt[2]sint (Sqrt[t] (Ai[[1]] + 2MP Ai[[4]]) -2 pt(Ai[[1]] + t Ai[[2]])  );
Amm = kt/Sqrt[2]sint (Sqrt[t] (Ai[[1]] + 2MP Ai[[4]]) +2 pt(Ai[[1]] + t Ai[[2]])  );
Apm = kt/Sqrt[2] (1-cost)(- (2MP Ai[[1]] -t Ai[[4]]) +2 pt(Ai[[1]] + Sqrt[t] Ai[[3]])  );
Amp = kt/Sqrt[2] (1+cost)(- (2MP Ai[[1]] -t Ai[[4]]) +2 pt(Ai[[1]] - Sqrt[t] Ai[[3]])  );

(*return*)
(Abs[App]^2 + Abs[Amm]^2 +Abs[Apm]^2 +Abs[Amp]^2 )/(64π MP^2plab^2)389.3
];

Model1[plab_,t_,g1_,g4_,α0_,αp_]:=
(* Return A1,A2,A3,A4 *)
Module[{S,kt,pt,Σ,cost,sint,ϕ,u,αv,ξv,A1,A2,A3,A4},
Σ = 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 = α0+αp t;(* trajectory *)
ξv =αv(αv+1)(1-Exp[-I π αv])/2/Sin[π αv] S^(αv-1); (* signature factor *)

A1 = g1 t ξv;
A2 = -g1 ξv;
A3 = 0;
A4  = g4 ξv;
{A1,A2,A3,A4}
];


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

Fit data with model 1

(* load data *)
dat = dat2fit ;datsig = dat[[All,1;;3]]; errsig = dat[[All,4]];
ander2 = AndersonPlot2; (* data without joining the points *)
(* 2 residues and one linear trajectory : 4 parameters *)
nlmγπ1=NonlinearModelFit[datsig,{DiffCrossSect1[pl,t,g1,g4,α0,αp]},
{{g4,20 },{g1,1},{α0,0.5},{αp,0.9}},{pl,t}, Weights->1/errsig^2,VarianceEstimatorFunction->(1&),AccuracyGoal->5];

 χ2 =nlmγπ1["EstimatedVariance",VarianceEstimatorFunction->Automatic];
Print["Resulting χ^2 = ",χ2]
nlmγπ1["ParameterTable"]
plotSigA = Show[{ander2,pAfit[nlmγπ1]},leg, FrameLabel->{"t ((GeV^2))", "dσ/dt μb/(GeV^2)"},framestyle]

(* The dimensions of the parameters are 
a4 GeV-3
a2 GeV-4
a1 GeV-2
αp GeV-2 
*)

Model 2

(* Compute the differential cross section in μ barn/Gev^2 for pseudoscalar-photoproduction processes*)
DiffCrossSect2[plab_,t_,g1_,g4_,α0_,αp_, g1c_,g4c_,α0c_,αpc_]:=
Module[{S,kt,pt,Σ,sint, cost,App, Amm, Apm, Amp,Ai},
Σ = 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];

Ai =Model2[plab,t,g1,g4,α0,αp, g1c,g4c,α0c,αpc];

App = kt/Sqrt[2]sint (Sqrt[t] (Ai[[1]] + 2MP Ai[[4]]) -2 pt(Ai[[1]] + t Ai[[2]])  );
Amm = kt/Sqrt[2]sint (Sqrt[t] (Ai[[1]] + 2MP Ai[[4]]) +2 pt(Ai[[1]] + t Ai[[2]])  );
Apm = kt/Sqrt[2] (1-cost)(- (2MP Ai[[1]] -t Ai[[4]]) +2 pt(Ai[[1]] + Sqrt[t] Ai[[3]])  );
Amp = kt/Sqrt[2] (1+cost)(- (2MP Ai[[1]] -t Ai[[4]]) +2 pt(Ai[[1]] - Sqrt[t] Ai[[3]])  );

(*return*)
(Abs[App]^2 + Abs[Amm]^2 +Abs[Apm]^2 +Abs[Amp]^2 )/(64π MP^2plab^2)389.3
];

Model2[plab_,t_,g1_,g4_,α0_,αp_, g1c_,g4c_,α0c_,αpc_]:=
(* Return A1,A2,A3,A4 *)
Module[{S,kt,pt,Σ,cost,sint,ϕ,u,αv,ξv,αc,ξc,A1,A2,A3,A4},
Σ = 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 = α0+αp t;(* vector trajectory *)
αc = α0c+αpc t;(* vector cut trajectory *)
(*
ξv =Gamma[1-αv](1-Exp[-I π αv])/2 S^(αv-1); (* signature factor *)
ξc =Gamma[1-αc](1-Exp[-I π αc])/2 S^(αc-1)/Log[S]; (* signature factor *)
*)

ξv =αv(αv+1)(1-Exp[-I π αv])/2/Sin[π αv] S^(αv-1); (* signature factor *)
ξc =αc(αc+1)(1-Exp[-I π αc])/2/Sin[π αc] S^(αc-1)/Log[S]; (* signature factor *)

A1 = g1 t ξv + g1c t ξc;
A2 = -g1 ξv-g1c ξc;
A3 = 0;
A4  = g4 ξv + g4c ξc;
{A1,A2,A3,A4}
];

Fit data with model 2

(* load data *)
dat = dat2fit ;datsig = dat[[All,1;;3]]; errsig = dat[[All,4]];
ander2 = AndersonPlot2; (* data without joining the points *)
(* 2x2 residues and 2 linear trajectories : 8 parameters *)
nlmγπ2=NonlinearModelFit[datsig,{DiffCrossSect2[pl,t,g1,g4,α0,αp,g1c,g4c,α0c,αpc]},
{{g4,-20 },{g1,1},{α0,0.5},{αp,0.9},{g4c,4 },{g1c,0},{α0c,0.5},{αpc,0.2}},{pl,t}, Weights->1/errsig^2,VarianceEstimatorFunction->(1&),AccuracyGoal->5];

 χ2 =nlmγπ2["EstimatedVariance",VarianceEstimatorFunction->Automatic];
Print["Resulting χ^2 = ",χ2]
nlmγπ2["ParameterTable"]
plotSigB = Show[{ander2,pAfit[nlmγπ2]},leg, FrameLabel->{"t ((GeV^2))", "dσ/dt μb/(GeV^2)"},framestyle]

(* The dimensions of the parameters are 
a4 GeV-3
a2 GeV-4
a1 GeV-2
αp GeV-2 
*)

Export Figures

SetDirectory["/Users/vincent/Dropbox/RTSchool"];
Export["plotAlpEff.pdf",plotAlpEff,ImageSize->Scaled[1.5]];
Export["Anderson1.pdf",plotCS,ImageSize->Scaled[1.5]];
Export["Model1.pdf",plotSigA,ImageSize->Scaled[1.5]];
Export["Model2.pdf",plotSigB,ImageSize->Scaled[1.5]];

Eta Protoproduction

SetDirectory["/Users/vincent/Dropbox/mes articles/Regge/data"];
Braun70 = Module[{sig4,sig6,dat,pp4,pp6,pp},
(* load the data *)
sig4 = Import["Braun70_Sig_PhotEta_4.txt","Table"];
sig6 = Import["Braun70_Sig_PhotEta_6.txt","Table"];
dat = sig4;
pp4 = ErrorListLogPlot[
Table[{{-dat[[i,1]],dat[[i,4]]}, ErrorBar[ dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->{Thick,Orange}];
dat = sig6;
pp6 = ErrorListLogPlot[
Table[{{-dat[[i,1]],dat[[i,4]]}, ErrorBar[ dat[[i,5]]]},{i,1,Dimensions[dat][[1]]}],PlotMarkers->Automatic,PlotStyle->{Thick,Darker[Green]}];
pp = {pp4,pp6};
{pp}];

leg = Graphics[Text[Column[{
Row[{Style["•",Orange,20],Style["4 GeV",20]},Spacer[10]],
Row[{Style["•",Darker[Green],20],Style["6 GeV",20]},Spacer[10]]
}],{-1,-1.7}]];
framestyle = {PlotRange->{{-1.5,0},All},Frame->True, FrameStyle->FontSize->25,Axes->None,GridLines->Automatic};
plotBraun=Show[Braun70,leg,framestyle,FrameLabel->{"t ((GeV^2))", "dσ/dt μb/(GeV^2)"}]

SetDirectory["/Users/vincent/Dropbox/RTSchool"];
Export["plotEta.pdf",plotBraun,ImageSize->Scaled[1.5]];

Sqrt[.5/0.9]