Input

Exit[];

<<NumericalDifferentialEquationAnalysis`;
Needs["NumericalCalculus`"];
Needs["ErrorBarPlots`"]


ACURACY=5;MAXREC=200;
ϵ=10.^(-4);
ϵ0=10.^(-6);
ϵ2=10.^(-5);

LamPhShift=1.3(*GeV*);(* for s>LamPhShift the phase smoothly approaches its asymptotic: made for comparision *)
LamOmnes=1.0;
LamDisp=1.0;
CaseOfstudy="Omega";(*Phi, Omega*)

JJ=1; 
(* Conformal mapping *)
si=LamDisp^2;
sE=0;

(* Experimental input *)
GammaOmegaTot=8.49*10^(-3);
GammaPhiTot=4.26*10^(-3);

GammaOmega3pi=0.892GammaOmegaTot;
GammaPhi3pi=0.153GammaPhiTot;

nOmn=100;
Which[CaseOfstudy=="Omega",M=Momega;Gamma3piExp=GammaOmega3pi;,CaseOfstudy=="Phi",M=Mphi;Gamma3piExp=GammaPhi3pi;];
SetDirectory[NotebookDirectory[]];
<<Styling.m;
<<Masses.m;
<<PiPi_Madrid.m;
<<Omnes.m;
(* delta[J][s]; t[J][s]; for E < 1.42 GeV*)
(* deltaFinal[J][s] - smoothly extrapolated at LamPhShift (J=0,1 only)*)
(* OmnesIϵ[±1][s] (above+1/below-1 cut), Omnes0[J][s] - elsewhere *)

s0=4mpi^2;
a0=1/2(M^2-mpi^2);
a=(M-mpi)^2;
b=(M+mpi)^2;

(*
kI[s_]:=(k0=1/s*Sqrt[Abs[s(s-4mpi^2)(s-(M-mpi)^2)(s-(M+mpi)^2)]]; k0*Which[4*mpi^2≤s<(M-mpi)^2,1,(M-mpi)^2≤s<(M+mpi)^2,I,(M+mpi)^2≤s<Infinity,-1]);
*)

kI[s_]:=Sqrt[(s-4mpi^2)/s] Sqrt[a-s]Sqrt[b-s];
tPI[s_]:=(M^2+3mpi^2-s)/2+kI[s]/2;
tMI[s_]:=(M^2+3mpi^2-s)/2-kI[s]/2;

RV=2/3 M(M-3mpi);
x[s_,t_]:=(t-(M^2+3mpi^2-s-t))/(Sqrt[3]*RV);
y[s_]:=(1/3(M^2+3mpi^2)-s)/RV;

(* Conformal mapping *)
w0[s_]:=(Sqrt[si-sE]-Sqrt[si-s])/(Sqrt[si-sE]+Sqrt[si-s]) ;
ss[w1_]:=(4 LamDisp^2 w1+sE(w1-1)^2)/(((w1+1))^2);
w[sig_][s_]:=If[s==0,(Sqrt[si-sE]-Sqrt[si])/(Sqrt[si-sE]+Sqrt[si]),Which[sig==0, w0[s],sig==1||sig==-1,w0[s+I sig 10^(-9)]]];

fig00[1]=Show[Plot[deltaFinal[0][par^2]*180/Pi,{par,2mpi,2},PlotPoints->50,MaxRecursion->0,PlotStyle->{Black,Thick}],DeltaExp[0],ImageSize->{250,250}];fig00[2]=Plot[{Re[t[0][par^2]],Im[t[0][par^2]]},{par,2mpi,1.4},PlotPoints->50,MaxRecursion->0,PlotStyle-> {{Blue,Thick},{Red,Thick}},PlotRange->All,ImageSize->{250,250}];
fig11[1]=Show[Plot[deltaFinal[1][par^2]*180/Pi,{par,2mpi,2},PlotPoints->50,MaxRecursion->0,PlotStyle->{Black,Thick}],DeltaExp[1],ImageSize->{250,250}];fig11[2]=Plot[{Re[t[1][par^2]],Im[t[1][par^2]]},{par,2mpi,1.2},PlotPoints->50,MaxRecursion->0,PlotStyle-> {{Blue,Thick},{Red,Thick}},PlotRange->All,ImageSize->{250,250}];
GraphicsRow[{fig00[1],fig00[2],fig11[1],fig11[2]}]

GraphicsRow[{Plot[{Re[OmnesPrime[1][par^2]],Re[OmnesEl[1][par^2]],Re[Omnes[1][par^2]]},{par,ϵ,LamDisp+0.1},PlotStyle->{{Gray,Dashed},{Brown,DotDashed},{Black}},MaxRecursion->0,PlotPoints->100,PlotRange->All,ImageSize->{250,250},FrameLabel->{"(Sqrt[s]) [GeV]","Re(Ω(s))"}],Plot[{Im[OmnesPrime[1][par^2]],Im[OmnesEl[1][par^2]],Im[Omnes[1][par^2]]},{par,ϵ,LamDisp+0.1},PlotStyle->{{Gray,Dashed},{Brown,DotDashed},{Black}},MaxRecursion->0,PlotPoints->100,PlotRange->All,ImageSize->{250,250},FrameLabel->{"(Sqrt[s]) [GeV]","Im(Ω(s))"}]}]

Solution of the integral equaiton

(*======*)
Nmax=1;(* Maximum power of w(s) in the conformal mapping expansion *)
(*======*)

(* Starting function *)
FF[0][nxi_][sig_][s_]:=Omnes[sig][s]*If[nxi==0,1,w[sig][s]^nxi];

pl1=Plot[{Re[FF[0][0][+1][par^2]],Im[FF[0][0][+1][par^2]]},{par,ϵ,1.1},PlotStyle->{{Gray},{Black,Dashed},{Black},{Black,Thick}},MaxRecursion->0,PlotPoints->100,ImageSize->{250,250},FrameLabel->{"(Sqrt[s]) [GeV]","(F^((0)))(s)[(Subscript[n, ξ])=0]"},GridLines->{{{LamOmnes,Dotted}},{}}];
pl2=Plot[{Re[FF[0][1][+1][par^2]],Im[FF[0][1][+1][par^2]]},{par,ϵ,1.1},PlotStyle->{{Gray},{Black,Dashed},{Black},{Black,Thick}},MaxRecursion->0,PlotPoints->100,ImageSize->{250,250},FrameLabel->{"(Sqrt[s]) [GeV]","(F^((0)))(s)[(Subscript[n, ξ])=1]"},GridLines->{{{LamOmnes,Dotted}},{}}];
GraphicsRow[{pl1,pl2}]

NumIter=3;

Do[
Do[
FFIntpl[iter-1][nxi][1]=(a1=4mpi^2;b1=0.95;Interpolation[Table[{s,FF[iter-1][nxi][1][s]},{s,a1,b1,Abs[b1-a1]/200}]]);
FFIntpl[iter-1][nxi][-1]=(a1=4mpi^2;b1=0.95;Interpolation[Table[{s,FF[iter-1][nxi][-1][s]},{s,a1,b1,Abs[b1-a1]/200}]]);

TempFun[s_,t_]=Module[{zs},zs=(2t+s-M^2-3mpi^2)/kI[s];3(1-zs^2)*1/kI[s]*(a-s)^(3/2)//Simplify];

FHat[iter_][nxi_][s_]:=Which[(s==s0),2FFIntpl[iter-1][nxi][+1][tMI[s0]](a-s0)^(3/2),(s0< s<a0),NIntegrate[TempFun[s,t]FFIntpl[iter-1][nxi][+1][t],{t,tMI[s],tPI[s]},AccuracyGoal->ACURACY,MaxRecursion->MAXREC],(a0≤ s≤a),NIntegrate[TempFun[s,t]FFIntpl[iter-1][nxi][-1][t],{t,tMI[s],4mpi^2-ϵ2},AccuracyGoal->ACURACY,MaxRecursion->MAXREC]+NIntegrate[TempFun[s,t]FFIntpl[iter-1][nxi][+1][t],{t,4mpi^2-ϵ2,tPI[s]},AccuracyGoal->ACURACY,MaxRecursion->MAXREC],
(a <s<b),NIntegrate[TempFun[s,t]FF[iter-1][nxi][0][t],{t,tMI[s],2tMI[(M+mpi)^2]},AccuracyGoal->ACURACY,MaxRecursion->MAXREC]+NIntegrate[TempFun[s,t]FF[iter-1][nxi][0][t],{t,2tPI[(M+mpi)^2],tPI[s]},AccuracyGoal->ACURACY,MaxRecursion->MAXREC],
(b≤s),NIntegrate[TempFun[s,t]FF[iter-1][nxi][0][t],{t,tMI[s],tPI[s]},AccuracyGoal->ACURACY,MaxRecursion->MAXREC]];

Mtilde[iter_][nxi_][sb_]:=Sin[deltaFinal[1][sb]]/Abs[Omnes[1][sb]]FHat[iter][nxi][sb];(* Needed only for real s'=[4m^2, LamDisp^2] -> interpolate *)


(**)
(* Store the results for Mtilde, once you did it - comment this block *)
Export["V3p_final_ListMtilde_"<>ToString[CaseOfstudy]<>ToString["_nXi"]<>ToString[nxi]<>ToString["_iter"]<>ToString[iter]<>"_s0_a.m",Table[{sb,Mtilde[iter][nxi][sb]},{sb,Flatten[{Table[i,{i,s0,a-0.011,(a-0.011-s0)/100}],Table[i,{i,a-0.01,a,0.001}]}]}]];
Export["V3p_final_ListMtilde_"<>ToString[CaseOfstudy]<>ToString["_nXi"]<>ToString[nxi]<>ToString["_iter"]<>ToString[iter]<>"_a_Lambda.m",Table[{sb,Mtilde[iter][nxi][sb]},{sb,Flatten[{Table[i,{i,a+ϵ,a+0.01,0.001}],Table[i,{i,a+0.011,LamDisp^2-0.006,(LamDisp^2-0.006-a-0.011)/100}],Table[i,{i,LamDisp^2-0.005,LamDisp^2,0.01/10}]}]}]];
(**)


Do[M1[iter][nxi][s_]=Interpolation[Import["V3p_final_ListMtilde_"<>ToString[CaseOfstudy]<>ToString["_nxi"]<>ToString[nxi]<>ToString["_iter"]<>ToString[iter]<>"_s0_a.m"],InterpolationOrder->3][s];
M2[iter][nxi][s_]=Interpolation[Import["V3p_final_ListMtilde_"<>ToString[CaseOfstudy]<>ToString["_nxi"]<>ToString[nxi]<>ToString["_iter"]<>ToString[iter]<>"_a_Lambda.m"],InterpolationOrder->3][s];
,{iter,1,NumIter}];
MtildeIntpl[iter_][nxi_][s_]:=Which[s0≤ s≤a,Chop[M1[iter][nxi][s]],a< s≤LamDisp^2,Chop[M2[iter][nxi][s]]];


FFtemp[iter_][nxi_][sig_][s_]:=Omnes[sig][s](If[nxi==0,1,w[sig][s]^nxi]+1/π NIntegrate[MtildeIntpl[iter][nxi][sb]/((a-sb+I ϵ*10)^(3/2)(sb-s-sig I ϵ)),{sb,s0,a,LamDisp^2},AccuracyGoal->ACURACY,MaxRecursion->MAXREC]-1/π MtildeIntpl[iter][nxi][LamDisp^2]/(a-LamDisp^2)^(3/2) Log[(LamDisp^2-s-sig I ϵ)/(LamDisp^2-4mpi^2)](*this term regularize the cutoff point*));

FF[iter_][nxi_][sig_][s_]:=Which[sig==0,(FFtemp[iter][nxi][0][s]),sig==1||sig==-1,(Which[s≤a-0.05||s≥a+0.05,FFtemp[iter][nxi][sig][s],a-0.05<s<a+0.05,Interpolation[{{a-0.055,FFtemp[iter][nxi][sig][a-0.055]},{a-0.05,FFtemp[iter][nxi][sig][a-0.05]},{a+0.05,FFtemp[iter][nxi][sig][a+0.05]},{a+0.055,FFtemp[iter][nxi][sig][a+0.055]}}][s]])];

,{iter,1,NumIter}];
,{nxi,0,Nmax}];

plot1=Plot[{Re[FF[0][(*nxi*)0][+1][par^2]],Re[FF[1][(*nxi*)0][+1][par^2]],Re[FF[2][(*nxi*)0][+1][par^2]],Re[FF[3][(*nxi*)0][+1][par^2]]},{par,2mpi,LamDisp+0.1},PlotStyle->{{Gray},{Black,Dotted},{Black,Dashed},{Black},Blue},MaxRecursion->0,PlotPoints->100,PlotRange->All,ImageSize->{250,250}];
plot2=Plot[{Im[FF[0][(*nxi*)0][+1][par^2]],Im[FF[1][(*nxi*)0][+1][par^2]],Im[FF[2][(*nxi*)0][+1][par^2]],Im[FF[3][(*nxi*)0][+1][par^2]]},{par,2mpi,LamDisp+0.1},PlotStyle->{{Gray},{Black,Dotted},{Black,Dashed},{Black},Blue},MaxRecursion->0,PlotPoints->100,PlotRange->All,ImageSize->{250,250}];

plot3=Plot[{Re[FF[0][(*nxi*)1][+1][par^2]],Re[FF[1][(*nxi*)1][+1][par^2]],Re[FF[2][(*nxi*)1][+1][par^2]],Re[FF[3][(*nxi*)1][+1][par^2]]},{par,2mpi,LamDisp+0.1},PlotStyle->{{Gray},{Black,Dotted},{Black,Dashed},{Black},Blue},MaxRecursion->0,PlotPoints->100,PlotRange->All,ImageSize->{250,250}];
plot4=Plot[{Im[FF[0][(*nxi*)1][+1][par^2]],Im[FF[1][(*nxi*)1][+1][par^2]],Im[FF[2][(*nxi*)1][+1][par^2]],Im[FF[3][(*nxi*)1][+1][par^2]]},{par,2mpi,LamDisp+0.1},PlotStyle->{{Gray},{Black,Dotted},{Black,Dashed},{Black},Blue},MaxRecursion->0,PlotPoints->100,PlotRange->All,ImageSize->{250,250}];

GraphicsRow[{plot1,plot2,plot3,plot4}]

Store results and then interpolate

iter=NumIter;(* Number of iteration in the amplitude *)

(**)
(* Store the results for F, once you did it - comment this block *)
Do[Export["V3p_final_ListF_"<>ToString[CaseOfstudy]<>ToString["_nxi"]<>ToString[nxi]<>ToString["_iter"]<>ToString[iter]<>".m",Table[{sb,FF[iter][nxi][1][sb]},{sb,0,(LamDisp+0.1)^2,Abs[(LamDisp+0.1)^2]/200}]];

Export["V3p_final_ListF_"<>ToString[CaseOfstudy]<>ToString["_nxi"]<>ToString[nxi]<>ToString["_iter"]<>ToString[0]<>".m",Table[{sb,FF[0][nxi][1][sb]},{sb,0,(LamDisp+0.1)^2,Abs[(LamDisp+0.1)^2]/200}]];
,{nxi,0,Nmax}];

Export["V3p_final_ListF_"<>ToString[CaseOfstudy]<>ToString["_iter"]<>ToString[iter]<>".dat",Table[{sb,Re[FF[iter][0(*nxi*)][+1][sb]],Im[FF[iter][0(*nxi*)][+1][sb]],Re[FF[iter][1(*nxi*)][+1][sb]],Im[FF[iter][1(*nxi*)][+1][sb]]},{sb,0,(LamDisp+0.1)^2,Abs[(LamDisp+0.1)^2]/200}]];
Export["V3p_final_ListF_"<>ToString[CaseOfstudy]<>ToString["_iter"]<>ToString[0]<>".dat",Table[{sb,Re[FF[0][0(*nxi*)][+1][sb]],Im[FF[0][0(*nxi*)][+1][sb]],Re[FF[0][1(*nxi*)][+1][sb]],Im[FF[0][1(*nxi*)][+1][sb]]},{sb,0,(LamDisp+0.1)^2,Abs[(LamDisp+0.1)^2]/200}]];

Do[Export["V3p_final_ListFHat_"<>ToString[CaseOfstudy]<>ToString["_nxi"]<>ToString[nxi]<>ToString["_iter"]<>ToString[iter]<>"_s0_a.m",Table[{sb,FHat[iter][nxi][sb]},{sb,Flatten[{Table[i,{i,s0,a-0.011,(a-0.011-s0)/100}],Table[i,{i,a-0.01,a,0.001}]}]}]];
Export["V3p_final_ListFHat_"<>ToString[CaseOfstudy]<>ToString["_nxi"]<>ToString[nxi]<>ToString["_iter"]<>ToString[iter]<>"_a_Lambda.m",Table[{sb,FHat[iter][nxi][sb]},{sb,Flatten[{Table[i,{i,a+ϵ,a+0.01,0.001}],Table[i,{i,a+0.011,LamDisp^2-0.006,(LamDisp^2-0.006-a-0.011)/100}],Table[i,{i,LamDisp^2-0.005,LamDisp^2,0.01/10}]}]}]];
,{nxi,0,Nmax}];
(**)
*
(* Interpolate the final solution *)
Do[
FIntpl[iter][nxi][1][s_]=Interpolation[Import["V3p_final_ListF_"<>ToString[CaseOfstudy]<>ToString["_nxi"]<>ToString[nxi]<>ToString["_iter"]<>ToString[iter]<>".m"],InterpolationOrder->3][s];
FIntpl[0][nxi][1][s_]=Interpolation[Import["V3p_final_ListF_"<>ToString[CaseOfstudy]<>ToString["_nxi"]<>ToString[nxi]<>ToString["_iter"]<>ToString[0]<>".m"],InterpolationOrder->3][s];

FHat1[iter][nxi][s_]=Interpolation[Import["V3p_final_ListFHat_"<>ToString[CaseOfstudy]<>ToString["_nxi"]<>ToString[nxi]<>ToString["_iter"]<>ToString[iter]<>"_s0_a.m"],InterpolationOrder->3][s];
FHat2[iter][nxi][s_]=Interpolation[Import["V3p_final_ListFHat_"<>ToString[CaseOfstudy]<>ToString["_nxi"]<>ToString[nxi]<>ToString["_iter"]<>ToString[iter]<>"_a_Lambda.m"],InterpolationOrder->3][s];
,{nxi,0,Nmax}];

FHatIntpl[iter][nxi_][s_]:=Which[s0≤ s≤a,Chop[FHat1[iter][nxi][s]],a< s≤LamDisp^2,Chop[FHat2[iter][nxi][s]]];
f1[nxi_][s_]:=FIntpl[iter][nxi][1][s]+FHatIntpl[iter][nxi][s]/(a-s)^(3/2);

(* Plots *)
plot1=Plot[{Re[FHatIntpl[iter][(*nxi*)0][s]],Im[FHatIntpl[iter][(*nxi*)0][s]]},{s,s0,LamDisp^2},PlotStyle->{{Gray},{Black,Dotted},{Black,Dashed},{Black},Blue},MaxRecursion->0,PlotPoints->100,PlotRange->All,ImageSize->{250,250}];
plot2=Plot[{Re[FIntpl[0][(*nxi*)0][+1][par^2]],Re[FIntpl[iter][(*nxi*)0][+1][par^2]]},{par,0,LamDisp+0.01},PlotStyle->{{Gray},{Black,Dotted},{Black,Dashed},{Black},Blue},MaxRecursion->0,PlotPoints->100,PlotRange->All,ImageSize->{250,250}];plot3=Plot[{Im[FIntpl[0][(*nxi*)0][+1][par^2]],Im[FIntpl[iter][(*nxi*)0][+1][par^2]]},{par,0,LamDisp+0.01},PlotStyle->{{Gray},{Black,Dotted},{Black,Dashed},{Black},Blue},MaxRecursion->0,PlotPoints->100,PlotRange->All,ImageSize->{250,250}];

(* nxi=0 *)
plot4=Plot[{Arg[FIntpl[0][(*nxi*)0][+1][par^2]],Arg[f1[(*nxi*)0][par^2]]},{par,2mpi,LamDisp},PlotStyle->{{Gray},{Black},{Black,Dashed}},MaxRecursion->0,PlotPoints->350,ImageSize->{250,250}];
plot5=Plot[{Abs[FIntpl[0][(*nxi*)0][+1][par^2]],Abs[f1[(*nxi*)0][par^2]]},{par,2mpi,LamDisp},PlotStyle->{{Gray},{Black},{Black,Dashed}},MaxRecursion->0,PlotPoints->350,ImageSize->{250,250}];

(* nxi=1 *)
plot6=Plot[{Arg[FIntpl[0][(*nxi*)1][+1][par^2]],Arg[f1[(*nxi*)1][par^2]]},{par,2mpi,LamDisp},PlotStyle->{{Gray},{Black},{Black,Dashed}},MaxRecursion->0,PlotPoints->350,ImageSize->{250,250}];
plot7=Plot[{Abs[FIntpl[0][(*nxi*)1][+1][par^2]],Abs[f1[(*nxi*)1][par^2]]},{par,2mpi,LamDisp},PlotStyle->{{Gray},{Black},{Black,Dashed}},MaxRecursion->0,PlotPoints->350,ImageSize->{250,250}];
GraphicsGrid[{{plot1,plot2,plot3},{plot4,plot5,plot6,plot7}}]

Dalitz plots

CaseOfstudy2="3b";(*3b, No3b*)

Ffinal[iter_][1][s_]:=Which[Nmax==0,N0 FIntpl[iter][0(*nxi*)][+1][s],Nmax==1,N0(FIntpl[iter][0(*nxi*)][+1][s]+param FIntpl[iter][1(*nxi*)][+1][s])];

Which[CaseOfstudy2=="No3b",F[s_,t_]:=Ffinal[0][+1][s]+Ffinal[0][+1][t]+Ffinal[0][+1][3mpi^2+M^2-s-t];,CaseOfstudy2=="3b",F[s_,t_]:=Ffinal[iter][+1][s]+Ffinal[iter][+1][t]+Ffinal[iter][+1][3mpi^2+M^2-s-t];];

P[s_,t_]:=1/4(s t (3mpi^2+M^2-s-t)-mpi^2(M^2-mpi^2)^2);
M2[s_,t_]:=ComplexExpand[Re[F[s,t]]]^2+ComplexExpand[Im[F[s,t]]]^2;(*Abs[F[s,t]]^2*)

E2st[s_]:=Sqrt[s]/2;E3st[s_]:=(M^2-s-mpi^2)/(2Sqrt[s]);
tmin[s_]:=(E2st[s]+E3st[s])^2-(Sqrt[E2st[s]^2-mpi^2]+Sqrt[E3st[s]^2-mpi^2])^2;
tmax[s_]:=(E2st[s]+E3st[s])^2-(Sqrt[E2st[s]^2-mpi^2]-Sqrt[E3st[s]^2-mpi^2])^2;

E1st[t_]:=(M^2-t-mpi^2)/(2*Sqrt[t]);
Ett[t_]:=Sqrt[t]/2;
smin[t_]:=(E1st[t]+Ett[t])^2-(Sqrt[E1st[t]^2-mpi^2]+Sqrt[Ett[t]^2-mpi^2])^2;
smax[t_]:=(E1st[t]+Ett[t])^2-(Sqrt[E1st[t]^2-mpi^2]-Sqrt[Ett[t]^2-mpi^2])^2;

(* In terms of X and Y *)
Ymin=-1;
Ymax=(M+3mpi)/(2M);

X1min[Y_]:=-((Sqrt[-9*mpi^2+M*(M-2*M*Y+6*mpi*Y)]*Sqrt[(M^2*(M-3*mpi)*(1+Y)*(M+3*mpi+M*Y-3*mpi*Y))/(3*mpi^2+M*(M-2*M*Y+6*mpi*Y))])/(Sqrt[3]*M*(M-3*mpi)));
X1max[Y_]:=((Sqrt[-9*mpi^2+M*(M-2*M*Y+6*mpi*Y)]*Sqrt[(M^2*(M-3*mpi)*(1+Y)*(M+3*mpi+M*Y-3*mpi*Y))/(3*mpi^2+M*(M-2*M*Y+6*mpi*Y))])/(Sqrt[3]*M*(M-3*mpi)));
(*x[s1[Y],tmin[s1[Y]]]*)

Xmin=X1min[Y1]/.Solve[D[X1min[Y1],Y1]==0][[2]];
Xmax=-Xmin;

s1[y_]:=1/3 (M^2(1-2 y)+6 M mpi y+3 mpi^2);
t1[x_,y_]:=1/3 (M^2 (Sqrt[3] x+y+1)-3 M mpi (Sqrt[3] x+y)+3 mpi^2);
P1[x_,y_]:=(M^2*(M-3*mpi)^2*(-9*mpi^2*(-1+x^2+y^2)+6*M*mpi*(1-3*x^2*y+y^3)-M^2*(-1+2*y)*(-3*x^2+(1+y)^2)))/108;


Y1min[x_]:=y1/.Solve[P1[x,y1]==0,y1][[2]];
Y1max[x_]:=y1/.Solve[P1[x,y1]==0,y1][[3]];


(* Boundaries of Dalitz Plot *)
plotBoundary=ListLinePlot[Join[Table[{X1min[Y],Y},{Y,Ymin,Ymax,0.01}],Reverse[Table[{X1max[Y],Y},{Y,Ymin,Ymax,0.01}]]],PlotStyle->{Gray,Thick}];

(**)
dGammaS[s_]:=NIntegrate[1/(2π)^3*1/(32 M^3)*1/3*P[s,t]*M2[s,t],{t,tmin[s],tmax[s]}];
dGammaT[t_]:=NIntegrate[1/(2π)^3*1/(32 M^3)*1/3*P[s,t]*M2[s,t],{s,smin[t],smax[t]}];

dGammaY[y_]:=NIntegrate[1/(2π)^3*1/(32 M^3)*1/3*P[s1[y],t1[x,y]]*M2[s1[y],t1[x,y]],{x,X1min[y],X1max[y]}];
dGammaX[x_]:=NIntegrate[1/(2π)^3*1/(32 M^3)*1/3*P[s1[y],t1[x,y]]*M2[s1[y],t1[x,y]],{y,Ymin,Ymax}];

TotalGammaST:=NIntegrate[dGammaS[s],{s,4mpi^2,(M-mpi)^2}];
TotalGammaTS:=NIntegrate[dGammaT[t],{t,4mpi^2,(M-mpi)^2}];


(* Differential decay width dΓ/dx*)
dGammaX2[x_]:=Module[{ny,wgy,yn,wgyn},
ny=25;
wgy=GaussianQuadratureWeights[ny,Y1min[x],Y1max[x]]; (*includes points and weights*)
yn=Table[wgy[[i,1]],{i,ny}];
wgyn=Table[wgy[[i,2]],{i,ny}];
Sum[1/(2π)^3*1/(32M^3)*1/3*P[s1[yn[[i]]],t1[x,yn[[i]]]]*M2[s1[yn[[i]]],t1[x,yn[[i]]]]wgyn[[i]],{i,1,ny}]];

(* Differential decay width dΓ/dx*)
dGammaY2[y_]:=Module[{nx,wgx,xn,wgxn},
nx=25;
wgx=GaussianQuadratureWeights[nx,X1min[y],X1max[y]]; (*includes points and weights*)
xn=Table[wgx[[i,1]],{i,nx}];
wgxn=Table[wgx[[i,2]],{i,nx}];
Sum[1/(2π)^3*1/(32M^3)*1/3*P[s1[y],t1[xn[[i]],y]]*M2[s1[y],t1[xn[[i]],y]]wgxn[[i]],{i,1,nx}]];


(* Differential decay width *)
dGammaS2[s_]:=Module[{nt,wgt,tn,wgtn},
nt=25;
wgt=GaussianQuadratureWeights[nt,tmin[s],tmax[s]]; (*includes points and weights*)
tn=Table[wgt[[i,1]],{i,nt}];
wgtn=Table[wgt[[i,2]],{i,nt}];
Sum[1/(2π)^3*1/(32M^3)*1/3*P[s,tn[[i]]]*M2[s,tn[[i]]]wgtn[[i]],{i,1,nt}]];

(* decay width *)
ns=25; 
wgs=GaussianQuadratureWeights[ns,4mpi^2,(M-mpi)^2]; (*include points and weights*)
sn=Table[wgs[[i,1]],{i,ns}];
wgsn=Table[wgs[[i,2]],{i,ns}];
TotalGamma2:=Sum[dGammaS2[sn[[i]]]wgsn[[i]],{i,1,ns}];

TotalGamma2temp[N0_,param_]=TotalGamma2//Simplify//Chop

Print["Fix (Subscript[a, 0]) coefficient from the decay width in case (Subscript[n, ξ])=0"]
Print[CaseOfstudy2"1 term in ξ-pansion: "]
Nfixed=Sqrt[Gamma3piExp/TotalGamma2]/.{N0->1,param->0}
Nfixed^2*TotalGamma2*10^3/.{N0->1,param->0}(*MeV*)(* Cross-Check: 7.57 - Omega; 0.65 - Phi *);
(**)
Export["V3p_final_aCoeff_"<>ToString[CaseOfstudy]<>ToString["_iter"]<>ToString[iter]<>".m",Nfixed];

Next figures plotted for one term in the conformal expnasion ((Subscript[a, 0])) which is fixed from the (Subscript[Γ, V-> 3π])

Print["  "CaseOfstudy<>", "<>CaseOfstudy2]

plot1=Plot[Which[CaseOfstudy2=="No3b",Re[Ffinal[0][1][par^2]]/.{N0->Nfixed,param->0},CaseOfstudy2=="3b",Re[Ffinal[iter][1][par^2]]]/.{N0->Nfixed,param->0},{par,0,1},GridLines->{{2mpi,(M-mpi)},{}},PlotStyle->{Black,Which[CaseOfstudy2=="No3b",Dashed,CaseOfstudy2=="3b",Thick]},FrameLabel->{"(Sqrt[s]) [GeV]","Re[F(s)]"},MaxRecursion->0,PlotPoints->100,PlotRange->Which[CaseOfstudy=="Omega",{-1500,1200},CaseOfstudy=="Phi",{-100,100}],ImageSize->{250,250}];
plot2=Plot[Which[CaseOfstudy2=="No3b",Im[Ffinal[0][1][par^2]]/.{N0->Nfixed,param->0},CaseOfstudy2=="3b",Im[Ffinal[iter][1][par^2]]]/.{N0->Nfixed,param->0},{par,0,1},GridLines->{{2mpi,(M-mpi)},{}},PlotStyle->{Black,Which[CaseOfstudy2=="No3b",Dashed,CaseOfstudy2=="3b",Thick]},FrameLabel->{"(Sqrt[s]) [GeV]","Im[F(s)]"},MaxRecursion->0,PlotPoints->100,PlotRange->Which[CaseOfstudy=="Omega",{-300,2200},CaseOfstudy=="Phi",{-50,150}],ImageSize->{250,250}];

plot3=Plot[dGammaS2[s]/.{N0->Nfixed,param->0},{s,4mpi^2,(M-mpi)^2},PlotStyle->{Black,Which[CaseOfstudy2=="No3b",Dashed,CaseOfstudy2=="3b",Thick]},FrameLabel->{"s [(GeV^2)]","dΓ/ds"},PlotPoints->40,MaxRecursion->0,PlotRange->Which[CaseOfstudy=="Omega",{0,0.05},CaseOfstudy=="Phi",{0,0.0025}],ImageSize->{250,250}];

plot4=Plot[dGammaX2[x]/.{N0->Nfixed,param->0},{x,Xmin,Xmax},PlotStyle->{Black,Which[CaseOfstudy2=="No3b",Dashed,CaseOfstudy2=="3b",Thick]},FrameLabel->{"x","dΓ/dx"},PlotPoints->40,MaxRecursion->0,PlotRange->Which[CaseOfstudy=="Omega",{0,0.3},CaseOfstudy=="Phi",{0,0.006}],ImageSize->{250,250}];

plot5=Plot[dGammaY2[y]/.{N0->Nfixed,param->0},{y,Ymin,Ymax},PlotStyle->{Black,Which[CaseOfstudy2=="No3b",Dashed,CaseOfstudy2=="3b",Thick]},FrameLabel->{"y","dΓ/dy"},PlotPoints->40,MaxRecursion->0,PlotRange->Which[CaseOfstudy=="Omega",{0,0.3},CaseOfstudy=="Phi",{0,0.006}],ImageSize->{250,250}];

GraphicsGrid[{{plot1,plot2,plot3},{plot4,plot5}}]

colorD=ColorData[{"PigeonTones","Reverse"}];

Dalitz[X_,Y_,N1_,param1_]:=M2[s1[Y],t1[X,Y]]/.{N0->N1,param->param1};
DalitzNorm[X_,Y_,N1_,param1_]:=Dalitz[X,Y,N1,param1]/Dalitz[0,0,N1,param1];

t=Flatten[Table[{X,Y,DalitzNorm[X,Y,Nfixed,100]},{Y,Ymin+ϵ0,Ymax-ϵ0,0.01},{X,X1min[Y]+ϵ0,X1max[Y]-ϵ0,0.01}],1];

{contourPlot,contourColors,contourRange}=reportColorRange[ListContourPlot[t,FrameLabel->{"x","y"}]];
densityPlot=ListDensityPlot[t,FrameLabel->{"x","y"},ColorFunction->colorD];

Block[{labelWidth=.15,labelHeight=0.8,aspectRatio=0.8,numberOfContours=5,plotRange},plotRange=Round[Table[contourRange[[1]]+i(contourRange[[2]]-contourRange[[1]])/(numberOfContours),{i,0,numberOfContours}],0.05];numberOfContours=Length[plotRange]-1;

display[{densityPlot(*contourPlot*)//at[{0,0},aspectRatio],colorLegend[contourColors,plotRange[[{1,-1}]],LabelStyle->Directive[12],Background->None,"ColorSwathes"->numberOfContours]//at[{1-labelWidth,0},{labelWidth,1}]},AspectRatio->aspectRatio,ImageSize->400,ColorFunction->colorD]]

(*Export["1.pdf",a1,ImageResolution->350];*)

*Omega, phi

Dalitz plot parameters

(* Polar coordinates z,ϕ *)
x2[z_,ϕ_]:=(Sqrt[z])Cos[ϕ];y2[z_,ϕ_]:=(Sqrt[z])Sin[ϕ];
P2[z_,ϕ_]:=(M^2*(M-3*mpi)^2*(6*M*mpi+M^2*(1-3*z)-9*mpi^2*(-1+z)+2*M*(M-3*mpi)*z^(3/2)*Sin[3*ϕ]))/108;
Zmax[ϕ_?NumberQ]:=z0/.Solve[z0==X1max[Sqrt[z0] Sin[ϕ]]^2/Cos[ϕ]^2,z0][[1,1]];
ZmaxIntpl=(a1=ϵ0;b1=2π-ϵ0;Interpolation[Table[{ϕ,Zmax[ϕ]},{ϕ,a1,b1,Abs[b1-a1]/201}]]);
Nd=NIntegrate[1,{ϕ,10^(-7),2π},{z,0,ZmaxIntpl[ϕ]}];

(*αtable=α0*10^3*)
Abs2Fpol[n_,α0_,β0_,γ0_,δ0_][z_,ϕ_]:=n^2(1+2α0*10.^(-3) z +2 β0 *10.^(-3)z^(3/2)Sin[3ϕ]+2γ0*10.^(-3) z^2+2δ0*10.^(-3) z^(5/2)Sin[3ϕ]);

chi2[n_,αt_,βt_,γt_,δt_]:=1/Nd*NIntegrate[(P2[z,ϕ]/P2[0,0])^2((Abs2Fpol[n,αt,βt,γt,δt][z,ϕ]-Abs[F[s1[y2[z,ϕ]],t1[x2[z,ϕ],y2[z,ϕ]]]]^2)/n^2)^2,{ϕ,ϵ0,2π-ϵ0},{z,0,ZmaxIntpl[ϕ]}]/.{N0->Nfixed,param->0};

(* Differential decay width *)
dchi2Gauss[ϕ_]:=Module[{nz,wgz,zn,wgzn},
nz=25;
wgz=GaussianQuadratureWeights[nz,0,ZmaxIntpl[ϕ]]; (*include points and weights*)
zn=Table[wgz[[i,1]],{i,nz}];
wgzn=Table[wgz[[i,2]],{i,nz}];Sum[1/Nd*(P2[zn[[i]],ϕ]/P2[0,0])^2((Abs2Fpol[n,αt ,βt,γt,δt][zn[[i]],ϕ]-Abs[F[s1[y2[zn[[i]],ϕ]],t1[x2[zn[[i]],ϕ],y2[zn[[i]],ϕ]]]]^2)/n^2)^2(**)wgzn[[i]],{i,1,nz}]]/.{N0->Nfixed,param->0};

chi2Gauss=Module[{nϕ,wgϕ,ϕn,wgϕn},
nϕ=50;
wgϕ=GaussianQuadratureWeights[nϕ,ϵ0,2π-ϵ0]; (*include points and weights*)
ϕn=Table[wgϕ[[i,1]],{i,nϕ}];
wgϕn=Table[wgϕ[[i,2]],{i,nϕ}];
Sum[dchi2Gauss[ϕn[[i]]]wgϕn[[i]],{i,1,nϕ}]
];
chi2Fit[n_,αt_,βt_,γt_,δt_]=chi2Gauss;

Print["  "CaseOfstudy<>", "<>CaseOfstudy2]
FindMinimum[(Sqrt[chi2Fit[Nf,αf,0,0,0]])*10^3,{{Nf,1440},{αf,90}}]
FindMinimum[(Sqrt[chi2Fit[Nf,αf,βf,0,0]])*10^3,{{Nf,1440},{αf,150},{βf,30}}]
FindMinimum[(Sqrt[chi2Fit[Nf,αf,βf,γf,0]])*10^3,{{Nf,1440},{αf,150},{βf,30},{γf,5}}]
FindMinimum[(Sqrt[chi2Fit[Nf,αf,βf,γf,δf]])*10^3,{{Nf,1440},{αf,150},{βf,30},{γf,5},{δf,10}}]

*Fhat=0

*full