aroot=NotebookDirectory[]

(* Set plot defaults *)
myPlotOptions={PlotRange->All,Frame->True,BaseStyle->{FontSize->15},ImageSize->{500}};
SetOptions[ListPlot,myPlotOptions];
SetOptions[Graphics,myPlotOptions];
SetOptions[ParametricPlot,myPlotOptions];
myPlotOptions={PlotRange->All,Frame->True,BaseStyle->{FontSize->15},ImageSize->{500}};
SetOptions[ListLogPlot,myPlotOptions];
SetOptions[Histogram,myPlotOptions];
SetOptions[ListDensityPlot,myPlotOptions];

(*Commonly used constants in SI units*)
cLightN=299792458.;
mec2N=0.510998928 (10^6);

Astra

From the Astra manual:

AstraColumn["x"] = 1;
AstraColumn["y"] =2;
AstraColumn["z"] = 3;
AstraColumn["px"] = 4;
AstraColumn["py"] = 5;
AstraColumn["pz"] = 6;
AstraColumn["t"] =7;
AstraColumn["Q"] = 8;
AstraColumn["index"] = 9;
AstraColumn["status"] = 10;

ImportAstraParticles[particlefile_, verbose_:False]:=Module[{dat, ref},
Check[dat=Import[particlefile, "Table"];, Print["Problem with "<>particlefile];Abort[]];
ref = dat[[1]];
If[ verbose, Print["Reference: ", ref]];
{#[[1]], #[[2]],ref[[3]] +  #[[3]], #[[4]], #[[5]], ref[[6]] +#[[6]], (10^(-9))(ref[[7]] +#[[7]] ),(10^(-9)) #[[8]], #[[9]], #[[10]]}&/@ + dat[[2;;]]

];

BmadParticleFromAstraParticle[{x_, y_, z_, cpx_, cpy_, cpz_, t_, Q_}, z0_,cp0_, t0_, mc2_:mec2N]:=Module[{Δz,Δt, β, cp, energy},
(*Drift by Δz to place particle exactly at z0*)
Δz = z0-z;
cp = Sqrt[(cpx^2)+(cpy^2)+(cpz^2)];
energy = Sqrt[(cpx^2)+(cpy^2)+(cpz^2) + (mc2^2)];
Δt = Δz energy / cpz / cLightN;

(* x, px/p0, y, py/p0, -βc(t-t0), δ = p/p0 -1, charge *)
{x + Δz (cpx/cpz), (cpx/cp0), y + Δz (cpy/cpz),  (cpy/cp0), ((-cp)/energy)cLightN( (*ignore particle's time: t + *)Δt),(cp/cp0)-1, Abs[Q], 1(* Status *)}
];
BmadBunchFromAstraBunch[abunch_, mc2_:mec2N]:=Module[{ref,z0, cp0, t0},
ref = Mean/@Transpose[abunch];
z0 = ref[[3]];
cp0 = Sqrt[((ref[[4]])^2) + ((ref[[5]])^2)+((ref[[6]])^2)];
t0 = ref[[7]];

BmadParticleFromAstraParticle[#[[1;;8]], z0, cp0, t0, mc2]&/@abunch


];




Clear[ImportAstraEmittance];
ImportAstraEmittance[emitfilebasename_]:=Module[{ex, ey, ez, fx, fy, fz},
fx=emitfilebasename<>".Xemit.001";
fy=emitfilebasename<>".Yemit.001";
fz=emitfilebasename<>".Zemit.001";
Check[ex=Import[fx, "Table"];, Print["Problem with "<>fx];Abort[]];
Check[ey=Import[fy, "Table"];, Print["Problem with "<>fy];Abort[]];Check[ez=Import[fz, "Table"];, Print["Problem with "<>fz];Abort[]];
Print["Imported: "<>fx];
Print["Imported: "<>fy];
Print["Imported: "<>fz];
{ex, ey, ez}
];

AstraEmitArray[property_, emitdat_]:=Module[{col, dim, units, scale},
{col, dim, units, scale}=Switch[property, 

(* From X file *)
"z", {1,1, "m", 1},
"t", {2,1, "ns", 1},
"mean_x", {3,1, "mm", 1},
"rms_x", {4,1,  "μm", 1000.},
"rms_px", {5, 1, "μrad", 1000.},
"emit_x", {6,1, "mm-mrad", 1}, 
"corr_xpx", {7, 1, "mrad (?)", 1}, 

(* From Y file *)
"mean_y", {3,2, "mm", 1},
"rms_y", {4, 2,  "μm", 1000.},
"rms_py", {5, 2,  "μrad", 1000.},
"emit_y", {6,2, "mm-mrad", 1}, 
"corr_ypy", {7, 2, "mrad (?)", 1},

(* From Z file *)
"kinetic_energy", {3, 3, "(Subscript[E, kinetic])(MeV)", 1},
"rms_z", {4, 3, "μm", 1000.},
_, Print["Bad property_: "<>property];Abort[];
];

{scale*emitdat[[dim]][[All, col]], units}



];


Clear[PlotAstraEmit];
PlotAstraEmit[yVar_, xVar_, emitdat_, opts:OptionsPattern[]]:=Module[{xlist, xunits, ylist, yunits, xlabel, ylabel},
{xlist, xunits} = AstraEmitArray[xVar, emitdat];
{ylist, yunits} = AstraEmitArray[yVar, emitdat];

xlabel = xVar<>" ("<>xunits<>")";
ylabel = yVar<>" ("<>yunits<>")";

ListPlot[Transpose[{xlist, ylist}],Evaluate[FilterRules[{opts}, Options[ListPlot]]],  Joined->True, FrameLabel->{xlabel, ylabel}, PlotStyle->Orange, PlotRange->All]

]



Clear[AstraGetScaleAndLabel];
AstraGetScaleAndLabel[PSvariable_String]:=Module[{},
Switch[PSvariable,
(*For particles *)
	"x", {(10^3), "x (mm)"}, 
	"px",{ 1, "c (Subscript[p, x]) (eV)"},
	"y",{(10^3), "y (mm)"},
	"py",{ 1, "c (Subscript[p, y]) (eV)"},
	"z",{1, "z (m)"},
	"pz", { (1/(10^6))1, "c (Subscript[p, z]) (MeV)"},
	"t",   {(10^9), "t (ns)"},
	"r", {(10^3), "r (mm)"}, 
	"pr",{ 1, "c (Subscript[p, r]) (eV)"},
	_, Print["Bad variable: "<>PSvariable];Abort[]
]
];


Clear[AstraArray];
(* Get a single summary column, or a function of columns*)
AstraArray[xVar_String, astradat_]:=Module[{f},
AstraGetScaleAndLabel[xVar]; (*just used to check*)
f = Switch[xVar,
(* Convert to from mcβγ to βγ *)
"px",(1)  #[[AstraColumn["px"]]]& ,
"py",(1)  #[[AstraColumn["py"]]]& ,
"pz",(1)  #[[AstraColumn["pz"]]]& ,
"r",(1)  Sqrt[((#[[AstraColumn["x"]]])^2)+((#[[AstraColumn["y"]]])^2)]& ,
"pr",(1)  Sqrt[((#[[AstraColumn["px"]]])^2)+((#[[AstraColumn["py"]]])^2)]& ,
_,  #[[AstraColumn[xVar]]]& (*defaut*)
];
f/@astradat

];



Clear[AstraPhaseSpace];
AstraPhaseSpace[xVar_,yVar_,astraparticles_, opts:OptionsPattern[]]:= Module[{xlist, ylist, yscale, ylabel, xscale,xlabel},
{xscale,xlabel}=AstraGetScaleAndLabel[xVar];
{yscale,ylabel}=AstraGetScaleAndLabel[yVar];


xlist =xscale AstraArray[xVar, astraparticles];
ylist =yscale  AstraArray[yVar, astraparticles];

ListPlot[Transpose[{xlist, ylist}],Evaluate[FilterRules[{opts}, Options[ListPlot]]],  Joined->False, FrameLabel->{xlabel, ylabel}, PlotStyle->Black, PlotRange->All]

];

BunchSigmaMatrix[bunch_]:=Module[{means, coords},
means = Mean[bunch[[All, 1;;6]]];
coords=#[[1;;6]]-means&/@bunch;
Covariance[coords,coords]
];

Dark Current Tracker Analysis DO NOT EDIT THIS, look in bmad_to_opal

Clear[BmadGetScaleAndLabel];
BmadGetScaleAndLabel[PSvariable_String]:=Module[{},
Switch[PSvariable,
	"x", {(10^3), "x (mm)"},
	"cp_x",{ 1, "c (Subscript[p, x]) (eV)"},
	"y",{(10^3), "y (mm)"},
	"cp_y",{1, "c (Subscript[p, y]) (eV)"},
	"s",{1, "s (m)"},
	"z",{1, "z (m)"},
	"cp_s", { (1/(10^6)), "c (Subscript[p, s]) (MeV)"},
	"cp_z", { (1/(10^6)), "c (Subscript[p, z]) (MeV)"},
	"t",   {(10^9), "t (ns)"},
	"E_x", { (10^(-6)), "(Subscript[E, x]) (MV/m)"},
	"E_y", { (10^(-6)), "(Subscript[E, y]) (MV/m)"},
	"E_s", { (10^(-6)), "(Subscript[E, s]) (MV/m)"},
	"E_z", { (10^(-6)), "(Subscript[E, z]) (MV/m)"},
	"B_x", { (10^3), "(Subscript[B, x]) (mT)"},
	"B_y", { (10^3), "(Subscript[B, y]) (mT)"},
	"B_s", { (10^3), "(Subscript[B, s]) (mT)"},
	"B_z", { (10^3), "(Subscript[B, z]) (mT)"},


	_, Print["Bad variable: "<>PSvariable];Abort[]
]
];

ImportBmadTimeParticles[particlefile_String]:=Module[{dat,ppos, timetracks, cnames, cunits},

Print["Importing Bmad time particles: "<>particlefile];
Check[dat=Import[particlefile, "Table"];, Print["Problem with: "<>particlefile];Abort[];];
cnames=dat[[1]]; (*names of columns*)
cunits=dat[[2]];(*Units of columns*)

Do[ix[cnames[[i]]]=i, {i,1,Length[cnames]}];
Do[units[cnames[[i]]]=cunits[[i]], {i,1,Length[cnames]}];

ppos=Position[dat, "particle"][[All,1]];
timetracks = Table[dat[[ppos[[i]]+1;;ppos[[i+1]]-1 ]], {i,1,Length[ppos]-1}];
timetracks  =timetracks~Join~{dat[[Last[ppos]+1;; ]]}; (*Append final particle*)
Print["Imported "<>ToString[Length[ppos]]<>" particles, "<>ToString[Length[dat]]<>" data points"];
timetracks
]

Clear[BmadPlot];
BmadPlot[xVar_String,yVar_String,timeTracks_, opts:OptionsPattern[]]:=Module[{yscale, ylabel, xscale,xlabel,xlist,  ylist,tracks, dat, fy},
{xscale,xlabel}=BmadGetScaleAndLabel[xVar];
{yscale,ylabel}=BmadGetScaleAndLabel[yVar];

tracks=timeTracks;

dat = Table[{xscale  #[[ix[xVar]]],yscale  #[[ix[yVar]]] }&/@tracks[[i]], {i,1,Length[tracks]}];

ListPlot[dat,
FrameLabel->{xlabel, ylabel},Evaluate[FilterRules[{opts}, Options[ListPlot]]], Joined->True]

];
Options[BmadPlot] =  {PlotRange->All,  Joined->True }

Options[BmadPlot]

Tests

GetATracks[]:=Module[{stepfiles, adat},
stepfiles = Flatten[ Import[aworkdir<>"steps.txt", "Table"]];
adat = Table[{f, ImportAstraParticles[aworkdir<>f]}, {f, stepfiles }];
Flatten[adat[[All, 2]], 1]
]

labels=Inset[
Grid[{
{Style["Astra", Red]},
{Style["Bmad", Blue]}
}, Alignment->Left],Scaled[{.9,.7}]];

apex_gun

aworkdir=aroot<>"test/apex_gun/"

atracks=GetATracks[];
g1=AstraPhaseSpace["z", "x", atracks, PlotStyle->Red]

btracks= ImportBmadTimeParticles[aworkdir<>"lat.out"]; 
g2=BmadPlot["z", "x",btracks, PlotStyle->Blue];
Show[g1,g2,Epilog->labels]

Show[AstraPhaseSpace["z", "pz", atracks, PlotStyle->Red],
 BmadPlot["z", "cp_z",btracks, PlotStyle->Blue]]

adat=ImportAstraParticles[aworkdir<>"test_astra.0035.001"];
g1=AstraPhaseSpace["z","x",adat, PlotStyle->Orange]

sigmat=BunchSigmaMatrix[adat];
(Sqrt[Diagonal[sigmat]])

AstraPhaseSpace["z","pz",adat, AspectRatio->1]

emitx=ImportAstraEmittance[aworkdir<>"test_astra"];
ListPlot[{#[[1]],#[[6]]}&/@emitx]

dcgun

aworkdir=aroot<>"test/dcgun/"

atracks=GetATracks[];
g1=AstraPhaseSpace["z", "x", atracks, PlotStyle->Red]

btracks= ImportBmadTimeParticles[aworkdir<>"lat.out"]; 
g2=BmadPlot["z", "x",btracks, PlotStyle->Blue];
Show[g1,g2,Epilog->labels]

adat=ImportAstraParticles[aworkdir<>"test_astra.0035.001"];
g1=AstraPhaseSpace["z","x",adat, PlotStyle->Orange]

sigmat=BunchSigmaMatrix[adat];
(Sqrt[Diagonal[sigmat]])

AstraPhaseSpace["z","pz",adat, AspectRatio->1]

emitx=ImportAstraEmittance[aworkdir<>"test_astra"];
ListPlot[{#[[1]],#[[6]]}&/@emitx]

dat=Import[aworkdir<>"input_beam.ini", "Table"];

ListPlot[dat[[All,4]]]

ListPlot[{#[[1]], #[[2]]}&/@dat]

apex_gun

aworkdir=aroot<>"test/apex_gun/"

stepfiles = Flatten[ Import[aworkdir<>"steps.txt", "Table"]];
adat = Table[{f, ImportAstraParticles[aworkdir<>f]}, {f, stepfiles }];
atracks = Flatten[adat[[All, 2]], 1];

g1=AstraPhaseSpace["z", "pz", atracks, PlotStyle->Red]

btracks= ImportBmadTimeParticles[aworkdir<>"lat.out"];

g2=BmadPlot["z", "cp_z",btracks, PlotStyle->Blue]

Show[g1,g2,Epilog->labels]

BmadPlot["z", "E_x",btracks, PlotStyle->Blue]

BmadPlot["x", "E_x",btracks, PlotStyle->Blue]

apex_buncher

aworkdir=aroot<>"test/apex_buncher/"

atracks=GetATracks[];
g1=AstraPhaseSpace["t", "pz", atracks, PlotStyle->Red]

btracks= ImportBmadTimeParticles[aworkdir<>"lat.out"];
g2=BmadPlot["t", "cp_z",btracks, PlotStyle->Blue];
Show[g1,g2,Epilog->labels]

BmadPlot["z", "E_x",btracks, PlotStyle->Blue]

apex_buncher

aworkdir=aroot<>"test/cavity9/"

atracks=GetATracks[];
g1=AstraPhaseSpace["t", "pz", atracks, PlotStyle->Red]

btracks= ImportBmadTimeParticles[aworkdir<>"lat.out"];
g2=BmadPlot["t", "cp_z",btracks, PlotStyle->Blue];
Show[g1,g2,Epilog->labels]

BmadPlot["z", "E_x",btracks, PlotStyle->Blue]

fodo

aworkdir=aroot<>"test/fodo/"

atracks=GetATracks[];
g1=AstraPhaseSpace["z", "x", atracks, PlotStyle->Red]

btracks= ImportBmadTimeParticles[aworkdir<>"lat.out"];
g2=BmadPlot["z", "x",btracks, PlotStyle->Blue];
Show[g1,g2,Epilog->labels]

BmadPlot["z", "B_y",btracks, PlotStyle->Blue]

Emit

aemit=ImportAstraEmittance[aworkdir<>"test_astra"];
(*aemitnomirror=ImportAstraEmittance[workdir<>"astra_sc_nomirror"];*)
 PlotAstraEmit[  "emit_x", "z", aemit]

PlotAstraEmit[  "rms_x", "t", aemit]

workdir<>"astra"

ListPlot[AstraEmitArray["rms_x", aemit][[1]]]

.

PlotAstraEmit[  "kinetic_energy", "t", aemit]

gastrasigx = PlotAstraEmit[  "rms_x", "t", aemit]

gastrasigy = PlotAstraEmit[  "rms_y", "t", aemit]

gastraemitx= PlotAstraEmit[  "emit_x", "t", aemit]

gastraemity= PlotAstraEmit[  "emit_y", "t", aemit]

Show[ggptsigx, gopalsigx ,gastrasigx,  Epilog->labels]

8.10316500000000062*750/350

Show[ggptemitx, gopalemitx,gastraemitx, PlotRange->{0,3}, Epilog->labels]

OPAL vs GPT vs Astra

{gptnosc, gptscdefault, gptsccathode} = ImportGPTSummary[workdir<>#]&/@{"summary_nosc.txt", "summary_sc_default.txt", "summary_sc_cathode.txt"};

{ostatnosc, ostat} = ImportOPALStat[workdir<>#]&/@{"opal_nosc.stat", "opal_sc_1.stat"};

{anosc, asc}=ImportAstraEmittance[aroot<>"injector_gun/"<>#]&/@{"astra", "astra_sc"};

labels=Inset[
Grid[{
{Style["GPT no sc", Green]},
{Style["GPT sc default", Red]},
{Style["GPT sc cathode",Black]},
{Style["OPAL", Purple]},
{Style["Astra", Orange]}
}, Alignment->Left],Scaled[{.2,.7}]];

gxVar="time";
gyVar="stdx";
oxVar = "t";
oyVar="rms_x";
ayVar = "rms_x";
Show[
GPTPlot[ gyVar,gxVar,gptnosc, PlotStyle->Green],
GPTPlot[ gyVar,gxVar,gptscdefault, PlotStyle->{Dashed,Red}],
GPTPlot[ gyVar,gxVar,gptsccathode, PlotStyle->Black],
OPALStatPlot[oyVar,oxVar,ostatnosc,  PlotStyle->{Dashed, Purple}],
OPALStatPlot[oyVar,oxVar,ostat,  PlotStyle->Purple],
PlotAstraEmit[ ayVar , "t",anosc, PlotStyle->{Dashed,Orange}],
PlotAstraEmit[ayVar , "t",  asc, Joined->False,PlotStyle->{Orange}],
 Epilog->labels, PlotRange->All]

{"nemixrms", "emit_x"},
{"rms_px", "rms_px"},
{"energy", "energy"},
{"stdz", "rms_s"}

res=Table[
gxVar="time";
gyVar=prop[[1]];
oxVar = "t";
oyVar=prop[[2]];
Show[
GPTPlot[ gyVar,gxVar,gptnosc, PlotStyle->Green],
GPTPlot[ gyVar,gxVar,gptscdefault, PlotStyle->{Dashed,Red}],
GPTPlot[ gyVar,gxVar,gptsccathode, PlotStyle->Black],
OPALStatPlot[oyVar,oxVar,ostatnosc,  PlotStyle->{Dashed, Purple}],
OPALStatPlot[oyVar,oxVar,ostat,  PlotStyle->Purple],
PlotAstraEmit[ oyVar , "t",anosc, PlotStyle->{Dashed,Orange}],
PlotAstraEmit[oyVar , "t",  asc,PlotStyle->{Orange}],
 Epilog->labels,  ImagePadding->60,
PlotRange->{All,Automatic}],

{prop, {
{"nemixrms", "emit_x"},
{"stdx", "rms_x"},
{"kinetic_energy","kinetic_energy"}
} }];
Column[res]