Initialization

(* Define some constants *)
wall="wall.out";
fileName="l1_n500_p200";
fileIn=fileName<>".dat";
file=fileName<>".out";
numHeaderLines=3 ;  (* 3 lines of header data *)
mc2=510998.928; (* Rest mass of an electron *)

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

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

ixIn["x"]=1;
ixIn["βxγx"]=2;
ixIn["y"]=3;
ixIn["βyγy"]=4;
ixIn["s"]=5;
ixIn["βsγs"]=6;
ixIn["phase"]=7;
ixIn["Q"]=8;

Make Lattice Wall

(* Import Data *)
wallDat=Import[wall,"Table"];

(* Make an array of last point of all data sets *)
myWallDat1=Transpose[{wallDat[[numHeaderLines+1;;All,ix["s"]]],wallDat[[numHeaderLines+1;;All,ix["x"]]]}];
myWallDat2=Transpose[{wallDat[[numHeaderLines+1;;All,ix["s"]]],-wallDat[[numHeaderLines+1;;All,ix["x"]]]}];

(* Plot wall *)
gWall1=ListPlot[myWallDat1,FrameLabel->{"s (m)","x (m)"},Joined->True,PlotStyle->Black];
gWall2=ListPlot[myWallDat2,FrameLabel->{"s (m)","x (m)"},Joined->True,PlotStyle->Black];
gWall=Show[gWall1,gWall2,PlotRange->{{0.3,1.4},All}]

(* Look at dark current origins *)

inDat=Import[fileIn,"Table"];

tablex={#[[ixIn["s"]]],(Sqrt[((#[[ixIn["x"]]])^2)+((#[[ixIn["y"]]])^2)])}&/@inDat[[2;;]];
origins=ListPlot[tablex];

Show[origins,gWall,FrameLabel->{"s (m)","x (m)"},PlotRange->{{0.4,1.3},All}]

Plot Dark Current Tracks

(* Import Data *)
trackDat=Import[file, "Table"];

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

(* Make an array of all data sets and one of phases *)
myTrackDat=trackDat[[begin[[#]];;end[[#]]]]&/@Range[Length[begin]];
phases=#[[1,ix["time"]]]&/@myTrackDat;
phases=(phases-Min[phases])/(Max[phases]-Min[phases]);

(* Make table of energies and scale to [0,1] *)
energies=Table[(Sqrt[(mc2^2)+(((#[[ix["cpx"]]])^2)+((#[[ix["cpy"]]])^2)+((#[[ix["cps"]]])^2))])&/@myTrackDat[[i]],{i,1,Length[myTrackDat]}];
energies =1-(energies-Min[energies])/(Max[energies]-Min[energies]);

(* Plot tracks, colored by energy *)
tableTrack=Table[{#[[ix["s"]]],#[[ix["x"]]]}&/@myTrackDat[[i]],{i,1,Length[myTrackDat]}];
s=Interpolation[tableTrack[[#,All,1]],InterpolationOrder->1]&;
x=Interpolation[tableTrack[[#,All,2]],InterpolationOrder->1]&;
e=Interpolation[energies[[#,All]],InterpolationOrder->1]&;

g=ParametricPlot[{s[#][t],x[#][t]},{t,1,Length[myTrackDat[[#]]]},ColorFunction->Function[{x,y,t},ColorData["Rainbow"][e[#][t]]],ColorFunctionScaling->False]&;
gTrack=Show[g/@Range[Length[myTrackDat]]]

(* Plot tracks, colored by phase *)
tableTrack=Table[{#[[ix["s"]]],#[[ix["x"]]]}&/@myTrackDat[[i]],{i,1,Length[myTrackDat]}];
gTrack=ListPlot[tableTrack,FrameLabel->{"s (m)","x (m)"},Joined->True,PlotStyle->Array[ColorData["RedBlueTones",phases[[#]]]&,Length[phases]],PlotRange->{{7,8.5},All}]

(* Analyze one particle track *)
gTrack1=ListPlot[{#[[ix["time"]]]*1.3*(10^9),#[[ix["x"]]]}&/@myTrackDat[[290]],FrameLabel->{"time (s)","x (m)"},Joined->True,PlotStyle->Array[ColorData["RedBlueTones",phases[[#]]]&,Length[phases]],AspectRatio->1/GoldenRatio]

Plot Dark Current Final Coordinates

(* Import Data *)
endDat=Import[file,"Table"];

(* Use this instead if we've imported trackDat *)
myEndDat=trackDat[[#]]&/@end;

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

(* Plot final coordinates *)
energies=Log[(Sqrt[(mc2^2)+(((#[[ix["cpx"]]])^2)+((#[[ix["cpy"]]])^2)+((#[[ix["cps"]]])^2))])]&/@myEndDat;
minE=Min[energies];
maxE=Max[energies];
tableEnd={#[[ix["s"]]],#[[ix["x"]]],1-(Log[Sqrt[mc2^2+((#[[3]])^2)+((#[[5]])^2)+((#[[7]])^2)]]-minE)/(maxE-minE)}&/@myEndDat;

gEnd=Graphics[{ColorData["Rainbow"][#[[3]]],Point[{#[[1]],#[[2]]}]}&/@tableEnd,Axes->True, PlotRange->All,FrameLabel->{"s (m)","x (m)"}]

Overlay Data with Wall

gTrackWall=Show[{gTrack,gWall},AspectRatio->Automatic,PlotRange->{All,All}]

gEndWall=Show[{gEnd,gWall},PlotRange->{All,All}]

Export[fileName<>"_Track.png",gTrackWall];

Export[fileName<>"_End.png",gEndWall];

Histograms and Density Plots

(* Run Plot Dark Current Final Coordinates first *)
energies=(Sqrt[(mc2^2)+((#[[ix["cps"]]])^2)+((#[[ix["cpx"]]])^2)+((#[[ix["cpy"]]])^2)])&;
angles=#[[ix["hitθ"]]]&/@myEndDat;
(* Create list of indices of particles that make it through lattice *)
passIndex=Intersection[Flatten[Position[myEndDat[[All,ix["hitθ"]]],0.0]],Complement[Range[Length[myEndDat]],Flatten[Position[myEndDat[[All,ix["s"]]],0.0]]]];
(* Create list of indices of particles that hit lattice walls *)
hitIndex=Complement[Range[Length[myEndDat]],Flatten[Position[myEndDat[[All,ix["hitθ"]]],0.0]]];
(* Reduce myEndDat to particles that hit lattice walls *)
hitDat=myEndDat[[#]]&/@hitIndex;

(* Histogram of charges to check fields *)
hQ=Histogram[Transpose[myEndDat][[ix["Q"]]]*(10^34),FrameLabel->{"Charge (C)"},AspectRatio->1/GoldenRatio]

(* Final x coordinates vs. energy, for Val *)
tablexE={#[[ix["x"]]],energies[#]}&/@myEndDat[[hitIndex]];
gxE=ListPlot[tablexE,FrameLabel->{"x (m)","Energy (eV)"},AspectRatio->1/GoldenRatio]

(* Final x coordinate vs. charge, for Val *)
tablexC={#[[ix["x"]]],#[[ix["Q"]]]}&/@myEndDat[[hitIndex]];
gxC=ListPlot[tablexC,FrameLabel->{"x (m)","Charge (C)"},AspectRatio->1/GoldenRatio]

(* Histogram of final x coordinate, for Val *)
hx=Histogram[Transpose[hitDat][[ix["x"]]],FrameLabel->{"x (m)"},AspectRatio->1/GoldenRatio]

(* Histogram of final energy *)
hE = Histogram[{energies[#]&/@myEndDat[[hitIndex]],energies[#]&/@myEndDat[[passIndex]]},FrameLabel->{"energy (eV)"},AspectRatio->1/5,ScalingFunctions->{"Log","Log"}]

(* Histogram of final s coordinate *)
hS = Histogram[Transpose[myEndDat][[ix["s"]]],FrameLabel->{"s (m)"},AspectRatio->1/GoldenRatio,ScalingFunctions->{"Log","Log"}]

(* Histogram of final angles *)
hA=Histogram[Transpose[myEndDat][[ix["hitθ"]]],FrameLabel->{"Hit angle (radians)"},AspectRatio->1/GoldenRatio]

(* Final s coordinate vs. energy *)
tableEs = {#[[ix["s"]]],Log[energies[#]]}&/@myEndDat;
gEs=ListPlot[tableEs,FrameLabel->{"s (m)","Log[energy] (Log[eV])"},Joined->False,AspectRatio->1/GoldenRatio]

(* Final position (x, y) colored by energy, for particles that exited *)
tableEnergies=energies[#]&/@myEndDat[[passIndex]];
minE=Min[tableEnergies];
maxE=Max[tableEnergies];
tableEnd = {#[[ix["y"]]],#[[ix["x"]]],(energies[#]-minE)/(maxE-minE)}&/@myEndDat[[passIndex]];
gExy=Graphics[{ColorData["SolarColors"][#[[3]]],Point[{#[[1]],#[[2]]}]}&/@tableEnd,Axes->True,AspectRatio->1/GoldenRatio]

(* Energy as a function of angle *)
tableAngle={#[[ix["hitθ"]]],energies[#]}&/@myEndDat;
gAngle=ListPlot[tableAngle,FrameLabel->{"Hit Angle (radians)", "Energy (eV)"},Joined->False,AspectRatio->1/GoldenRatio]

Binning by Hand for Weighted Histograms and Density Plots

(* Setup *)
energies=(Sqrt[(mc2^2)+((#[[ix["cps"]]])^2)+((#[[ix["cpx"]]])^2)+((#[[ix["cpy"]]])^2)])&;
(* Create list of indices of particles that make it through lattice *)
passIndex=Intersection[Flatten[Position[myEndDat[[All,ix["hitθ"]]],0.0]],Complement[Range[Length[myEndDat]],Flatten[Position[myEndDat[[All,ix["s"]]],0.0]]]];
(* Create list of indices of particles that hit lattice walls *)
hitIndex=Complement[Range[Length[myEndDat]],Flatten[Position[myEndDat[[All,ix["hitθ"]]],0.0]]];

testDat={#[[6]],energies[#],#[[9]]}&/@myEndDat;

(* Weighted 1D Histogram of energy *)

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);
(*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_,nbins_: 250,color_: Red]:=Module[{xC,energyC,weightC,xlabel,ylabel,xscale,yscale,pbins,filename},If[Length[dat]<5,Abort[]];
energyC=2;
weightC=3;
{xC,xlabel,xscale,ylabel,yscale}=Switch[dimension,"energy",{energyC,"E (MeV)",(10^(-6)) (*eV/MeV*),"I/E (A/MeV)",(10^6)(*MeV/eV*)}];
pbins={xscale #1,yscale #2}&@@@BinParticles[xC,weightC,nbins,dat];
(*pbins=Select[pbins,#[[2]]>0&];*)
ListLogPlot[pbins,Filling->Axis,Frame->True,PlotRange->{All,All},BaseStyle->{FontSize->15},FrameLabel->{xlabel,ylabel},ImageSize->1600/Sqrt[2],FillingStyle->color,PlotStyle->Black,ImagePadding->1/Sqrt[2] {{80,60},{100,20}},Axes->False,AspectRatio->1/2,Joined->True,PlotRange->{All,{(10^(-24)),(10^(-14))}}]]

weightedE=BunchHistogram[testDat,"energy"]

Export[fileName<>"_weightedE.png",weightedE];

(* Weighted 2D Histogram of energy & s *)

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*)]

dat=BinParticles2D[{1,2,3},{250,250},testDat];

dat=Select[dat,#[[3]]≠0&];

ListDensityPlot[dat,PlotRange->All]