Bmad to GPT test

root=NotebookDirectory[];
scratchdir="~/Scratch/";

(* 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);

Dark Current Tracker Analysis

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, dat},
{xscale,xlabel}=BmadGetScaleAndLabel[xVar];
{yscale,ylabel}=BmadGetScaleAndLabel[yVar];

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

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

];


Clear[BmadPlot3D];
BmadPlot3D[xVar_String,yVar_String,zVar_String,timeTracks_, opts:OptionsPattern[]]:=Module[{yscale, ylabel, xscale,xlabel, zscale, zlabel,tracks, dat},
{xscale,xlabel}=BmadGetScaleAndLabel[xVar];
{yscale,ylabel}=BmadGetScaleAndLabel[yVar];
{zscale,zlabel}=BmadGetScaleAndLabel[yVar];

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

(*ListPointPlot3D[dat,
AxesLabel->{xlabel, ylabel,zlabel},FilterRules[{opts}, Options[ListPointPlot3D]]]*)

Graphics3D[Line/@dat, BoxRatios->{2,1,1}, ViewPoint->{1.4,-2.4,0}]

];

BmadPSList= {"x", "y", "z","cp_x", "cp_y", "cp_z"};
BmadEBList={"E_x", "E_y", "E_z", "B_x", "B_y", "B_z"};

Clear[BmadPlotAll];
BmadPlotAll[ bmadTracks_, xVar_String, yVarList_, opts:OptionsPattern[]]:=Module[{},
Table[BmadPlot[xVar, yVar,  bmadTracks,opts], {yVar, yVarList}]
]

GPT Analysis (taken from gpt_vs_opal.nb)

ImportGPTSummary[summaryfile_String]:=Module[{gptoutput, gptlabelline},

Check[gptoutput=Import[summaryfile, "Table"];, Print["Problem with "<>summaryfile];Abort[]];

gptlabelline=gptoutput[[1]];
Do[GPTcolumn[ gptlabelline[[i]]] = i, {i,1,Length[gptlabelline]}];
(*gptdat=gptoutput[[2;;-2]];*)

(*Show sample*)
Print[gptoutput[[1;;1]]];
gptoutput[[2;;-2]]

];

Clear[GPTGetScaleAndLabel];
GPTGetScaleAndLabel[PSvariable_String]:=Module[{},
Switch[PSvariable,
(*For particles.txt *)
	"x", {(10^3), "x (mm)"}, 
	"px",{ mec2N, "c (Subscript[p, x]) (eV)"},
	"y",{(10^3), "y (mm)"},
	"py",{ mec2N, "c (Subscript[p, y]) (eV)"},
	"z",{1, "z (m)"},
	"pz", { (1/(10^6))mec2N, "c (Subscript[p, z]) (MeV)"},
	"r", {(10^3), "r (mm)"}, 
	"pr",{ mec2N, "c (Subscript[p, r]) (eV)"},
	"Ex", { (10^(-6)), "(Subscript[E, x]) (MV/m)"},
	"Ey", { (10^(-6)), "(Subscript[E, y]) (MV/m)"},
	"Ez", { (10^(-6)), "(Subscript[E, z]) (MV/m)"},
	"Bx", { (10^3), "(Subscript[B, x]) (mT)"},
	"By", { (10^3), "(Subscript[B, y]) (mT)"},
	"Bz", { (10^3), "(Subscript[B, z]) (mT)"},
	"Er", { (10^(-6)), "|(Subscript[E, r])| (MV/m)"},
	"Br", { (10^3), "|(Subscript[B, r])| (mT)"},
(*For summary.txt *)
	"t",   {(10^9), "t (ns)"},
	"time",   {(10^9), "t (ns)"},
	"Q",   {-(10^12), "-Q (pC)"},
	"avgx",   {(10^6), "<x> (μm)"},
	"avgy",   {(10^6), "<y> (μm)"},
	"avgz",   {1, "<z> (m)"},
	"stdx",   {(10^6), "(Subscript[σ, x]) (μm)"},
	"stdy",   {(10^6), "(Subscript[σ, y]) (μm)"},
	"stdz",   {(10^6), "(Subscript[σ, z]) (μm)"},
	"stdBx",   {1, "(Subscript[σ, βx]) (1)"},
	"stdBy",   {1, "(Subscript[σ, βx]) (1)"},
	"stdBz",   {1, "(Subscript[σ, βx]) (1)"},
	"stdG",   {1, "(Subscript[σ, γ]) (1)"},
	"rms_px",   {mec2N, "c*<γ>(Subscript[σ, βx]) (eV)"},
	"rms_py",   {mec2N, "c*<γ>(Subscript[σ, βy]) (eV)"},
	"rms_ps",   {mec2N, "c*<γ>(Subscript[σ, βz]) (eV)"},
	"nemixrms", {(10^6), "\!\(\*SubscriptBox[\(ϵ\), \(N, \\ x\)]\) (mm-mrad)"},
	"nemiyrms", {(10^6), "\!\(\*SubscriptBox[\(ϵ\), \(N, \\ y\)]\) (mm-mrad)"},
	"nemizrms", {(10^6), "nemizrms (?)"},
	"energy",   {(10^(-6)), "Energy (MeV)"},
	"kinetic_energy",   {(10^(-6)), "Kinetic Energy (MeV)"},
	_, Print["Bad GPT variable: "<>PSvariable];Abort[]
]
];


GPTTimeShift=0;
Clear[GPTArray];
(*Get a single summary column, or a function of columns*)
GPTArray[xVar_String, gptdat_]:=Module[{f},
GPTGetScaleAndLabel[xVar]; (*just used to check*)
f = Switch[xVar,
"energy", mec2N #[[GPTcolumn["avgG"]]]&,
"kinetic_energy", mec2N (#[[GPTcolumn["avgG"]]]-1)&,
"rms_px", #[[GPTcolumn["stdBx"]]]*#[[GPTcolumn["avgG"]]]&,
"rms_py", #[[GPTcolumn["stdBy"]]]*#[[GPTcolumn["avgG"]]]&,
"rms_ps", #[[GPTcolumn["stdG"]]]&,
"t",GPTTimeShift+ #[[GPTcolumn["time"]]]&,
_,  #[[GPTcolumn[xVar]]]& (*defaut*)
];
f/@gptdat

];


Clear[GPTSliceArray];
(*Makes phase space array from slice*)
GPTSliceArray[PSvariable_String, gptslice_]:=Module[{f},
GPTGetScaleAndLabel[PSvariable]; (*just used to check*)
f = Switch[PSvariable,
"x",#[[GPTslicecolumn["x"]]]&,
"y",#[[GPTslicecolumn["y"]]]&,
"z",#[[GPTslicecolumn["z"]]]&,
"r",Sqrt[ ((#[[GPTslicecolumn["x"]]])^2)+((#[[GPTslicecolumn["y"]]])^2)]&,
"pr",#[[GPTslicecolumn["G"]]]*Sqrt[ ((#[[GPTslicecolumn["Bx"]]])^2)+((#[[GPTslicecolumn["By"]]])^2)]&,
"px",#[[GPTslicecolumn["G"]]]*#[[GPTslicecolumn["Bx"]]]&,
"py",#[[GPTslicecolumn["G"]]]*#[[GPTslicecolumn["By"]]]&,
"pz",#[[GPTslicecolumn["G"]]]*#[[GPTslicecolumn["Bz"]]]&,
"Ex",  #[[GPTslicecolumn["fEx"]]]&,
"Ey",  #[[GPTslicecolumn["fEy"]]]&,
"Ez",  #[[GPTslicecolumn["fEz"]]]&,
"Bx",  #[[GPTslicecolumn["fBx"]]]&,
"By",  #[[GPTslicecolumn["fBy"]]]&,
"Bz",  #[[GPTslicecolumn["fBz"]]]&,
"Er",Sqrt[ ((#[[GPTslicecolumn["fEx"]]])^2)+((#[[GPTslicecolumn["fEy"]]])^2)]&,
"Br",Sqrt[ ((#[[GPTslicecolumn["fBx"]]])^2)+((#[[GPTslicecolumn["fBy"]]])^2)]&,
_,   Print["Bad variable: "<>PSvariable];Abort[]
];
f/@ gptslice

];


Clear[GPTPhaseSpace];
GPTPhaseSpace[xVar_,yVar_, gptslice_, opts:OptionsPattern[]]:= Module[{xlist, ylist, yscale, ylabel, xscale,xlabel},
{xscale,xlabel}=GPTGetScaleAndLabel[xVar];
{yscale,ylabel}=GPTGetScaleAndLabel[yVar];


xlist =xscale GPTSliceArray[xVar, gptslice];
ylist =yscale  GPTSliceArray[yVar, gptslice];

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

];

Clear[GPTPlot];
GPTPlot[yVar_String, xVar_String,gptdat_, opts:OptionsPattern[]]:=Module[{yscale, ylabel, xscale,xlabel,xlist,  ylist, datx, daty},
{xscale,xlabel}=GPTGetScaleAndLabel[xVar];
{yscale,ylabel}=GPTGetScaleAndLabel[yVar];

datx=GPTArray[xVar, gptdat];
daty=GPTArray[yVar, gptdat];

ListPlot[Transpose[{xscale datx, yscale daty}],
FrameLabel->{xlabel, ylabel},Evaluate[FilterRules[{opts}, Options[ListPlot]]], Joined->True]


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

Clear[GPTParticlesToBmadAndOPAL];

GPTParticlesToBmadAndOPAL[gptparticlefile_,screennumber_Integer, bmadoutfile_, opaloutfile_]:=
(*Makes .time_particles and .opal_particles style files from a gpt particles.txt style file*)
Module[{screen,nextscreen,nparticles, dat, header,hcol,  pdat, γ,x, βx, y, βy, z, βz, t, bmadline, opalline, time},
Check[dat = Import[gptparticlefile, "Table"];, Print["Problem with : "<>gptparticlefile];Abort[]];
screen=Position[dat, "time"][[screennumber,1]];
nextscreen=Position[dat, "time"][[screennumber+1,1]];

time=dat[[screen]][[2]];

header=dat[[screen+1]];

pdat = dat[[screen+2;;nextscreen-2]];
nparticles=Length[pdat ];
Print["Found "<>ToString[nparticles]<>" particles"];
Print["Time: "<>ToString[FortranForm[time]]];


Do[hcol[header[[i]]] = i, {i,1,Length[header]}];
(*"x","y","z","rxy","Bx","By","Bz","G","t","m","q","nmacro","rmacro","ID"*)
(*Convert to t, x, cpx, y, cpy, z, cpz*)
bmadline[gptline_]:=(x=gptline[[hcol["x"]]];
 βx=gptline[[hcol["Bx"]]]; 
y=gptline[[hcol["y"]]] ;
βy=gptline[[hcol["By"]]];
  z=gptline[[hcol["z"]]] ;
βz=gptline[[hcol["Bz"]]];
γ=gptline[[hcol["G"]]];
{ x, mec2N βx γ, y, mec2N βy γ,z, mec2N βz γ, time, 0}); 
opalline[gptline_]:=(x=gptline[[hcol["x"]]];
 βx=gptline[[hcol["Bx"]]]; 
y=gptline[[hcol["y"]]] ;
βy=gptline[[hcol["By"]]];
  z=gptline[[hcol["z"]]] ;
βz=gptline[[hcol["Bz"]]];
γ=gptline[[hcol["G"]]];
{ x, γ βx ,  y, γ βy ,z, γ βz}); 


Export[bmadoutfile, {{nparticles}}~Join~(bmadline/@pdat), "Table"];
Export[opaloutfile, {{nparticles}}~Join~(opalline/@pdat), "Table"];

];

Trajectory Utilities

ImportGPTTrajectories[trajectoryfile_, tshift_:0]:=Module[{dat,IDpos,pdat, npart,γ,x, βx, y, βy, z, βz, t, temp },
Check[dat=Import[trajectoryfile, "Table"], Print["Problem with: "<>trajectoryfile];Abort[]];
Print[dat[[1;;2]]];
IDpos=Position[dat, "ID"][[All,1]];
npart = Length[IDpos];
pdat=Table[dat[[IDpos[[i]]+3;;IDpos[[i+1]]-2]], {i,1,npart-1}]~Join~{dat[[Last[IDpos]+3;;-2]]};

(*Convert to t, x, cpx, y, cpy, z, cpz*)
Table[
(x=#[[1]]; βx=#[[2]]; y=#[[3]] ;βy=#[[4]];  z=#[[5]] ;βz=#[[6]];γ=#[[7]];t=#[[8]];

Flatten[{ x, mec2N βx γ, y, mec2N βy γ,z, mec2N βz γ, t+ tshift, #[[9;;14]]}])&/@pdat[[i]], {i,1,Length[pdat]}]
]

Particle Utilities

Clear[ImportGPTParticles]
ImportGPTParticles[particlefile_String, type_:"time"]:=Module[{dat,ppos, cnames, gptslices},

Print["Importing: "<>particlefile];
Check[dat=Import[particlefile, "Table"];, Print["Problem with: "<>particlefile];Abort[];];


ppos=Position[dat, type][[All,1]];
gptslicelabels = dat[[ppos[[1]]+1]];

Do[GPTslicecolumn[gptslicelabels[[i]]]=i, {i,1,Length[gptslicelabels]}];

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

];


Clear[GPTPlotSlice];
GPTPlotSlice[ix_Integer, c1_String, c2_String,opts:OptionsPattern[]]:=Module[{slice, scale1, scale2, slicename,particles, avgz, avgx,maxx,varx, avgx2, avgBx, avgBx2},
slice=gptslices[[ix]];
slicename=slice[[1]];
particles=slice[[3;;-2]];

scale1=OptionValue[scalex];
scale2=OptionValue[scaley];

avgz=Mean[#[[ GPTslicecolumn["z"] ]] &/@ particles ];
maxx=Max[#[[ GPTslicecolumn["x"] ]] &/@ particles ];
avgx = Mean[#[[ GPTslicecolumn["x"] ]] &/@ particles ];
varx= Variance[(#[[ GPTslicecolumn["x"] ]]) &/@ particles ];
avgx2 = Mean[(#[[ GPTslicecolumn["x"] ]])^2 &/@ particles ];
avgBx = Mean[#[[ GPTslicecolumn["Bx"] ]] &/@ particles ];
avgBx2 = Mean[(#[[ GPTslicecolumn["Bx"] ]])^2 &/@ particles ];

Grid[{
{ListPlot[  {scale1 *#[[ GPTslicecolumn[c1] ]], scale2*#[[ GPTslicecolumn[c2] ]]}&/@ particles, ImageSize->500,Frame->True, PlotRange->All,PlotStyle->Black,BaseStyle->{FontSize->15},FrameLabel->{scale1 c1,scale2 c2},  AspectRatio->1]},
{"average z = ", avgz},
{"max x = ", maxx},
{"average x = " ,avgx},
{"RMS x = " ,(Sqrt[avgx2])},
{"sqr(Var(x)) = " ,(Sqrt[varx])},
{"average Bx = " ,avgBx},
{"\!\(\*SubscriptBox[\(ϵ\), \(n, \\ x\)]\) = ", (Sqrt[avgx2 avgBx2])}
},Frame->True]

];
Options[GPTPlotSlice] = 
{scalex->1, scaley->1}~Join~Options[ListPlot]~Join~{PlotRange->All};

Clear[GPTPlotSlice3D];
GPTPlotSlice3D[ix_Integer, c1_String, c2_String, c3_String, opts:OptionsPattern[]]:=Module[{slice, scale1, scale2, slicename,dat, scale3, particles, avgx,maxx,varx, avgx2, avgBx, avgBx2},
slice=gptslices[[ix]];
slicename=slice[[1]];
particles=slice[[3;;-2]];

scale1=OptionValue[scalex];
scale2=OptionValue[scaley];
scale3=OptionValue[scalez];



 dat={scale1 *#[[ GPTslicecolumn[c1] ]],
 scale2*#[[ GPTslicecolumn[c2] ]],
 scale3*#[[ GPTslicecolumn[c3] ]]
}&/@ particles;

Grid[{
{ix},
{ListPointPlot3D[dat , ImageSize->500, PlotStyle->Black,BaseStyle->{FontSize->15},AxesLabel->{scale1 c1,scale2 c2, scale3 c3},  AspectRatio->1,
Evaluate[FilterRules[{opts}, Options[ListPointPlot3D]]] ]
}

},Frame->True]

]
Options[GPTPlotSlice3D] = 
{scalex->1, scaley->1, scalez->1}~Join~Options[ListPointPlot3D]~Join~{PlotRange->All};


Clear[GPTPhaseSpaceSlice];
GPTPhaseSpaceSlice[ix_Integer]:=Module[{slice, slicename, particles,dat},
slice=gptslices[[ix]];
slicename=slice[[1]];
particles=slice[[3;;-2]];
 dat={
#[[ GPTslicecolumn["x"] ]],
#[[ GPTslicecolumn["Bx"] ]],
#[[ GPTslicecolumn["y"] ]],
#[[ GPTslicecolumn["By"] ]],
#[[ GPTslicecolumn["z"] ]],
#[[ GPTslicecolumn["Bz"] ]]
}&/@ particles;
dat
];


GPTSigmaMatrixTerm[i_Integer,j_Integer,gptbunch_]:=Module[{charge,wC,averagei,averagej,nparticles},
averagei=Mean[#[[i]]&/@gptbunch];
averagej=Mean[#[[j]] &/@gptbunch];
(*Subtract off Centroid for Covariance*)
Mean[(#[[i]]-averagei) (#[[j]]-averagej)&/@gptbunch] 
]

GPTBunchStats[gptbunch_]:=Module[{ averages, sigmamatrix, ϵNx},
averages=Table[Mean[#[[i]]&/@gptbunch], {i,1,6}];

sigmamatrix=Table[ GPTSigmaMatrixTerm[i, j, gptbunch], {i,1,6}, {j,1,6}];

{averages, sigmamatrix}

];
DisplayGPTBunchStats[gptbunch_]:=Module[{res, averages, sigmamatrix, avgγ, θ},

res=GPTBunchStats[gptbunch];
averages=res[[1]];
sigmamatrix=res[[2]]; 

θ = ArcTan[((averages[[2]])/(averages[[6]]))];

avgγ=1/(Sqrt[1-(averages[[2]]^2+averages[[3]]^2+averages[[6]]^2)]);

Grid[{{"avgx", "avgBx", "avgy", "avgBy", "avgz", "avgBz"},
averages,
{},
{"stdx", "stdBx", "stdBy", "stdBy", "stdz", "stdBz"},
(Sqrt[Diagonal[sigmamatrix]]), 
{},
{"ϵNx", ," ϵNy", ,"(Subscript[θ, xz])(deg)","<γ>"  },
{avgγ Sqrt[Det[sigmamatrix[[1;;2,1;;2]]]], ,avgγ Sqrt[Det[sigmamatrix[[3;;4,3;;4]]]], ,θ *(180/π),avgγ  }
}
,Frame->True]


]



GPTCenterBunch[gptbunch_]:=Module[{avgx, avgy, avgz},
avgx=Mean[#[[1]]&/@gptbunch];
avgy=Mean[#[[3]]&/@gptbunch];
avgz=Mean[#[[5]]&/@gptbunch];

{
#[[1]] - avgx,
#[[2]],
#[[3]] - avgy,
#[[4]], 
#[[5]] - avgz,
#[[6]]
}&/@gptbunch
]

GPTRealignXZ[gptbunch_]:=Module[{newgptbunch, avgBx, avgBz,avgx, avgz, xC=1, yC=3, zC=5, BxC=2, ByC=4, BzC=6, θ},

newgptbunch = GPTCenterBunch[gptbunch];

avgBx=Mean[#[[BxC]]&/@newgptbunch ];

avgBz=Mean[#[[BzC]]&/@newgptbunch ];
θ=ArcTan[(avgBx/avgBz)];



{#[[xC]] Cos[θ] - #[[zC]]Sin[θ],
#[[BxC]] Cos[θ] - #[[BzC]]Sin[θ], 
#[[yC]],
 #[[ByC]],
#[[xC]] Sin[θ] + #[[zC]]Cos[θ],
#[[BxC]] Sin[θ] + #[[BzC]]Cos[θ]
} &/@newgptbunch

]

Line test

workdir = root<>"dc_gun/"

Zcut = 10;

Bmad

bmadTracksFull=ImportBmadTimeParticles[workdir<>"lat.out"];

bmadTracks=Table[Select[bmadTracksFull[[i]], #[[ix["z"]]]≤ Zcut&], {i,1,Length[bmadTracksFull]}];


xVar="z";
gBmad=BmadPlotAll[bmadTracks, xVar, Join[BmadPSList,BmadEBList]] ;
gBmad

BmadPlot["z","E_z",bmadTracks, PlotStyle->{Red, PointSize[Tiny]}]

gref = BmadPlot["t","cp_z",bmadTracks, PlotStyle->{Red, PointSize[Tiny]}]

GPT

tshift =0;
gptTracksFull=ImportGPTTrajectories[workdir<> "gpt_trajectories.txt", tshift];
gptTracks=Table[Select[gptTracksFull[[i]], #[[ix["z"]]]≤ Zcut&], {i,1,Length[gptTracksFull]}];
Show[BmadPlot["z","cp_z", gptTracks, PlotStyle->{Red,PointSize[Tiny]}]]

ddat=Import[workdir<>"1D_fieldmap_1_ASCII.gpt", "Table"][[2;;]];
fEz = Interpolation[ddat];
Plot[fEz[z], {z,Min[ddat[[All,1]]], Max[ddat[[All,1]]]}, PlotRange->All]

volt = NIntegrate[fEz[z]Exp[I 2 π 1.3 (10^9)/299792458 z], {z,0, 0.5}]

Arg[volt]

fEz[z]

Table[BmadPlot["z",yVar, gptTracks, PlotStyle->{Red, PointSize[Tiny]}], {yVar,  BmadPSList}]

Comparison

ggpt=Table[BmadPlot["z",yVar, gptTracks, PlotStyle->{Green,PointSize[Tiny]}], {yVar,  BmadPSList}];
gbmad=Table[BmadPlot["z",yVar, bmadTracks, PlotStyle->{Blue, PointSize[Tiny]}], {yVar,  BmadPSList}];
gcomp = Table[Show[ gbmad[[i]],ggpt[[i]]], {i,1,Length[ggpt]}]

Export[workdir<>"comparison3D.pdf", gcomp]

ll

BmadPlot3D["z", "x", "y",gptTracks]

Show[BmadPlot3D["z", "x", "y", bmadTracks], BmadPlot3D["z", "x", "y",gptTracks]]

Scratch

gptsummary = ImportGPTSummary[workdir<>"gpt_summary.txt"];

GPTPlot[ "avgz","time",gptsummary, PlotStyle->Green]

ggptenergy = GPTPlot[ "energy","time",gptsummary, PlotStyle->Green]

ggptsigx = GPTPlot[ "stdx","time",gptsummary, PlotStyle->Green]

ggptsigy = GPTPlot[ "stdy","time",gptsummary, PlotStyle->Green]

ggptsigpx =GPTPlot[ "rms_px","time",gptsummary, PlotStyle->Green]

ggptsigpy =GPTPlot[ "rms_py","time",gptsummary, PlotStyle->Green]

ggptsigz= GPTPlot[ "stdz","time",gptsummary, PlotStyle->Green]

ggptemitx= GPTPlot[ "nemixrms","time",gptsummary, PlotStyle->Green]

ggptemity= GPTPlot[ "nemiyrms","time",gptsummary, PlotStyle->Green]

GPTPlot[ "nemizrms","time",gptsummary, PlotStyle->Green]