Dark Current Analysis

root=NotebookDirectory[]

groot="~/Desktop/"
scratchdir="~/Scratch/";

Basic Definitions

(* Define some constants *)

numHeaderLines=1 ;  (* 1 lines of header data (particle)*)
mec2N=510998.928; 
cLightN=299792458.;
frequency = 1.3 (10^9)(*Hz*);
(* Rest mass of an electron *)

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

(* Set indexed variables for file data and fileIn data *)
ix["x"]=1;
ix["cpx"]=2;
ix["y"]=3;
ix["cpy"]=4;
ix["s"]=5;
ix["cps"]=6;
ix["t"]=7;
ix["charge"]=8;
ix["hitθ"]=9;
ix["s_origin"]=10;
ix["ix_ele"]=11;

ix["Ex"] = 10;
ix["Ey"] = 11;
ix["Es"] = 12;
ix["Bx"] = 13;
ix["By"] = 14;
ix["Bs"] = 15;

Emagnitudef=Sqrt[#[[ix["Ex"]]]^2+#[[ix["Ey"]]]^2+#[[ix["Es"]]]^2]&;

radiusf = Sqrt[((#[[ix["x"]]])^2)+((#[[ix["y"]]])^2)]&;
chargef = #[[ix["charge"]]]&;



pctotf=Sqrt[(#[[ix["cpx"]]]^2+#[[ix["cpy"]]]^2+#[[ix["cps"]]]^2)]&;
energyf=Sqrt[mec2N^2+(#[[ix["cpx"]]]^2+#[[ix["cpy"]]]^2+#[[ix["cps"]]]^2)]&;

TrackDomain[trackdat_, colname_String:"s"]:={ Min[Table[Min[t[[All, ix[colname]]]], {t,trackdat}]], Max[Table[Max[t[[All, ix[colname]]]], {t, trackdat}]]}

ImportTracks[particlefile_]:=Module[{begin, end, trackDat},

Check[trackDat=Import[particlefile,"Table"];,Abort[]];
(* Find where each data set begins and ends *)
begin=Position[trackDat, "particle"][[All,1]]+1;
end=begin[[#+1]]-2&/@Range[Length[begin]-1];
end=Append[end,Length[trackDat]];

(* Make an array of all data sets and one of phases *)
trackDat[[begin[[#]];;end[[#]]]]&/@Range[Length[begin]]
]

Clear[WallPlot];
WallPlot[wallout_String, opts:OptionsPattern[]]:=Module[{wdat, cnames, wallix, innerwalls, outerwalls, angleixmax,topwalls, bottomwalls},
(* Import Data *)
Check[wdat=Import[wallout, "Table"];,Abort[]];

(*Get Column names*)
cnames=wdat[[2]];
Do[wallix[cnames[[i]]] = i, {i,1,Length[cnames]}];
(*
They should be:
x normal_x y normal_y z normal_z ix_ele angle_index *)
wdat=wdat[[4;;]];
angleixmax=Max[wdat[[All, wallix["angle_index"]]]];

(*Split by element index*)
wdat = SplitBy[wdat, #[[wallix["ix_ele"]]]&];

(* top view, or top and side view *)
Switch[ angleixmax ,
1,
(*Get inner and outer walls by selecting the angle_index. This is so that the points can be joined properly in the plot*)
innerwalls=Table[{ #[[wallix["z"] ]], #[[wallix["x"] ]]}&/@Select[wdat[[i]], #[[wallix["angle_index"]]]==0&], {i,1,Length[wdat]}];
outerwalls=Table[{ #[[wallix["z"] ]], #[[wallix["x"] ]]}&/@Select[wdat[[i]], #[[wallix["angle_index"]]]==1&], {i,1,Length[wdat]}];


ListPlot[{Flatten[innerwalls,1], Flatten[outerwalls,1]},Evaluate[FilterRules[{opts}, Options[ListPlot]]],Joined->True,PlotStyle->Black, AspectRatio->Automatic],
3,
innerwalls=Table[{ #[[wallix["z"] ]], #[[wallix["x"] ]]}&/@Select[wdat[[i]], #[[wallix["angle_index"]]]==0&], {i,1,Length[wdat]}];
outerwalls=Table[{ #[[wallix["z"] ]], #[[wallix["x"] ]]}&/@Select[wdat[[i]], #[[wallix["angle_index"]]]==2&], {i,1,Length[wdat]}];
topwalls=Table[{ #[[wallix["z"] ]], #[[wallix["y"] ]]}&/@Select[wdat[[i]], #[[wallix["angle_index"]]]==1&], {i,1,Length[wdat]}];
bottomwalls=Table[{ #[[wallix["z"] ]], #[[wallix["y"] ]]}&/@Select[wdat[[i]], #[[wallix["angle_index"]]]==3&], {i,1,Length[wdat]}];

{ListPlot[{Flatten[innerwalls,1], Flatten[outerwalls,1]},Evaluate[FilterRules[{opts}, Options[ListPlot]]],Joined->True,PlotStyle->Black, AspectRatio->Automatic],
ListPlot[{Flatten[topwalls,1], Flatten[bottomwalls,1]},Evaluate[FilterRules[{opts}, Options[ListPlot]]], Joined->True,PlotStyle->Black, AspectRatio->Automatic]},
_,
Print["error: can't draw this many angles: "<>ToString[angleixmax+1]];Abort[]
]

];

Clear[WallPlot3D];
WallPlot3D[wallout_String, plotrange_:Automatic, color_:Black]:=Module[{wdat, cnames, wallix, innerwalls, outerwalls, angleix},
(* Import Data *)
wdat=Import[wallout, "Table"];
(*Get Column names*)
cnames=wdat[[2]];
Do[wallix[cnames[[i]]] = i, {i,1,Length[cnames]}];

angleix = wdat[[4;;,wallix["angle_index"] ]];
angleix = Union[angleix,angleix];

(*wdat = Select[wdat[[4;;]],869< #[[wallix["ix_ele"]]]<911&];*)
wdat=wdat[[4;;]];
(*Split by s*)
wdat = SplitBy[wdat, #[[wallix["s"]]]&];

(*
They should be:
x normal_x y normal_y z normal_z ix_ele angle_index s*)

Graphics3D[Flatten[{color}]~Join~Table[Line[{ #[[wallix["z"] ]], #[[wallix["x"] ]],  #[[wallix["y"] ]]}&/@Join[w, {w[[1]]}]], {w,wdat}],  BoxRatios->{2,Automatic,Automatic}, PlotRange->plotrange]
];

WallNormalPlot[wallout_String, vecscale_]:=Module[{wallDat, myWallDat1, myWallDat2, gWall1, gWall2, arrows, gArrows},
(* Import Data *)
Check[wallDat=Import[wallout,"Table"][[4;;]];, Print["Problem with: "<>wallout]Abort[]];
(* Make an array of last point of all data sets *)
arrows=Arrow[ {{#[[ix["s"]]], #[[ix["x"]]]}, {#[[ix["s"]]]+vecscale#[[ix["cps"]]] , #[[ix["x"]]]+vecscale#[[ix["cpx"]]] }}] &/@wallDat;
gArrows=Graphics[{Arrowheads[Small], arrows}]
]

EnergyTrackPlot[myTrackDat_,opts:OptionsPattern[]]:=Module[{sortedTrackDat,origins,gorigins,energies, tableTrack, gTrackEnergy, Δenergies},

origins=myTrackDat[[All,1]];
gorigins=ListPlot[{#[[ix["s"]]],#[[ix["x"]]]}&/@origins, PlotStyle->Black];

sortedTrackDat=Sort[myTrackDat,energyf[Last[#1]]<energyf[Last[#2]]&];
energies=Table[energyf/@sortedTrackDat[[i]],{i,1,Length[sortedTrackDat]}];
Δenergies=(Max[energies]-Min[energies]);
energies =If[Δenergies >1.0(*eV*), 1-(energies-Min[energies])/Δenergies, 1];

tableTrack=Table[{#[[ix["s"]]],#[[ix["x"]]]}&/@sortedTrackDat[[i]],{i,1,Length[sortedTrackDat]}];
gTrackEnergy=ListPlot[tableTrack,Evaluate[FilterRules[{opts}, Options[ListPlot]]],FrameLabel->{"s (m)","x (cm)"},Joined->True,PlotStyle->Array[ColorData["BlueGreenYellow",Last[energies[[#]]]]&,Length[energies]],PlotRange->{All,All}];
Show[gTrackEnergy,gorigins,gWall,FrameTicks->{Automatic,Table[{i,100 i}, {i,-.1,.1,0.05}],None,None}, Axes->False]
];
Clear[ChargeTrackPlot];
ChargeTrackPlot[myTrackDat_,opts:OptionsPattern[]]:=Module[{sortedTrackDat,origins,gorigins, tableTrack, gTrackEnergy,charges, plane, wall},
plane=OptionValue[Plane];
origins=myTrackDat[[All,1]];
gorigins=ListPlot[{#[[ix["s"]]],#[[ix[plane]]]}&/@origins, PlotStyle->Black];

sortedTrackDat=Sort[myTrackDat,chargef[Last[#1]]<chargef[Last[#2]]&];
charges=Table[chargef/@sortedTrackDat[[i]],{i,1,Length[sortedTrackDat]}];
charges = charges/Max[charges];

tableTrack=Table[{#[[ix["s"]]],#[[ix[plane]]]}&/@sortedTrackDat[[i]],{i,1,Length[sortedTrackDat]}];
gTrackEnergy=ListPlot[tableTrack,Evaluate[FilterRules[{opts}, Options[ListPlot]]],FrameLabel->{"s (m)",plane<>" (cm)"},Joined->True,PlotStyle->Array[ColorData["DarkRainbow",Last[charges[[#]]]]&,Length[charges]],PlotRange->{All,All}];

If[plane=="x", wall=gWall];
If[plane=="y", wall=gWall2];
Show[gTrackEnergy,gorigins,wall,FrameTicks->{Automatic,Table[{i,100 i}, {i,-.1,.1,0.05}],None,None}, Axes->False]
];
Options[ChargeTrackPlot] = {Plane->"x"}~Join~Options[ListPlot];

EnergyTrackPlotRadius[myTrackDat_,opts:OptionsPattern[]]:=Module[{sortedTrackDat,origins,gorigins,energies, tableTrack, gTrackEnergy, Δenergies},

origins=myTrackDat[[All,1]];
gorigins=ListPlot[{#[[ix["s"]]],Sqrt[#[[ix["x"]]]^2 +#[[ix["y"]]]^2]}&/@origins, PlotStyle->Black];

sortedTrackDat=Sort[myTrackDat,energyf[Last[#1]]<energyf[Last[#2]]&];
energies=Table[energyf/@sortedTrackDat[[i]],{i,1,Length[sortedTrackDat]}];
Δenergies=(Max[energies]-Min[energies]);
energies =If[Δenergies >1.0(*eV*), 1-(energies-Min[energies])/Δenergies, 1];

tableTrack=Table[{#[[ix["s"]]],Sqrt[#[[ix["x"]]]^2 +#[[ix["y"]]]^2]}&/@sortedTrackDat[[i]],{i,1,Length[sortedTrackDat]}];
gTrackEnergy=ListPlot[tableTrack,Evaluate[FilterRules[{opts}, Options[ListPlot]]],FrameLabel->{"s (m)","r (cm)"},Joined->True,PlotStyle->Array[ColorData["BlueGreenYellow",Last[energies[[#]]]]&,Length[energies]],PlotRange->{All,All}];
Show[gTrackEnergy,gorigins,gWall,FrameTicks->{Automatic,Table[{i,100 i}, {i,-.1,.1,0.05}],None,None}, Axes->False]
];

ChargeTrackPlotRadius[myTrackDat_,opts:OptionsPattern[]]:=Module[{sortedTrackDat,origins,gorigins,charges, tableTrack, gTrackEnergy},

origins=myTrackDat[[All,1]];
gorigins=ListPlot[{#[[ix["s"]]],Sqrt[#[[ix["x"]]]^2 +#[[ix["y"]]]^2]}&/@origins, PlotStyle->Black];

sortedTrackDat=Sort[myTrackDat,chargef[Last[#1]]<chargef[Last[#2]]&];
charges=Table[chargef/@sortedTrackDat[[i]],{i,1,Length[sortedTrackDat]}];
charges = charges/Max[charges];

tableTrack=Table[{#[[ix["s"]]],Sqrt[#[[ix["x"]]]^2 +#[[ix["y"]]]^2]}&/@sortedTrackDat[[i]],{i,1,Length[sortedTrackDat]}];
gTrackEnergy=ListPlot[tableTrack,Evaluate[FilterRules[{opts}, Options[ListPlot]]],FrameLabel->{"s (m)","r (cm)"},Joined->True,PlotStyle->Array[ColorData["Rainbow",Last[charges[[#]]]]&,Length[charges]],PlotRange->{All,All}];
Show[gTrackEnergy,gorigins,gWall,FrameTicks->{Automatic,Table[{i,100 i}, {i,-.1,.1,0.05}],None,None}, Axes->False]
];

(*TODO*)
FowlerNordheimCurrent[EdotN_, βFN_, AFN_]:=Module[{},

(3.85 (10^(-7)))AFN (((βFN*EdotN))^2) Exp[ -5.464 (10^10) / (βFN*EdotN) ] 
];


CheckFieldEmissionCharge[singletrack_]:=Module[{X, r,charge,chargeFN, Efield},
X = singletrack[[1]];
charge=X[[ix["charge"]]];
r = (Sqrt[((X[[ix["x"]]])^2)+((X[[ix["y"]]])^2)]);
Efield = Emagnitudef[X];
 chargeFN=FowlerNordheimCurrent[Efield, 100.(*βFN*),2 π r ((1.57310318452660/100)(*(Lwall/numSources)*))((1000 (*num emitters/m2*))/(50(*phases*)))(1/(1.3 (10^9))) 1. (10^(-15))] ;
{X[[ix["s"]]],0.1charge/chargeFN}

];

(* 2D Force vectors for visualization*)
EForceArrow[particle_, size_, cutoff_]:=Module[{s1=size*particle[[ix["Es"]]], x1= size*particle[[ix["Ex"]]]},
If[Sqrt[(s1^2)+(x1^2)]>cutoff,
Arrow[{{particle[[ix["s"]]], particle[[ix["x"]]]},
	 {particle[[ix["s"]]] +  s1,particle[[ix["x"]]] +  x1} }],
{}]
];


BForceArrow[particle_, size_, cutoff_]:=Module[{Etot, vx, vs,x1,s1},
Etot=energyf[particle];
vx =cLightN*particle [[ix["cpx"]]]/Etot;
vs =cLightN*particle [[ix["cps"]]]/Etot;

s1 = size*vx*particle[[ix["By"]]];
x1=- size*vs*particle[[ix["By"]]];
If[Sqrt[(s1^2)+(x1^2)]>cutoff,
Arrow[{{particle[[ix["s"]]], particle[[ix["x"]]]},
	 {particle[[ix["s"]]] + s1,particle[[ix["x"]]] +x1} }]
,{}]

];

EBVectors[track_, norm_, cutoff_]:=Module[{Efields, Bfields},

Efields={Red}~Join~(EForceArrow[#,norm ,cutoff]&/@track);
Bfields={Blue}~Join~(BForceArrow[#,norm, cutoff ]&/@track);
Graphics[Join[{Arrowheads[Tiny]},Bfields,Efields ]]
]

Histograms Procedures SHOULD BE REWRITTEN

LoadEndParticles[outfile_]:=Module[{endDat, begin,end, myEndDat},
(* Import Data *)
endDat=Import[outfile,"Table"][[3;;]];

(* Find where each data set begins and ends *)
begin=Position[endDat,"particle"][[All,1]]+numHeaderLines;
end=begin[[#+1]]-numHeaderLines-1&/@Range[Length[begin]-1];
end=Append[end,Length[endDat]];

(* Make an array of last point of all data sets *)
myEndDat=endDat[[#]]&/@end;

(*Complement output data with the energy *)
ix["energy"] = ix["charge"]+1;

Flatten[{ #, (Sqrt[(mec2N^2)+((#[[ix["cps"]]])^2)+((#[[ix["cpx"]]])^2)+((#[[ix["cpy"]]])^2)])}]&/@myEndDat

]

Basic Procedures

BinParticles[ibin_Integer (*column index for x positions*),iweight_Integer(*column inxed for weight*),nbins_Integer,dat_]:=Module[{max,min,bins,Δ},bins=Table[0,{nbins}];
min=Min[#[[ibin]]&/@dat];
max=Max[#[[ibin]]&/@dat];
Δ=(max-min)/(nbins-1);If [Δ==0,Print["not enough to bin. Aborting..."]; Abort[];];
(*min and max data points are centered in the end bins*)Scan[bins[[Ceiling[(#[[ibin]]-min)/Δ+0.5]]]+=#[[iweight]]&,dat];
Table[{min+(b-1) Δ,1/Δ bins[[b]]},{b,1,nbins}] (*bin center,counts/binwidth*)
];


Clear[BunchHistogram];
BunchHistogram[dat_,dimension_,totalcharge_, nbins_: 250,plotstyle_: Red, plotrange_:{All,All}]:=Module[{xC,weightC,totweight, xlabel,ylabel,xscale,yscale,pbins,filename},If[Length[dat]<5,Abort[]];
weightC= ix["charge"];
totweight = totalcharge;(* Total[dat[[All, ix["charge"]]]];*)
{xC,xlabel,xscale,ylabel,yscale}=Switch[dimension,
"energy",{ix["energy"],"E ( MeV )",(10^(-6)) (*MeV/eV*),"(dI/I)/dE ((MeV^(-1)))",(1/totweight)  (10^6)(*eV/MeV*)},
"x",{ix["x"],"x ( cm )",(10^2) (*cm/m*),"(dI/I)/dL ((cm^(-1)))",(1/totweight)*(10^(-2))(*m/cm*)},
"y",{ix["y"],"y ( cm )",(10^2) (*cm/m*),"(dI/I)/dL ((cm^(-1)))",(1/totweight)*(10^(-2))(*m/cm*)},
"s",{ix["s"],"s ( m )",1 (*m/m*),"(dI/I)/dL ((cm^(-1)))",(1/totweight)*(10^(-2))(*m/cm*)},
"hitθ",{ix["hitθ"],"θ (rad)",1 (* *),"(dI/I)/dθ ((rad^(-1)))",(1/totweight)(* *)}

];
pbins={xscale #1,yscale #2}&@@@BinParticles[xC,weightC,nbins,dat];
(*pbins=Select[pbins,#[[2]]>0&];*)
ListLogPlot[pbins,PlotRange->plotrange,FrameLabel->{xlabel,ylabel},PlotStyle->plotstyle,Axes->False,AspectRatio->1/2, Joined->True]

]

SetupHistograms[pdat_]:=Module[{Lbuffer, currentTotal, currentL0, currentL1, currentC, minL, maxL},
minL = Min[pdat[[All, ix["s"] ]] ];
maxL = Max[pdat[[All, ix["s"] ]] ];
(*Get particles within Lbuffer or the ends*)
Lbuffer = 0.01;
datL0=Select[pdat, #[[ix["s"]]]< minL + Lbuffer&];
datL1=Select[pdat, #[[ix["s"]]]> maxL - Lbuffer&];
(*Particles lost in the center (bulk) of the cryomodule *)
datC=Select[pdat, minL + Lbuffer <  #[[ix["s"]]] <  maxL - Lbuffer&];

{
currentTotal =  Total[#[[ix["charge"]]] &/@pdat],
currentL0 = Total[#[[ix["charge"]]] &/@datL0],
currentL1 = Total[#[[ix["charge"]]] &/@datL1],
currentC = Total[#[[ix["charge"]]] &/@datC]
};
totalcharge=currentTotal;


Grid[{
{"Total current:", (10^9)frequency *currentTotal, "nA"},
{"Current Lost"},
{"% Left", "% center", "% Right"},
(100/currentTotal){currentL0, currentC, currentL1}}, Frame->True, Alignment->"."]


]

BinParticles2D[{ibin_Integer, jbin_Integer, iweight_Integer},{xbins_Integer,ybins_Integer}, dat_]:=Module[{xmax,xmin,ymax, ymin, bins, Δx, Δy},
bins=Table[0,{xbins},{ybins} ];
xmin=Min[#[[ibin]]&/@dat];
xmax=Max[#[[ibin]]&/@dat];
ymin=Min[#[[jbin]]&/@dat];
ymax=Max[#[[jbin]]&/@dat];

Δx = ((xmax-xmin)/(xbins-1));Δy = ((ymax-ymin)/(ybins-1));
(* min and max data points are centered in the end bins *)

Scan[bins[[Ceiling[ ((#[[ibin]] - xmin)/(Δx))+0.5], Ceiling[ ((#[[jbin]] - ymin)/(Δy))+0.5]  ]] += #[[iweight]]& , dat];
Flatten[Table[{xmin + (bx-1)Δx ,ymin + (by-1)Δy  ,   (1/(Δx Δy)) bins[[bx, by]]}, {bx,1,xbins},  {by,1,ybins}] ,1](* bin center x, y, counts/binarea *)
]

Clear[BunchPhaseSpaceDensity];
BunchPhaseSpaceDensity[bunch_, {xaxis_String,yaxis_String}, {xbins_, ybins_},opts:OptionsPattern[]]:=Module[{ixx, jxx, bins2d,xlabel, ylabel,xscale, yscale,  bins2dmax, scaledbins2d, GetAxis},


GetAxis[axis_]:=Switch[axis,
"x", {ix["x"],"x ( μm )",(10^6)},
"cpx", {ix["cpx"],"(Subscript[cp, x]) ( eV )",1},
"y", {ix["y"],"y ( μm )",(10^6)},
"cpy", {ix["cpy"],"(Subscript[cp, y]) ( eV )",1},
"s", {ix["s"],"s ( m )",1},
"cps", {ix["cps"],"(Subscript[cp, s]) ( MeV )",(10^(-6))},
"t", {ix["t"],"t ( ps )",(10^12)}
];

{ixx,xlabel,xscale}=GetAxis[xaxis];

{jxx,ylabel,yscale}=GetAxis[yaxis];

bins2d=BinParticles2D[{ixx,jxx,ix["charge"]}, {xbins,ybins},bunch];
bins2dmax=Max[#[[3]]&/@bins2d];
scaledbins2d={xscale#1, yscale  #2, (1/bins2dmax) #3}&@@@bins2d;
ListDensityPlot[scaledbins2d,Evaluate[FilterRules[{opts}, Options[ListDensityPlot]]], ColorFunction->(ColorData["AvocadoColors"][#1]&),
BaseStyle->{White, FontSize->18},Background->Black,ColorFunctionScaling->False, FrameLabel->{xlabel, ylabel}, ImageSize->500,PlotRange->All ]

]

Tests

solenoid

workdir=root<>"solenoid/"

gWall=Show[WallPlot[workdir<>"wall.out"], ImageSize->800, AspectRatio->Automatic]

p0c=55 (10^6);
s0=0.5;
myt= 0;
myr= 0.005;
freq=1.3 (10^9);
tmax = 1/freq;
t = 0;
dat=
(*x px y py s ps t Q*)
Table[ {r, 0, 0, 0, s0,p0c, t,1.0},{r,0,0.001, 0.0002}];

file="ReferenceParticles.dat";
Export[workdir<>file, {{Length[dat]}}~Join~dat,"Table"];
Length[dat]

myTrackDat=ImportTracks[workdir<>"ReferenceParticles.out"];
Show[EnergyTrackPlot[myTrackDat], AspectRatio->Automatic, ImageSize->1000, BaseStyle->{FontSize->20}, PlotRange->All]

myf={#[[ix["s"]]],#[[ix["x"]]]}&;
ListPlot[Table[myf/@myTrackDat[[i]], {i,Length[myTrackDat]}], Joined->True ]

myf={#[[ix["s"]]],#[[ix["y"]]]}&;
ListPlot[Table[myf/@myTrackDat[[i]], {i,Length[myTrackDat]}], Joined->True ]

Lzf= #[[ix["x"]]] #[[ix["cpy"]]]-  #[[ix["y"]]] #[[ix["cpx"]]]&;
myf={#[[ix["s"]]],Lzf[#]}&;
ListPlot[Table[myf/@myTrackDat[[i]], {i,Length[myTrackDat]}], Joined->True , FrameLabel->{"s", "Lz"}]

rf = Sqrt[ ((#[[ix["x"]]])^2)+ ((#[[ix["y"]]])^2)]&;
Bzf= #[[ix["Bs"]]]&;
myf={#[[ix["s"]]],Bzf[#]((rf[#])^2)}&;
ListPlot[Table[myf/@myTrackDat[[i]], {i,Length[myTrackDat]}], Joined->True , FrameLabel->{"s", "Lz"}]

/

e_gun

workdir=root<>"e_gun/"

gWall=Show[WallPlot[workdir<>"wall.out"], ImageSize->800, AspectRatio->Automatic]

p0c=0;
s0=0.000001;
myt= 0;
myr= 0.005;
freq=1.3 (10^9);
tmax = 1/freq;

dat=
(*x px y py s ps t Q*)
Flatten[Table[ {r, 0, 0, 0, s0,0, t,1.0},{r,0,0.01, 0.001}, {t,-180tmax/360.,180tmax/360, 10tmax/360}], 1];

file="ReferenceParticles.dat";
Export[scratchdir<>file, {{Length[dat]}}~Join~dat,"Table"];
Run["mv "<>scratchdir<>file<>" "<>workdir<>file];
Length[dat]

myTrackDat=ImportTracks[workdir<>"ReferenceParticles.out"];
Show[EnergyTrackPlot[myTrackDat], AspectRatio->Automatic, ImageSize->1000, BaseStyle->{FontSize->20}, PlotRange->All]

cavity2

workdir=root<>"cavity2/"

gWall=WallPlot[workdir<>"wall.out"];
Show[gWall, ImageSize->500]

Show[WallNormalPlot[workdir<>"wall.out", 0.01], gWall, ImageSize->500]

monitordat=Import[workdir<>"l1_n100_p50.monitor"][[3;;]];

mbeg=Select[monitordat, #[[ix["s"]]]<0.2&];Length[mbeg]

mend=Select[monitordat, #[[ix["s"]]]>0.2&];Length[mend]

myTrackDat=ImportTracks[workdir<>"l1_n100_p50.out"];

Show[EnergyTrackPlot[myTrackDat], AspectRatio->Automatic, ImageSize->1000, BaseStyle->{FontSize->20}]

injector_crmod

workdir=root<>"injector_crmod/"

gWall=WallPlot[workdir<>"wall.out"];
Show[gWall, ImageSize->500]

Show[WallNormalPlot[workdir<>"wall.out", 0.0001], gWall, ImageSize->500]

monitordat=Import[workdir<>"l1_n100_p50.monitor"][[3;;]];

mbeg=Select[monitordat, #[[ix["s"]]]<0.2&];Length[mbeg]

5

mend=Select[monitordat, #[[ix["s"]]]>0.2&];Length[mend]

myTrackDat=ImportTracks[workdir<>"l1_n100_p50.out"];

Show[EnergyTrackPlot[myTrackDat], AspectRatio->1/3, ImageSize->1500]

cavity7

workdir=root<>"cavity7/"

gWall=Show[WallPlot[workdir<>"wall.out"], ImageSize->800]

Show[WallNormalPlot[workdir<>"wall.out", 0.01], gWall, ImageSize->800]

monitordat=Import[workdir<>"l1_n100_p50.monitor"][[3;;]];

mbeg=Select[monitordat, #[[ix["s"]]]<1&];Length[mbeg]

mend=Select[monitordat, #[[ix["s"]]]>1&];Length[mend]

myTrackDat=ImportTracks[workdir<>"l1_n100_p50.out"];

Show[EnergyTrackPlot[myTrackDat],ImageSize->1000, AspectRatio->Automatic]

mbeg

cebaf_cavity7

workdir=root<>"cebaf_cavity7/"

gWall=Show[WallPlot[workdir<>"wall.out"], ImageSize->800]

Show[WallNormalPlot[workdir<>"wall.out", 0.01], gWall, ImageSize->800]

monitordat=Import[workdir<>"l1_n100_p50.monitor"][[3;;]];

mbeg=Select[monitordat, #[[ix["s"]]]<1&];Length[mbeg]

mend=Select[monitordat, #[[ix["s"]]]>1&];Length[mend]

myTrackDat=ImportTracks[workdir<>"l1_n100_p50.out"];
Dimensions[myTrackDat]

Show[EnergyTrackPlot[myTrackDat], AspectRatio->Automatic, ImageSize->1200]

cavity9

workdir=root<>"cavity9/"

gWall=Show[WallPlot[workdir<>"wall.out"], ImageSize->800]

Show[WallNormalPlot[workdir<>"wall.out", 0.01], gWall, ImageSize->800]

monitordat=Import[workdir<>"l1_n100_p50.monitor"][[3;;]];

mbeg=Select[monitordat, #[[ix["s"]]]<1&];Length[mbeg]

mend=Select[monitordat, #[[ix["s"]]]>1&];Length[mend]

myTrackDat=ImportTracks[workdir<>"l1_n200_p50.out"];
Dimensions[myTrackDat]

Show[EnergyTrackPlot[myTrackDat], AspectRatio->Automatic, ImageSize->1200]

htc

workdir=root<>"htc/"

gWall=WallPlot[workdir<>"wall.out"];Show[gWall, ImageSize->800]

Show[WallNormalPlot[workdir<>"wall.out", 0.01], gWall, ImageSize->800]

endDat=LoadEndParticles[workdir<>"l1_n100_p50.out"];Length[endDat]

SetupHistograms[endDat]

prange={{0, All},All};
Show[BunchHistogram[datC,"energy",totalcharge,100,Black,prange],BunchHistogram[Join[datL0,datL1],"energy",totalcharge,20,Red,prange],ImageSize->500]

prange={All,All};
Show[BunchHistogram[datC,"s",totalcharge,200,Black,prange],ImageSize->500]

ListPlot[{ #[[ix["s"]]],  #[[ix["x"]]]}&/@datC, PlotStyle->{PointSize[Tiny], Black}, AspectRatio->Automatic]

myTrackDat=ImportTracks[workdir<>"l1_n100_p50.out"];

Show[EnergyTrackPlot[myTrackDat], AspectRatio->Automatic, ImageSize->1000]

.

crmod

workdir=root<>"crmod/"

gWall=WallPlot[workdir<>"wall.out"];Show[gWall, ImageSize->1000]

Show[WallNormalPlot[workdir<>"wall.out", 0.001], gWall, ImageSize->500, PlotRange->{{2,2.1},{0,0.1}}]

monitordat=Import[workdir<>"l1_n1000_p500.monitor", "Table"][[3;;]];

mbeg=Select[monitordat, #[[ix["s"]]]<1&];Length[mbeg]

mend=Select[monitordat, #[[ix["s"]]]>1&];Length[mend]

ListLogPlot[ {
{#[[ix["s_origin"]]],#[[ix["charge"]]]}&/@mbeg,
{#[[ix["s_origin"]]],#[[ix["charge"]]]}&/@mend}, PlotStyle->{Red, Blue}]

ListLogPlot[ {
{#[[ix["cps"]]],#[[ix["charge"]]]}&/@mbeg,
{#[[ix["cps"]]],#[[ix["charge"]]]}&/@mend
}, PlotStyle->{Red, Blue}, PlotRange->{All,All}]

myTrackDat=ImportTracks[workdir<>"l1_n100_p50.out"];
myTrackDat=Sort[myTrackDat,energyf[Last[#1]]<energyf[Last[#2]]&];
Dimensions[myTrackDat]

Show[EnergyTrackPlot[myTrackDat], AspectRatio->Automatic,  ImageSize->1000]

lefttracks=Select[myTrackDat, #[[-1, ix["s"]]]<0.1&];
righttracks=Select[myTrackDat, #[[-1, ix["s"]]]>8.8&];

Show[EnergyTrackPlot[lefttracks], AspectRatio->1/4, ImageSize->1400, BaseStyle->{FontSize->20}]

Show[EnergyTrackPlot[righttracks], AspectRatio->1/4, ImageSize->1400, BaseStyle->{FontSize->20}]

drift

workdir=root<>"drift/"

gWall=Show[WallPlot[workdir<>"wall.out"], ImageSize->800]

small=.000001;
dat={
{2},
{0, 0, 0, 0, 0, 1+small, 0, 1},
{0, 0, 0, 0, 1.0-small, -1, 0, 1}};
Export[workdir<>"ReferenceParticles.dat", dat,"Table"];

p0c=15 (10^6);
s0=0.5;
dat=
Table[ {0, p0c Sin[θ], 0, 0, s0, p0c  Cos[θ], 0,1},{θ,-π, π,10 π/180.}];

Export[workdir<>"ReferenceParticles.dat", {{Length[dat]}}~Join~dat,"Table"];
Length[dat]

myTrackDat=ImportTracks[workdir<>"ReferenceParticles.out"];
Show[EnergyTrackPlot[myTrackDat],AspectRatio->1  ,ImageSize->500]

tableTrack=Table[{#[[ix["t"]]],#[[ix["s"]]]}&/@myTrackDat[[i]],{i,1,Length[myTrackDat]}];

ListPlot[tableTrack, Joined->True,PlotRange->All]

mask

workdir=root<>"mask/"

gWall=Show[WallPlot[workdir<>"wall.out"], ImageSize->800]

small=.000001;
dat={
{2},
{0, 0, 0, 0, 0, 1+small, 0, 1},
{0, 0, 0, 0, 1.0-small, -1, 0, 1}};
Export[workdir<>"ReferenceParticles.dat", dat,"Table"];

p0c=-15 (10^6);
s0=1.99;
mm = .001;
dat=
Flatten[Table[ {x ,0, y, 0, s0, p0c , 0,1},{x,-40 mm,100mm,2mm}, {y, -30mm,30mm, 2mm}],1];

Export[workdir<>"ReferenceParticles.dat", {{Length[dat]}}~Join~dat,"Table"];
Length[dat]

myTrackDat=ImportTracks[workdir<>"ReferenceParticles.out"];
Show[EnergyTrackPlot[myTrackDat],AspectRatio->1  ,ImageSize->500]

xyz=Table[Line[{#[[ix["s"]]],#[[ix["x"]]], #[[ix["y"]]]}&/@t],{t, myTrackDat}];
Graphics3D[xyz, Axes->True]

multipass

workdir=root<>"multipass/"

gWall=Show[WallPlot[workdir<>"wall.out"], ImageSize->800]

small=.000001;
dat={
{2},
{0, 0, 0, 0, 0, 1+small, 0, 1},
{0, 0, 0, 0, 1.0-small, -1, 0, 1}};
Export[workdir<>"ReferenceParticles.dat", dat,"Table"];

p0c=15 (10^6);
s0=2.000000001;
myt= ((Random[])/(1.3 (10^9)));
θdat = Table[θ, {θ, -π,π,π/1800.}];
dat=
Table[ {0, p0c Sin[θ], 0, 0, s0, p0c  Cos[θ], myt,1},{θ,θdat }];

file="ReferenceParticles.dat";
Export[scratchdir<>file, {{Length[dat]}}~Join~dat,"Table"];
Run["cp "<>scratchdir<>file<>" "<>workdir<>file];
Length[dat]

myTrackDat=ImportTracks[workdir<>"ReferenceParticles.out"];
Show[EnergyTrackPlot[myTrackDat],AspectRatio->Automatic  ,ImageSize->500, PlotRange->All]

Field query

myt =0/(1.3 (10^9));
myx = 0.05;
mys=0.99999;
query=Table[ {myt,x,0,mys}, {x,-(10^(-7)), (10^(-7)),(10^(-7))/2.}];
file="querypoints.dat";
Export[scratchdir<>file, {{Length[query]}}~Join~query,"Table"];
Run["cp "<>scratchdir<>file<>" "<>workdir<>file];

dat = Import[workdir <> "querypoints.out", "Table"][[2 ;;]];
ListPlot[{#[[2]], #[[5]]} & /@ dat, PlotStyle -> {Red, Black}]

dat = Import[workdir <> "querypoints.out", "Table"][[2 ;;]]; 1
ListPointPlot3D[{#[[2]], #[[4]], #[[5]]} & /@ dat, PlotStyle -> {Red, Black}]

ion

workdir=root<>"ion/"

gWall=Show[WallPlot[workdir<>"wall.out"], ImageSize->800]

small=.000001;
dat={
{2},
{0, 0, 0, 0, 0, 1+small, 0, 1},
{0, 0, 0, 0, 1.0-small, -1, 0, 1}};
Export[workdir<>"ReferenceParticles.dat", dat,"Table"];

p0c= (10^5);
x0 = 0.05;
s0=0.5;
dat=
Table[ {x0, p0c Sin[θ], 0, 0, s0, p0c  Cos[θ], 0,1},{θ,{π/4}}];

Export[workdir<>"ReferenceParticles.dat", {{Length[dat]}}~Join~dat,"Table"];
Length[dat]

myTrackDat=ImportTracks[workdir<>"ReferenceParticles.out"];
Show[EnergyTrackPlot[myTrackDat],AspectRatio->1  ,ImageSize->500]

tableTrack=Table[{#[[ix["t"]]],#[[ix["s"]]]}&/@myTrackDat[[i]],{i,1,Length[myTrackDat]}];

ListPlot[tableTrack, Joined->True,PlotRange->All]

sbend

workdir=root<>"sbend/"

gWall=Show[WallPlot[workdir<>"wall.out"], ImageSize->800, AspectRatio->Automatic]

p0c=90 (10^6);
s0=0.01;
myt= 0;
θdat = Table[θ, {θ, {0.0}}];
dat=
Table[ {0, p0c Sin[θ], 0, 0, s0, p0c  Cos[θ], myt,1},{θ,θdat }];

file="ReferenceParticles.dat";
Export[scratchdir<>file, {{Length[dat]}}~Join~dat,"Table"];
Run["cp "<>scratchdir<>file<>" "<>workdir<>file];
Length[dat]

myTrackDat=ImportTracks[workdir<>"ReferenceParticles.out"];
Show[EnergyTrackPlot[myTrackDat], AspectRatio->Automatic, ImageSize->1000, BaseStyle->{FontSize->20}, PlotRange->All]

myTrackDat

energyf[#]&/@myTrackDat

lcavity

A bmad Lcavity with no field grid.

workdir=root<>"lcavity/"

gWall=Show[WallPlot[workdir<>"wall.out"], ImageSize->800, AspectRatio->Automatic]

p0c=1000 (10^6);
s0=0;
myt= 0;
θdat = Table[θ, {θ, {0}}];
dat=
Table[ {0, p0c Sin[θ], 0, 0, s0, p0c  Cos[θ], myt,1},{θ,θdat }];

file="ReferenceParticles.dat";
Export[scratchdir<>file, {{Length[dat]}}~Join~dat,"Table"];
Run["cp "<>scratchdir<>file<>" "<>workdir<>file];
Length[dat]

myTrackDat=ImportTracks[workdir<>"ReferenceParticles.out"];
Show[EnergyTrackPlot[myTrackDat], AspectRatio->Automatic, ImageSize->1000, BaseStyle->{FontSize->20}, PlotRange->All]

myTrackDat

energyf[#]&/@myTrackDat

linac_short

workdir=root<>"linac_short/"

sEnd=Import[workdir<>"wall.out", "Table"][[-1, ix["s"]]]

gWall=WallPlot[workdir<>"wall.out"]; Show[gWall, ImageSize->1500]

myTrackDat=ImportTracks[workdir<>"l1_n100_p50.out"];
myTrackDat=Sort[myTrackDat,energyf[Last[#1]]<energyf[Last[#2]]&];
Dimensions[myTrackDat]

Show[EnergyTrackPlot[myTrackDat], AspectRatio->1/4, ImageSize->1500, BaseStyle->{FontSize->20}, PlotRange->{{-1,21},{-.11,.11}} ]

lefttracks=Select[myTrackDat, #[[-1, ix["s"]]]<0.01&];
righttracks=Select[myTrackDat, #[[-1, ix["s"]]]>sEnd-0.01&];

{Length[lefttracks], Length[righttracks]}

Total[#[[-1, ix["cps"]]]&/@lefttracks]

ListPlot[{
{#[[1, ix["s"]]],(10^(-6))#[[-1, ix["cps"]]]}&/@lefttracks,
{#[[1, ix["s"]]],(10^(-6))#[[-1, ix["cps"]]]}&/@righttracks}, PlotStyle->{Red, Blue}, PlotRange->150{-1,1}, FrameLabel->{"Original s (m)", "Exiting c (Subscript[p, s]) (MeV)"}, AspectRatio->1/3, BaseStyle->{FontSize->20}, ImageSize->800]

Total[#[[-1, ix["cps"]]]&/@righttracks]

gLeft=Show[EnergyTrackPlot[lefttracks], AspectRatio->1/4, ImageSize->1500, BaseStyle->{FontSize->20}, PlotRange->{{-1,21},{-.11,.11}}]

gRight=Show[EnergyTrackPlot[righttracks], AspectRatio->1/4, ImageSize->1500, BaseStyle->{FontSize->20}, PlotRange->{{-1,21},{-.11,.11}}]

Export[groot<>"crmod2_left_bad.pdf", gLeft];
Export[groot<>"crmod2_right_bad.pdf", gRight];

lefttracksfromright=Select[myTrackDat, (#[[-1, ix["s"]]]<10) &&  #[[1, ix["s"]]]>10&];Length[lefttracksfromright]

righttracksfromleft=Select[myTrackDat, (#[[-1, ix["s"]]]>10) &&  #[[1, ix["s"]]]<10&];Length[righttracksfromleft]

testtracks=Select[myTrackDat, (#[[-1, ix["s"]]]<10) &&  #[[1, ix["s"]]]>10&];Length[testtracks]

Show[EnergyTrackPlot[lefttracksfromright], AspectRatio->1/4, ImageSize->1500, BaseStyle->{FontSize->21}]

Show[EnergyTrackPlot[righttracksfromleft], AspectRatio->1/4, ImageSize->1500, BaseStyle->{FontSize->20}]

MomentumEvolution[track_]:=Module[{},
{#[[ix["s"]]],(10^(-6)) #[[ix["cps"]]]}&/@track
]

g1=ListPlot[MomentumEvolution/@righttracksfromleft[[1;;]], Joined->True, FrameLabel->{"s", "c (Subscript[p, s]) (MeV)"}, PlotStyle->Blue];
g2=ListPlot[MomentumEvolution/@lefttracksfromright[[1;;]], Joined->True, FrameLabel->{"s", "c (Subscript[p, s]) (MeV)"}, PlotStyle->Red];
Show[g1,g2, ImageSize->800 , BaseStyle->{FontSize->20}]

g1=ListPlot[MomentumEvolution/@lefttracks, Joined->True, FrameLabel->{"s (m)", "c (Subscript[p, s]) (MeV)"}, PlotStyle->Blue];
g2=ListPlot[MomentumEvolution/@righttracks, Joined->True, FrameLabel->{"s (m)", "c (Subscript[p, s]) (MeV)"}, PlotStyle->Red];
Show[g1,g2, ImageSize->800 , BaseStyle->{FontSize->20}, PlotRange->{-160,160}]

linac_short_spacer

workdir=root<>"linac_short_spacer/"

sEnd=Import[workdir<>"wall.out", "Table"][[-1, ix["s"]]]

gWall=WallPlot[workdir<>"wall.out"]; Show[gWall, ImageSize->1500]

myTrackDat=ImportTracks[workdir<>"l1_n100_p50.out"];
myTrackDat=Sort[myTrackDat,energyf[Last[#1]]<energyf[Last[#2]]&];
Dimensions[myTrackDat]

Show[EnergyTrackPlot[myTrackDat], AspectRatio->1/8, ImageSize->1500, BaseStyle->{FontSize->20}, PlotRange->{{-1,21},{-.11,.11}} ]

lefttracks=Select[myTrackDat, #[[-1, ix["s"]]]<0.01&];
righttracks=Select[myTrackDat, #[[-1, ix["s"]]]>sEnd-0.01&];

{Length[lefttracks], Length[righttracks]}

Total[#[[-1, ix["cps"]]]&/@lefttracks]

ListPlot[{
{#[[1, ix["s"]]],(10^(-6))#[[-1, ix["cps"]]]}&/@lefttracks,
{#[[1, ix["s"]]],(10^(-6))#[[-1, ix["cps"]]]}&/@righttracks}, PlotStyle->{Red, Blue}, PlotRange->150{-1,1}, FrameLabel->{"Original s (m)", "Exiting c (Subscript[p, s]) (MeV)"}, AspectRatio->1/3, BaseStyle->{FontSize->20}, ImageSize->800]

Total[#[[-1, ix["cps"]]]&/@righttracks]

gLeft=Show[EnergyTrackPlot[lefttracks], AspectRatio->1/4, ImageSize->1500, BaseStyle->{FontSize->20}, PlotRange->{{-1,21},{-.11,.11}}]

gRight=Show[EnergyTrackPlot[righttracks], AspectRatio->1/4, ImageSize->1500, BaseStyle->{FontSize->20}, PlotRange->{{-1,21},{-0.11,.11}}]

Export[groot<>"crmod2_left_good.pdf", gLeft];
Export[groot<>"crmod2_right_good.pdf", gRight];

lefttracksfromright=Select[myTrackDat, (#[[-1, ix["s"]]]<10) &&  #[[1, ix["s"]]]>10&];Length[lefttracksfromright]

righttracksfromleft=Select[myTrackDat, (#[[-1, ix["s"]]]>10) &&  #[[1, ix["s"]]]<10&];Length[righttracksfromleft]

testtracks=Select[myTrackDat, (#[[-1, ix["s"]]]<10) &&  #[[1, ix["s"]]]>10&];Length[testtracks]

Show[EnergyTrackPlot[lefttracksfromright], AspectRatio->1/4, ImageSize->1500, BaseStyle->{FontSize->20}]

Show[EnergyTrackPlot[righttracksfromleft], AspectRatio->1/4, ImageSize->1500, BaseStyle->{FontSize->20}]

MomentumEvolution[track_]:=Module[{},
{#[[ix["s"]]],(10^(-6)) #[[ix["cps"]]]}&/@track
]

g1=ListPlot[MomentumEvolution/@righttracksfromleft[[1;;]], Joined->True, FrameLabel->{"s", "c (Subscript[p, s]) (MeV)"}];
g2=ListPlot[MomentumEvolution/@lefttracksfromright[[1;;]], Joined->True, FrameLabel->{"s", "c (Subscript[p, s]) (MeV)"}];
Show[g1,g2, ImageSize->800 , BaseStyle->{FontSize->20}, PlotRange->{-150,150}]

g1=ListPlot[MomentumEvolution/@lefttracks, Joined->True, FrameLabel->{"s (m)", "c (Subscript[p, s]) (MeV)"}, PlotStyle->Blue];
g2=ListPlot[MomentumEvolution/@righttracks, Joined->True, FrameLabel->{"s (m)", "c (Subscript[p, s]) (MeV)"}, PlotStyle->Red];
Show[g1,g2, ImageSize->800 , BaseStyle->{FontSize->20}, PlotRange->{All, {-160,160}}]

sbend_3D

workdir=root<>"sbend_3D/"

gWall=Show[WallPlot[workdir<>"wall.out"], ImageSize->800]

small=.000001;
dat={
{2},
{0, 0, 0, 0, 0, 1+small, 0, 1},
{0, 0, 0, 0, 1.0-small, -1, 0, 1}};
Export[workdir<>"ReferenceParticles.dat", dat,"Table"];

p0c=10 (10^6);
s0=0.1000000001;
dat=
Table[ {0, p0c Sin[θ], 0, 0, s0, p0c  Cos[θ], 0,1},{θ,-π/1000,π/1000,π/2000.}];

Export[workdir<>"ReferenceParticles.dat", {{Length[dat]}}~Join~dat,"Table"];
Length[dat]

myTrackDat=ImportTracks[workdir<>"ReferenceParticles.out"];
Show[EnergyTrackPlot[myTrackDat],AspectRatio->1  ,ImageSize->500]

tableTrack=Table[{#[[ix["s"]]],#[[ix["By"]]]}&/@myTrackDat[[i]],{i,1,Length[myTrackDat]}];
ListPlot[tableTrack, Joined->True,PlotRange->All]

tableTrack=Table[{#[[ix["s"]]],#[[ix["x"]]]}&/@myTrackDat[[i]],{i,1,Length[myTrackDat]}];
ListPlot[tableTrack, Joined->True,PlotRange->All]

Linac Monitors

Load Data: linac_bad_quadsoff

workdir=root<>"../../../../erl/DarkCurrent/real_linac/linac_bad_quadsoff/";

Nmonitors=10;

header=Import[workdir<>"l1_n500_p250.monitor","Table"][[1;;2]]//Grid

monitordat=Flatten[Table[Import[workdir<>"l"<>ToString[i]<>"_n500_p250.monitor","Table","HeaderLines"->2],{i,1,Nmonitors}],1];Dimensions[monitordat]

Load Data: linac_bad

workdir=root<>"../../../../erl/DarkCurrent/real_linac/linac_bad/"

Nmonitors=10;

header=Import[workdir<>"l1_n2000_p500.monitor","Table"][[1;;2]]//Grid

monitordat=Flatten[Table[Import[workdir<>"l"<>ToString[i]<>"_n2000_p500.monitor","Table","HeaderLines"->2],{i,1,Nmonitors}],1];Dimensions[monitordat]

Load Data: linac_good

workdir=root<>"../../../../erl/DarkCurrent/real_linac/linac_good/";

Nmonitors=10;

header=Import[workdir<>"l1_n4000_p1000.monitor","Table"][[1;;2]]//Grid

monitordat=Flatten[Table[Import[workdir<>"l"<>ToString[i]<>"_n4000_p1000.monitor","Table","HeaderLines"->2],{i,1,Nmonitors}],1];Dimensions[monitordat]

Plots

(*Find monitors*)
slist=monitordat[[All,ix["s"]]];
slistno0=Select[slist,#>1&];
slistAll={0}~Join~Intersection[slistno0,slistno0][[1;;]];ListPlot[slistAll]

crmoddat=Table[Select[monitordat,(Abs[#[[ix["s"]]]-slistAll[[i]]]<1)&],{i,1,Length[slistAll]}];
ListPlot[Length[#]&/@crmoddat, PlotRange->All]

Min[#[[ix["cps"]]]&/@monitordat]

Max[#[[ix["cps"]]]&/@monitordat]

g1=ListPlot[
Table[{#[[ix["s_origin"]]],10^-6 #[[ix["cps"]]]}&/@crmoddat[[i]],
{i,1,Length[crmoddat]}],AspectRatio->1/2,ImageSize->500,PlotRange->All,FrameLabel->{"(Subscript[s, origin]) (m)","c \!\(\*SubscriptBox[\(p\), \(s\)]\) exit (MeV)"}, PlotStyle->{Red,Blue}];Show[g1,ImageSize->1000, BaseStyle->{FontSize->25}, PlotRange->{{0,380},{-3000,3000}}]

(*mod with the RF period*)
mbeg=Flatten[{Mod[(#[[ix["t"]]]), 1/frequency], #[[2;;]]}]&/@First[crmoddat];
mend=Flatten[{Mod[(#[[ix["t"]]]), 1/frequency], #[[2;;]]}]&/@Last[crmoddat];

(*Complement output data with the energy *)
ix["energy"] = ix["s_origin"]+1;
mbeg=Flatten[{ #, (Sqrt[(mec2N^2)+((#[[ix["cps"]]])^2)+((#[[ix["cpx"]]])^2)+((#[[ix["cpy"]]])^2)])}]&/@mbeg;
mend=Flatten[{ #, (Sqrt[(mec2N^2)+((#[[ix["cps"]]])^2)+((#[[ix["cpx"]]])^2)+((#[[ix["cpy"]]])^2)])}]&/@mend;

mendhighenergy=Select[mend, #[[ix["energy"]]] > 2 (10^9)&];
mendlowenergy=Select[mend, #[[ix["energy"]]] < 2 (10^9)&];



cp0Beg = 15. (10^6);
cp0End = 2739. (10^6);
trefend =Mod[1.1835 (10^(-6)), 1/frequency];

Grid[{
{"Beg current (nA): ", (10^12)frequency*Total[mbeg[[All, ix["charge"]]]]},
{"End current (nA): ", (10^12)frequency*Total[mend[[All, ix["charge"]]]]},
{"End high en. (nA): ", (10^12)frequency*Total[mendhighenergy[[All, ix["charge"]]]]}}, Frame->True]

Show[
BunchHistogram[mbeg, "energy", (*tot charge*)frequency*Total[mbeg[[All, ix["charge"]]]],50,Red],
BunchHistogram[mend, "energy", (*tot charge*)frequency*Total[mend[[All, ix["charge"]]]],50,Blue],
 ImageSize->700, BaseStyle->{FontSize->20}, PlotRange->{{0,3000}, All}]

Show[BunchHistogram[mbeg, "energy", (*tot charge*)frequency*Total[mbeg[[All, ix["charge"]]]],100,Red], ImageSize->700, BaseStyle->{FontSize->20}, PlotRange->{{0,400},Automatic}]

Show[BunchHistogram[mendlowenergy, "energy", (*tot charge*)frequency*Total[mbeg[[All, ix["charge"]]]],50,Blue], ImageSize->700, BaseStyle->{FontSize->20}, PlotRange->{{0,160},Automatic}]

Show[BunchHistogram[mendhighenergy, "energy", (*tot charge*)frequency*Total[mbeg[[All, ix["charge"]]]],20,Blue], ImageSize->700, BaseStyle->{FontSize->20}, PlotRange->{All,Automatic}]

ListPlot[{ #[[ix["x"]]],  #[[ix["cpx"]]]}&/@mend, AspectRatio->1, ImageSize->500, Frame->True, PlotRange->All, PlotStyle->{PointSize[Tiny],Black}, PlotStyle->Black, FrameLabel->{"x (m)", "c (Subscript[p, x]) (eV)"}]

mPlotXY[mdat_, color_]:=Module[{},
ListPlot[{ (10^3)#[[ix["x"]]], (10^3) #[[ix["y"]]]}&/@mdat, AspectRatio->1, ImageSize->500, Frame->True, PlotRange->All, PlotStyle->{PointSize[Tiny],color}, PlotStyle->Black, FrameLabel->{"x (mm)", "y (mm)"}, BaseStyle->{FontSize->20}, Axes->False]];
Grid[{{mPlotXY[mbeg, Red],mPlotXY[mend, Blue]}}]

Show[mPlotXY[mendlowenergy,Orange],mPlotXY[mendhighenergy, Purple]]

ListPlot[ {(10^12)#[[ix["t"]]], (10^(-6))#[[ix["energy"]]]}& /@mend, FrameLabel->{"t (ps)", "Energy (MeV)"}, PlotStyle->{PointSize[Tiny],Blue}, ImageSize->600 ,BaseStyle->{FontSize->20}]

BunchPhaseSpaceDensity[ mend, {"t", "cps"}, {200,200}]

BunchPhaseSpaceDensity[ mendhighenergy, {"t", "cps"}, {200,200}]

ListPlot[ {(10^12)#[[ix["t"]]], (10^(-6))#[[ix["energy"]]]}& /@mendhighenergy, FrameLabel->{"t (ps)", "Energy (MeV)"}, PlotStyle->{PointSize[Tiny], Black}]

ListPlot[ {(10^6)#[[ix["x"]]], (10^6)((#[[ix["cpx"]]])/cp0End)}& /@mendlowenergy, FrameLabel->{"x (μm)", "(Subscript[p, x])/(Subscript[p, 0]) (μrad)"}, PlotStyle->{PointSize[Tiny], Black}]

BunchPhaseSpaceDensity[ mendhighenergy, {"x", "y"}, {50,50}]

BunchPhaseSpaceDensity[ mendlowenergy, {"x", "cpx"}, {50,50}]