Tracking¶
A. Petrenko (Novosibirsk, 2019)
In [1]:
import numpy as np
import pandas as pd
import holoviews as hv
import os
from IPython.display import Latex
from IPython.display import Image
hv.extension("bokeh")
In [2]:
%opts Curve Spread [width=700 height=300 show_grid=True \
default_tools=['box_zoom','pan','wheel_zoom','reset']]
Elegant lattice file:
Plot sdds-beam¶
In [56]:
w = "../debuncher/results/w3.sdds"
In [58]:
png = "results/image.png"
!sddsplot -lay=2,2 -device=lpng -output=$png \
-col=x,y -graph=dot,type=2 $w -endPanel \
-col=yp,y -graph=dot,type=2 $w -endPanel \
-col=x,xp -graph=dot,type=2 $w -endPanel \
-col=yp,xp -graph=dot,type=2 $w
Image(png)
Out[58]:
In [59]:
!sddsplot -device=lpng -output=$png -lay=2,2 "-title=" \
-col=t,p -factor=yMult=0.511,xMult=-3e11 \
"-xlabel=z (mm)" "-ylabel=p (MeV/c)" -graph=dot,type=2 $w -endPanel \
-col=xp,p -factor=yMult=0.511,xMult=1e3 \
"-xlabel=xp (mrad)" "-ylabel=p (MeV/c)" -graph=dot,type=2 $w -endPanel \
-col=t,x -factor=xMult=-3e11,yMult=1e3 \
"-xlabel=z (mm)" "-ylabel=x (mm)" -graph=dot,type=2 $w -endPanel \
-col=xp,x -factor=yMult=1e3,xMult=1e3 \
"-xlabel=xp (mrad)" "-ylabel=x (mm)" -graph=dot,type=2 $w -endPanel
Image(png)
Out[59]:
Run sddsanalyzebeam to check the beam parameters:
In [60]:
with open("machine.lte", "r") as f:
print(f.read())
In [61]:
with open("track.ele", "r") as f:
print(f.read())
Run Elegant
In [9]:
#!elegant track.ele
In [69]:
import subprocess
import sys
def run_cmd(cmd):
sub_process = subprocess.Popen(cmd, shell=True,
stdout=subprocess.PIPE, stderr=subprocess.PIPE)
while sub_process.poll() is None:
out = sub_process.stdout.read(10)
sys.stdout.write(out)
sys.stdout.flush()
In [71]:
#run_cmd("elegant track.ele")
#or run "mpirun -np 36 Pelegant track.ele"
In [83]:
turn = 500
w = "results/page.sdds"
In [84]:
!sddsconvert "results/w0.sdds" $w -keepPages=$turn
In [85]:
!sddsplot -lay=2,2 -device=lpng -output=$png \
-col=x,y -graph=dot,type=2 $w -endPanel \
-col=yp,y -graph=dot,type=2 $w -endPanel \
-col=x,xp -graph=dot,type=2 $w -endPanel \
-col=yp,xp -graph=dot,type=2 $w
Image(png)
Out[85]:
In [86]:
!sddsplot -device=lpng -output=$png -lay=2,2 "-title=" \
-col=dt,p -factor=yMult=0.511,xMult=-3e11 \
"-xlabel=z (mm)" "-ylabel=p (MeV/c)" -graph=dot,type=2 $w -endPanel \
-col=xp,p -factor=yMult=0.511,xMult=1e3 \
"-xlabel=xp (mrad)" "-ylabel=p (MeV/c)" -graph=dot,type=2 $w -endPanel \
-col=dt,x -factor=xMult=-3e11,yMult=1e3 \
"-xlabel=z (mm)" "-ylabel=x (mm)" -graph=dot,type=2 $w -endPanel \
-col=xp,x -factor=yMult=1e3,xMult=1e3 \
"-xlabel=xp (mrad)" "-ylabel=x (mm)" -graph=dot,type=2 $w -endPanel
Image(png)
Out[86]:
In [87]:
!sddsconvert "results/w3.sdds" $w -keepPages=$turn
In [88]:
!sddsplot -lay=2,2 -device=lpng -output=$png \
-col=x,y -graph=dot,type=2 $w -endPanel \
-col=yp,y -graph=dot,type=2 $w -endPanel \
-col=x,xp -graph=dot,type=2 $w -endPanel \
-col=yp,xp -graph=dot,type=2 $w
Image(png)
Out[88]:
In [89]:
!sddsplot -device=lpng -output=$png -lay=2,2 "-title=" \
-col=dt,p -factor=yMult=0.511,xMult=-3e11 \
"-xlabel=z (mm)" "-ylabel=p (MeV/c)" -graph=dot,type=2 $w -endPanel \
-col=xp,p -factor=yMult=0.511,xMult=1e3 \
"-xlabel=xp (mrad)" "-ylabel=p (MeV/c)" -graph=dot,type=2 $w -endPanel \
-col=dt,x -factor=xMult=-3e11,yMult=1e3 \
"-xlabel=z (mm)" "-ylabel=x (mm)" -graph=dot,type=2 $w -endPanel \
-col=xp,x -factor=yMult=1e3,xMult=1e3 \
"-xlabel=xp (mrad)" "-ylabel=x (mm)" -graph=dot,type=2 $w -endPanel
Image(png)
Out[89]:
Plot beam charge vs turn:
In [90]:
!sddsplot -device=lpng -output=$png "-title=" -scale=0,0,0,4e-9 \
-par=Pass,Charge -graph=line,type=1,thick=3 "results/w3.sdds"
Image(png)
Out[90]:
Now let's plot the results with holoviews¶
In [91]:
os.listdir("results")
Out[91]:
Reading beamline.mag
In [92]:
out = !sdds2stream results/beamline.mag -col=ElementName,s,Profile -pipe=out
In [93]:
out[:3]
Out[93]:
In [94]:
print("\n".join(out[:3]))
In [95]:
from io import StringIO
In [96]:
DATA = StringIO("\n".join(out))
In [97]:
df = pd.read_csv(DATA, names=['ElementName', 's', 'Profile'], delim_whitespace=True)
In [98]:
df.head(3)
Out[98]:
In [99]:
df.tail(3)
Out[99]:
In [100]:
dim_s = hv.Dimension('s', unit='m', label="s") # range=(-80,0))
In [101]:
#hv.help(hv.Curve)
In [102]:
from bokeh.models import HoverTool
#hover = HoverTool(tooltips=[("Name", "@ElementName")])
hover = HoverTool(tooltips="@ElementName")
%opts Curve.mag [height=70 show_frame=False show_title=False \
xaxis=None yaxis=None tools=['xbox_zoom, xpan', hover]] (color='black', alpha=0.3)
mag = hv.Curve(df, kdims=dim_s, vdims=['Profile', 'ElementName'], group='mag')
mag
Out[102]:
Reading track.sig file
In [103]:
#!sddsquery results/track.sig
In [104]:
out = !sdds2stream results/track.sig -col=ElementName,s,minimum1,maximum1,minimum3,maximum3,Sx,Sy -pipe=out
DATA = StringIO("\n".join(out))
df = pd.read_csv(DATA, names=['ElementName','s','xmin','xmax','ymin','ymax','sigma_x','sigma_y'],
delim_whitespace=True)
In [105]:
df.tail(3)
Out[105]:
In [106]:
df['xmin'] = df['xmin']*1e3; df['xmax'] = df['xmax']*1e3 # mm
df['ymin'] = df['ymin']*1e3; df['ymax'] = df['ymax']*1e3 # mm
df['sigma_x'] = df['sigma_x']*1e3; df['sigma_y'] = df['sigma_y']*1e3 # mm
In [107]:
#df.tail(3)
In [108]:
dim_x = hv.Dimension('x', unit='mm', label="x", range=(0,None))
dim_y = hv.Dimension('y', unit='mm', label="y", range=(0,None))
%opts Curve.fx (color='red', alpha=0.7, line_width=3)
%opts Curve.fy (color='blue', alpha=0.7, line_width=3)
In [109]:
Sx = hv.Curve((df.s, df.sigma_x), label='σx', kdims=dim_s, vdims=dim_x, group='fx')
Sy = hv.Curve((df.s, df.sigma_y), label='σy', kdims=dim_s, vdims=dim_y, group='fy')
In [110]:
(Sx*Sy + mag).cols(1)
Out[110]:
Reading track.cen file
In [111]:
#!sddsquery results/track.cen
In [112]:
out = !sdds2stream results/track.cen -col=ElementName,s,Particles,pCentral,Cx,Cy,Charge -pipe=out
DATA = StringIO("\n".join(out))
df_cen = pd.read_csv(DATA, names=['ElementName','s','Particles','p','Cx','Cy','Charge'],
delim_whitespace=True)
df_cen['p'] = 0.511*df_cen['p'] # MeV/c
df_cen['Cx'] = 1e3*df_cen['Cx'] # mm
df_cen['Cy'] = 1e3*df_cen['Cy'] # mm
In [113]:
df_cen.tail(3)
Out[113]:
In [114]:
dim_Q = hv.Dimension('Q', range=(0,None), unit='nC') # range=(-80,0))
Charge = hv.Curve((df_cen.s, df_cen.Charge*1e9), kdims=dim_s, vdims=dim_Q)
In [115]:
(Charge + mag).cols(1)
Out[115]:
Reading twiss.sdds file
In [116]:
out = !sdds2stream results/twiss.twi \
-col=ElementName,s,betax,betay,alphax,alphay,etax,etay,pCentral0,xAperture,yAperture -pipe=out
DATA = StringIO("\n".join(out))
df_twi = pd.read_csv(DATA, names=['ElementName','s','betax','betay','alphax','alphay',
'etax','etay','p','xAperture','yAperture'],
delim_whitespace=True)
In [117]:
df_twi['xAperture'] = 1e3*df_twi['xAperture'] # mm
df_twi['yAperture'] = 1e3*df_twi['yAperture'] # mm
In [118]:
df_twi.tail(3)
Out[118]:
In [119]:
dim_x = hv.Dimension('x', unit='mm', label="x", range=(-16,+16))
dim_y = hv.Dimension('y', unit='mm', label="y", range=(-16,+16))
dim_p = hv.Dimension('p', unit='MeV/c', label="p", range=(0,None))
In [120]:
%opts Curve.Aper (color='gray', alpha=0.5, line_width=3)
dim_sigma = hv.Dimension('sigma', unit='mm', label="2σ")#, range=(0,None))
dim_aper = hv.Dimension('aper', unit='mm', label="Aperture")#, range=(0,None))
Ax = hv.Curve((df_twi.s, df_twi.xAperture), label='Aperture', kdims=dim_s, vdims=dim_sigma, group='Aper')
Ax_pos = Ax
Ax_neg = hv.Curve((df_twi.s, -df_twi.xAperture), label='Aperture', kdims=dim_s, vdims=dim_sigma, group='Aper')
In [121]:
#Ax_pos*Ax_neg
In [122]:
#s_p = hv.Curve((df_cen.s, df_cen.p), kdims=dim_s, vdims=dim_p)
Cx = hv.Spread((df_cen.s, df_cen.Cx, 2*df.sigma_x), label='x', \
kdims=dim_s, vdims=[dim_x,dim_aper])
Cy = hv.Spread((df_cen.s, df_cen.Cy, 2*df.sigma_y), label='y', \
kdims=dim_s, vdims=[dim_y,dim_aper])
In [123]:
(Ax_pos*Ax_neg*Cy*Cx + mag).cols(1)
Out[123]:
In [124]:
s_Cx = hv.Curve((df_cen.s, df_cen.Cx), kdims=dim_s, vdims=dim_x, group='fx', label='x')
s_Cy = hv.Curve((df_cen.s, df_cen.Cy), kdims=dim_s, vdims=dim_x, group='fy', label='y')
In [125]:
(s_Cx*s_Cy*Ax_pos*Ax_neg + mag).cols(1)
Out[125]:
Resulting beam paramaters:
In [126]:
beam_analyzed = 'results/beam_analyzed.sdds'
#w = "results/w3.sdds"
!sddsanalyzebeam $w $beam_analyzed
In [127]:
items = np.array([
# name in file symbol factor units
['enx', '\\epsilon_{nx}', 1.0e6, 'mm \\cdot mrad'],
['betax', '\\beta_x', 1.0, 'm' ],
['alphax', '\\alpha_x', 1.0, '' ],
['Sx', '\\sigma_x', 1.0e3, 'mm' ],
['Sxp', "\\sigma_{x'}", 1.0e3, 'mrad' ],
['eny', '\\epsilon_{ny}', 1.0e6, 'mm \\cdot mrad'],
['betay', '\\beta_y', 1.0, 'm' ],
['alphay', '\\alpha_y', 1.0, '' ],
# ['Sy', '\\sigma_y', 1.0e3, 'mm' ],
# ['Syp', "\\sigma_{y'}", 1.0e3, 'mrad' ],
['St', '\\sigma_t', 1e12, 'ps' ],
['St', '\\sigma_z', 3e11, 'mm' ],
['Sdelta', '\\sigma_{\\Delta p/p}', 100.0, '\\%' ],
['pAverage', 'p_{average}', 0.511, 'MeV/c' ],
['Charge', 'Q', 1.0e9, 'nC' ],
])
names = items[:,0]
cols_str = "-col=" + ",".join(names)
cols_str
Out[127]:
In [128]:
out = !sdds2stream $beam_analyzed $cols_str -pipe=out
values = out[0].split()
values = np.array(values)
values = values.astype(np.float)
values = dict(zip(names, values))
In [129]:
latex_str = r"\begin{aligned}"
for name,symbol,factor,units in items:
value = values[name]*factor.astype(float)
latex_str = latex_str + r"%s &= %.3f~\rm{%s} \\" % (symbol, value, units)
latex_str = latex_str + r"\end{aligned}"
In [130]:
from IPython.display import Latex
In [131]:
Latex(latex_str)
Out[131]:
In [ ]: