This program compares the measured Orbit Response and compares it to the MADX model. This program was originally written by A. Petrenko and modified by C.Y. Tan.
The required libraries are:
import numpy as np
import pandas as pd
from bokeh.plotting import figure, output_file, show
from bokeh.io import output_notebook
from bokeh.resources import INLINE
output_notebook(resources=INLINE)
Define functions to read in the twiss.out file generated by MADX
import shlex;
def read_twiss(fname):
DATA=open(fname, "r");
twiss_pars = pd.Series()
for line in DATA:
if line.startswith('@'):
itms = shlex.split(line);
twiss_pars[itms[1]] = itms[-1];
elif line.startswith('*'):
cols = shlex.split(line)[1:];
break;
DATA.readline();
df = pd.read_table(DATA, names=cols, delim_whitespace=True);
DATA.close();
return df, twiss_pars;
Set up directory where the data resides
dir="/Users/cytan/expt/booster/13Mar2017/check/";
df_twiss, twiss_pars = read_twiss(dir+"twiss.outx");
Display the first few lines of df
df_twiss.head()
and the twiss parameters
twiss_pars.index.values
First create the MADX elements graphics
from bokeh.models import (HoverTool, ColumnDataSource,
CustomJS)
from bokeh.layouts import gridplot
def pic_h(row):
if row.ANGLE != 0: return 6
if row.NAME.startswith('BPM'): return 4
return 3
def pic_L(row):
if row.L < 0.01:
return 0.01
else:
return row.L
def pic_color(row):
if row.ANGLE != 0:
if row.K1L > 0: return 'red'
return 'blue'
if row.NAME.startswith('BPM'): return 'green'
return 'white'
df = df_twiss
df["pic_h"] = df.apply(pic_h, axis = 1)
df["pic_L"] = df.apply(pic_L, axis = 1)
df["pic_color"] = df.apply(pic_color, axis = 1)
all_elements_src = ColumnDataSource(data=dict(
NAME=df_twiss.NAME.values,
s = df_twiss.S.values,
L = df_twiss.pic_L.values,
h = df_twiss.pic_h.values,
color = df_twiss.pic_color.values,
))
p = figure(tools="save,xbox_zoom,xpan,reset", toolbar_location="above",
logo="grey", plot_width=800, plot_height=100, active_drag="xbox_zoom",
y_axis_type=None)
p.x_range.start = df.S.values.min() - 7
p.x_range.end = df.S.values.max() + 7
r1 = p.rect('s', 0, width='L', height='h', fill_color='color',
fill_alpha=0.5, line_width=1.0, line_alpha=0.2,
line_color='black', source=all_elements_src)
tips = [
('Element', '@NAME'),
]
hover = HoverTool(tooltips=tips)
p.add_tools(hover)
#hover = p.select_one(HoverTool)
#hover.tooltips = tips
#hover.renderers=[c1,c2,r1,r2]
hover.renderers=[r1]
p.y_range.start = -7
p.y_range.end = +7
p.xaxis.axis_label='s (m)'
p_layout = p
#show(p_layout)
Then the twiss and dispersion functions
from bokeh.models import LinearAxis, Range1d;
p = figure(tools="save,xbox_zoom,xpan,reset", x_range=p_layout.x_range, toolbar_location="above",
logo="grey", plot_width=800, plot_height=300, active_drag="xbox_zoom");
p.line(df_twiss.S, df_twiss.BETX, color='red', line_width=2, line_alpha=0.5, legend='\u03b2x');
p.line(df_twiss.S, df_twiss.BETY, color='blue', line_width=2, line_alpha=0.5, legend='\u03b2y');
p.legend.background_fill_alpha = 0.5
p.y_range.start=-5; p.y_range.end = round(df_twiss[['BETX','BETY']].values.max() + 3);
#p.xaxis.axis_label='s (m)'
p.xaxis.major_label_text_alpha = 0;
p.xaxis.major_label_text_font_size = '1pt';
p.yaxis.axis_label='\u03b2 (m)';
p.extra_y_ranges={"disp": Range1d(start=0, end=15)};
p.add_layout(LinearAxis(y_range_name="disp", axis_label="Dispersion (m)"), "right");
p.line(df_twiss.S, df_twiss.DISPX, color='green', line_width=2, line_alpha=0.5, legend='Dx', y_range_name="disp");
lay=gridplot([
[p],
[p_layout]
], toolbar_location='right');
show(lay);
The orbit responses have to be calculated with transport matrices. The transport matrices are generated using the results from MADX from shifted initial coordinates.
The MADX generated files for the following initial particle coordinates:
dX = 1e-4; dPX = 0.5e-5; dY = 1e-4; dPY = 1e-5; dPT = 0.5e-4; # from the booster.madx
df_dX, pars = read_twiss(dir+"r_dX.outx");
df_dY, pars = read_twiss(dir+"r_dY.outx");
df_dPX, pars = read_twiss(dir+"r_dPX.outx");
df_dPY, pars = read_twiss(dir+"r_dPY.outx");
df_dPT, pars = read_twiss(dir+"r_dPT.outx");
Check
df_dPT.head()
df_twiss['R11'] = ( df_dX['X'] - df_twiss['X'] )/dX
df_twiss['R21'] = ( df_dX['PX'] - df_twiss['PX'])/dX
df_twiss['R31'] = ( df_dX['Y'] - df_twiss['Y'] )/dX
df_twiss['R41'] = ( df_dX['PY'] - df_twiss['PY'])/dX
df_twiss['R12'] = (df_dPX['X'] - df_twiss['X'] )/dPX
df_twiss['R22'] = (df_dPX['PX'] - df_twiss['PX'])/dPX
df_twiss['R32'] = (df_dPX['Y'] - df_twiss['Y'] )/dPX
df_twiss['R42'] = (df_dPX['PY'] - df_twiss['PY'])/dPX
df_twiss['R13'] = ( df_dY['X'] - df_twiss['X'] )/dY
df_twiss['R23'] = ( df_dY['PX'] - df_twiss['PX'])/dY
df_twiss['R33'] = ( df_dY['Y'] - df_twiss['Y'] )/dY
df_twiss['R43'] = ( df_dY['PY'] - df_twiss['PY'])/dY
df_twiss['R14'] = (df_dPY['X'] - df_twiss['X'] )/dPY
df_twiss['R24'] = (df_dPY['PX'] - df_twiss['PX'])/dPY
df_twiss['R34'] = (df_dPY['Y'] - df_twiss['Y'] )/dPY
df_twiss['R44'] = (df_dPY['PY'] - df_twiss['PY'])/dPY
GAMMA = float(twiss_pars['GAMMA'])
BETA = np.sqrt(1 - 1/(GAMMA*GAMMA))
df_twiss['R16'] = BETA*(df_dPT['X'] - df_twiss['X'] )/dPT
df_twiss['R26'] = BETA*(df_dPT['PX'] - df_twiss['PX'])/dPT
df_twiss['R36'] = BETA*(df_dPT['Y'] - df_twiss['Y'] )/dPT
df_twiss['R46'] = BETA*(df_dPT['PY'] - df_twiss['PY'])/dPT
Check
df_twiss.tail()
df_twiss.tail(1)
Bow, we can calculate the betatron tunes from the transport matrix and compare with the MADX calculated tunes: Q1, Q2 values. The last row of the df_twiss table gives us the full turn matrix.
from numpy import linalg as LA
[M11, M12, M13, M14,
M21, M22, M23, M24,
M31, M32, M33, M34,
M41, M42, M43, M44] = df_twiss.tail(1)[[
'R11', 'R12', 'R13', 'R14',
'R21', 'R22', 'R23', 'R24',
'R31', 'R32', 'R33', 'R34',
'R41', 'R42', 'R43', 'R44'
]].values[0]
M = np.array([
[M11, M12, M13, M14],
[M21, M22, M23, M24],
[M31, M32, M33, M34],
[M41, M42, M43, M44]
])
print('1-turn transport matrix M =')
print(M)
w, v = LA.eig(M)
tunes = 1 - np.angle(w)/(2*np.pi)
Q1 = tunes[0]
Q2 = tunes[2]
print( 'Coupled Q1 = {0:.7f}'.format(Q1) )
print( 'twiss.outx Q1 = {0:.7f}'.format(float(twiss_pars.Q1)) )
print( 'Coupled Q2 = {0:.7f}'.format(Q2))
print( 'twiss.outx Q2 = {0:.7f}'.format(float(twiss_pars.Q2)) )
First create the functions that finds the orbit response of coasting beam without RPOS feedback
def coasting_xy_response_to(dipole):
# Reusing the dX variable (we don't need dX from MAD-X any more):
if dipole.startswith('H'):
dX = np.array([
[0],
[1],
[0],
[0]
])
else:
dX = np.array([
[0],
[0],
[0],
[1]
])
# transport matrix from S=0 to corrector:
df_Rc = df_twiss[df_twiss.NAME == dipole]
[R11, R12, R13, R14,
R21, R22, R23, R24,
R31, R32, R33, R34,
R41, R42, R43, R44] = df_Rc[[
'R11', 'R12', 'R13', 'R14',
'R21', 'R22', 'R23', 'R24',
'R31', 'R32', 'R33', 'R34',
'R41', 'R42', 'R43', 'R44'
]].values[0]
Rc = np.array([
[R11, R12, R13, R14],
[R21, R22, R23, R24],
[R31, R32, R33, R34],
[R41, R42, R43, R44]
])
Sc = df_Rc.S.values[0]
# 1 turn matrix at the location of the corrector: Mc = Rc*M0*inv(Rc)
Mc = np.dot(Rc, np.dot(M, LA.inv(Rc)) )
#Xc = inv(eye(4)-Mc)*Mc*dX # solution to the equation Xc = Mc*(Xc + dX)
Xc = np.dot(LA.inv(np.eye(4)-Mc), np.dot(Mc, dX))
#X0 = inv(Rc)*Xc; # solution of Rc*X0 = Xc -- X at s = 0
X0 = np.dot(LA.inv(Rc), Xc)
dX0 = np.dot(LA.inv(Rc), dX)
# Now we can find orbit distortion around the ring:
s_c = df_Rc.S.values[0]
dx_dkick = []
dy_dkick = []
for i, row in df_twiss.iterrows():
[R11, R12, R13, R14,
R21, R22, R23, R24,
R31, R32, R33, R34,
R41, R42, R43, R44] = row[[
'R11', 'R12', 'R13', 'R14',
'R21', 'R22', 'R23', 'R24',
'R31', 'R32', 'R33', 'R34',
'R41', 'R42', 'R43', 'R44'
]].values
Rs = np.array([
[R11, R12, R13, R14],
[R21, R22, R23, R24],
[R31, R32, R33, R34],
[R41, R42, R43, R44]
])
#print('S = {0}'.format(row.S))
if row.S <= Sc:
# [x x' y y']' at the current location:
Xs = np.dot(Rs,X0)
else:
Xs = np.dot(Rs, X0 + dX0);
dx_dkick.append(Xs[0,0])
dy_dkick.append(Xs[2,0])
return dx_dkick, dy_dkick
Now using the remaining $R_{i6}$ elements we'll find the coupled dispersion function in $x$ and $y$ and take into account the RPOS feedback.
Dispersion function is defined as the closed beam orbit with initail coordinates $X_0 = (x, x', y, y', \Delta l, \Delta p/p)^T$ all satisfying the periodic condition except for $\Delta l$:
$$ \left( \begin{array}{c} x\\ x'\\ y\\ y'\\ \Delta l_1\\ \Delta p/p \end{array} \right) = \begin{pmatrix} M_{11} & M_{12} & M_{13} & M_{14} & 0 & M_{16} \\ M_{21} & M_{22} & M_{23} & M_{24} & 0 & M_{26} \\ M_{31} & M_{32} & M_{33} & M_{34} & 0 & M_{36} \\ M_{41} & M_{42} & M_{43} & M_{44} & 0 & M_{46} \\ M_{51} & M_{52} & M_{53} & M_{54} & 1 & M_{56} \\ 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} \left( \begin{array}{c} x\\ x'\\ y\\ y'\\ \Delta l\\ \Delta p/p \end{array} \right) $$or
$$ \left( \begin{array}{c} x\\ x'\\ y\\ y'\end{array} \right) = \begin{pmatrix} M_{11} & M_{12} & M_{13} & M_{14} \\ M_{21} & M_{22} & M_{23} & M_{24} \\ M_{31} & M_{32} & M_{33} & M_{34} \\ M_{41} & M_{42} & M_{43} & M_{44} \end{pmatrix} \left( \begin{array}{c} x\\ x'\\ y\\ y'\end{array} \right) + \begin{pmatrix} M_{16} \\ M_{26} \\ M_{36} \\ M_{46} \\ \end{pmatrix} \frac{\Delta p}{p}. $$This matrix equation $X_0 = M X_0 + M_{i6} \Delta p/p$ has the following solution: $$ X_0 = (I - M)^{-1} M_{i6} \frac{\Delta p}{p}. $$
Hence the vector dispersion function can be found from propagating the $X_0$ further along the ring:
$$ X(s) = R(s) X_0 + R_{i6} \frac{\Delta p}{p}. $$[M16,
M26,
M36,
M46] = df_twiss.tail(1)[[
'R16',
'R26',
'R36',
'R46']].values[0]
Mi6 = np.array([
[M16],
[M26],
[M36],
[M46]
])
X0 = np.dot( LA.inv(np.eye(4)-M), Mi6) # solution to the equation X0 = M*X0 + Mi6*dp/p, where dp/p = 1
#print(M)
#print(X0)
#np.dot(M, X0) + Mi6
Dx = []
Dy = []
for i, row in df_twiss.iterrows():
[R11, R12, R13, R14, R16,
R21, R22, R23, R24, R26,
R31, R32, R33, R34, R36,
R41, R42, R43, R44, R46] = row[[
'R11', 'R12', 'R13', 'R14', 'R16',
'R21', 'R22', 'R23', 'R24', 'R26',
'R31', 'R32', 'R33', 'R34', 'R36',
'R41', 'R42', 'R43', 'R44', 'R46'
]].values
Rs = np.array([
[R11, R12, R13, R14],
[R21, R22, R23, R24],
[R31, R32, R33, R34],
[R41, R42, R43, R44]
])
Ri6 = np.array([
[R16],
[R26],
[R36],
[R46]
])
Xs = np.dot(Rs, X0) + Ri6
Dx.append(Xs[0,0])
Dy.append(Xs[2,0])
df_twiss['Dx'] = Dx
df_twiss['Dy'] = Dy
p = figure(tools="save,box_zoom,xpan,reset", toolbar_location="above", x_range=p_layout.x_range,
logo="grey", plot_width=800, plot_height=300, active_drag="box_zoom")
p.line(df_twiss.S, df_twiss.DISPX, color='red', line_width=6, line_alpha=0.5, legend='DISPX')
p.line(df_twiss.S, df_twiss.DISPY, color='blue', line_width=6, line_alpha=0.5, legend='DISPY')
p.line(df_twiss.S, df_twiss.Dx, color='black', line_width=2, line_alpha=0.7, legend='Dx')
p.line(df_twiss.S, df_twiss.Dy, color='black', line_width=2, line_alpha=0.7, legend='Dy', line_dash = [5, 2])
p.legend.background_fill_alpha = 0.5
#p.y_range.start = df_twiss[['DISPX','DISPY']].values.min() - 0.2
#p.y_range.end = df_twiss[['DISPX','DISPY']].values.max() + 0.5
#p.xaxis.axis_label='s (m)'
p.yaxis.axis_label='DISPX, DISPY (m)'
#p.xaxis.axis_label='s (m)'
p.xaxis.major_label_text_alpha = 0
p.xaxis.major_label_text_font_size = '1pt'
lay=gridplot([
[p],
[p_layout]
], toolbar_location='right')
show(lay)
Now we can add in RPOS feedback. The RPOS feedback loop adjusts the energy of the beam so that the beam at the RPOS bpm does not change.
I have to multiply by rfactor for the correction to rpos or else the response offset is wrong. This is very peculiar and needs to be investigated further
rfactor = 1.0;
def RPOS_term_of_xy_response_to(dipole):
df = df_twiss[df_twiss['NAME']=='RPOS']
Dx_RPOS = df.Dx.values[0]
x_RPOS_resp = df['dx_d'+dipole + 'coast'].values[0]
# x_RPOS_resp * kick + Dx_RPOS*delta = 0 =>
# delta/kick = - x_RPOS_resp / Dx_RPOS
# dx/dkick = dx/dkick (coast) - x_RPOS_resp Dx / Dx_RPOS
x_term = - df_twiss['Dx'] * rfactor*x_RPOS_resp / Dx_RPOS
y_term = - df_twiss['Dy'] * rfactor*x_RPOS_resp / Dx_RPOS
return x_term, y_term
dipole = "HL02";
xresp, yresp = coasting_xy_response_to(dipole);
df_twiss['dx_d' + dipole + 'coast'] = xresp;
df_twiss['dy_d' + dipole + 'coast'] = yresp;
p = figure(tools="save,box_zoom,pan,reset", toolbar_location="above", x_range=p_layout.x_range,
logo="grey", plot_width=800, plot_height=300, active_drag="box_zoom")
p.line(df_twiss.S, df_twiss['dx_d' + dipole + 'coast'], color='red', line_width=2, line_alpha=0.5, legend='dx/dkick')
p.line(df_twiss.S, df_twiss['dy_d' + dipole + 'coast'], color='blue', line_width=2, line_alpha=0.5, legend='dy/dkick')
p.legend.background_fill_alpha = 0.5
#p.xaxis.axis_label='s (m)'
p.title.text = 'Example orbit response for the coasting beam'
p.yaxis.axis_label='dx,y/dkick (mm/mrad)'
#p.xaxis.axis_label='s (m)'
p.xaxis.major_label_text_alpha = 0
p.xaxis.major_label_text_font_size = '1pt'
lay=gridplot([
[p],
[p_layout]
], toolbar_location='right')
show(lay)
dipole = "HL02";
x_term, y_term = RPOS_term_of_xy_response_to(dipole);
df_twiss['dx_d' + dipole] = df_twiss['dx_d' + dipole + 'coast'] + x_term;
df_twiss['dy_d' + dipole] = df_twiss['dy_d' + dipole + 'coast'] + y_term;
from bokeh.models import LabelSet
# Plot the response to see if everything is fine:
p = figure(tools="save,box_zoom,xpan,reset", toolbar_location="above", x_range=p_layout.x_range,
logo="grey", plot_width=800, plot_height=300, active_drag="box_zoom")
p.line(df_twiss.S, df_twiss['dx_d' + dipole], color='red', line_width=2, line_alpha=0.5, legend='dx/dkick')
p.line(df_twiss.S, df_twiss['dy_d' + dipole], color='blue', line_width=2, line_alpha=0.5, legend='dy/dkick')
s_RPOS = df_twiss[df_twiss['NAME']=='RPOS'].S.values[0]
p.line([s_RPOS, s_RPOS], [-1, +1], color='black', line_width=1, line_alpha=0.7)
RPOS_src = ColumnDataSource(data=dict(
NAME=['RPOS'],
s = [s_RPOS],
))
lbl = LabelSet(x='s', y=+1, text='NAME', x_offset=0, y_offset=0,
text_font_size="10pt", text_color="black",
text_align='center', text_baseline='bottom', source=RPOS_src)
p.add_layout(lbl)
p.legend.background_fill_alpha = 0.5
#p.xaxis.axis_label='s (m)'
p.title.text = 'Example orbit response to ' + dipole + ' corrector'
p.yaxis.axis_label='dx,y/dkick (mm/mrad)'
#p.xaxis.axis_label='s (m)'
p.xaxis.major_label_text_alpha = 0
p.xaxis.major_label_text_font_size = '1pt'
lay=gridplot([
[p],
[p_layout]
], toolbar_location='right')
show(lay)
The measured ORM is compared to the model. This requires that ORM_viewer.ipynb be run first to create the required slope and error files
dir_m = "/Users/cytan/expt/booster/13Mar2017/"
df_ORM = pd.read_csv(dir_m+"13Mar2017_ORM_nb.txt", delim_whitespace=True )
df_ORM.head()
Now we can calculate the model orbit response
for dipole in df_ORM.Dipole.values:
print('\rWorking on ' + dipole + '...', end='')
xresp, yresp = coasting_xy_response_to(dipole)
df_twiss['dx_d' + dipole + 'coast'] = xresp
df_twiss['dy_d' + dipole + 'coast'] = yresp
x_term, y_term = RPOS_term_of_xy_response_to(dipole)
df_twiss['dx_d' + dipole] = df_twiss['dx_d' + dipole + 'coast'] + x_term
df_twiss['dy_d' + dipole] = df_twiss['dy_d' + dipole + 'coast'] + y_term
print('\nDone.')
df_twiss.head()
According to J. DiMarco et al. <a href=https://accelconf.web.cern.ch/accelconf/e08/papers/wepc144.pdf>Test Results of the AC Field Measurements of Fermilab Booster Corrector Magnets</a> // Proceedings of EPAC'2008 the corrector strength is $$ \frac{B L} {I} = 0.000366~\mathrm{\frac{T \cdot m}{A}}. $$ Angle given to the beam by the corrector dipole is $$ \theta = \frac{B L}{B \rho}, $$ where $B \rho$ is the rigidity proportional to the particle momentum. Let's calculate $B \rho$ from proton $\gamma$ which is available as twiss_pars.GAMMA (and we've used it already in dispersion calculations). In CGS units $$ \rho = \frac{pc}{eB} \Rightarrow $$ $$ B\rho~[\mathrm{Gs \cdot cm}] = \frac{pc}{e} = \frac{\gamma m_p v c} {e} = \gamma \beta \frac{m_p c^2} {e}. $$ Since $m_p c^2 \approx 0.938~\mathrm{GeV}$ and $\beta = \sqrt{1-1/\gamma^{2}}$, the rigidity can be expressed as $$ B\rho~[\mathrm{Gs \cdot cm}] = \sqrt{\gamma^2-1} \frac{0.938\cdot 10^{9} e \frac{1}{299.8} \mathrm{stat V}} {e} = 3.13\cdot 10^{6} \sqrt{\gamma^2-1}~[\mathrm{Gs \cdot cm}]. $$ or $$ B\rho~[\mathrm{T \cdot m}] = 3.13\cdot 10^{6} \sqrt{\gamma^2-1}~[\mathrm{10^{-4}~T \cdot 10^{-2}~m}] = 3.13 \sqrt{\gamma^2-1}~[\mathrm{T \cdot m}]. $$ The final expression for the corrector kick vs current: $$ \frac{\theta}{I} \left [ \mathrm{\frac{mrad}{A}} \right ] = \frac{0.366} {3.13\sqrt{\gamma^2-1}} = \frac{1} {8.55\sqrt{\gamma^2-1}}. $$
GAMMA = float(twiss_pars.GAMMA)
print('MAD-X Twiss file GAMMA = {0:.2f} => proton kinetic energy = {1:.0f} MeV'.format(GAMMA, GAMMA*938 - 938))
A_per_mrad = 8.55*np.sqrt(GAMMA*GAMMA - 1) # Ampere per mrad
Setup for the plot
from bokeh.models.widgets import Panel, Tabs, RadioGroup
from bokeh.models import LabelSet
cols = df_ORM.columns.values
BPMs = [col[:-5] for col in cols if col.endswith('slope')]
xBPMs = [itm for itm in BPMs if itm.startswith('H')]
yBPMs = [itm for itm in BPMs if itm.startswith('V')]
BPMs = [itm[1:] for itm in xBPMs] # names of BPMs without H,V.
s_BPMs = [df_twiss[df_twiss.NAME == 'BPM'+itm].S.values[0] for itm in BPMs]
Dipoles = df_ORM.Dipole.values
xDipoles = [itm for itm in Dipoles if itm.startswith('H')]
yDipoles = [itm for itm in Dipoles if itm.startswith('V')]
def plot_responses(dipole_names, ylabel):
dipole_s_list = []
dipole_list = []
text_y_list = []
x_resp_list = []
x_resp_model_list = []
x_err_list = []
y_resp_list = []
y_resp_model_list = []
y_err_list = []
x_slope_cols = [bpm+'slope' for bpm in xBPMs]
x_err_cols = [bpm+'std_err' for bpm in xBPMs]
y_slope_cols = [bpm+'slope' for bpm in yBPMs]
y_err_cols = [bpm+'std_err' for bpm in yBPMs]
for dipole in dipole_names:
df = df_ORM[df_ORM.Dipole == dipole]
dipole_list.append(dipole[1:])
dipole_s_list.append(df.s.values[0])
if dipole[1] =='S':
text_y_list.append(+1.9)
else:
text_y_list.append(-1.9)
x_slopes = list(df[x_slope_cols].values[0]*A_per_mrad)
x_slopes = [float(np.round(val, 2)) for val in x_slopes]
x_resp_list.append(x_slopes)
x_slopes_model = list(df_twiss['dx_d'+dipole].values)
x_slopes_model = [float(np.round(val, 2)) for val in x_slopes_model]
x_resp_model_list.append(x_slopes_model)
x_errs = list(df[x_err_cols].values[0]*A_per_mrad)
x_errs = [float(np.round(val*4, 3)) for val in x_errs]
x_err_list.append(x_errs)
y_slopes = list(df[y_slope_cols].values[0]*A_per_mrad)
y_slopes = [float(np.round(val, 2)) for val in y_slopes]
y_resp_list.append(y_slopes)
y_slopes_model = list(df_twiss['dy_d'+dipole].values)
y_slopes_model = [float(np.round(val, 2)) for val in y_slopes_model]
y_resp_model_list.append(y_slopes_model)
y_errs = list(df[y_err_cols].values[0]*A_per_mrad)
y_errs = [float(np.round(val*4, 3)) for val in y_errs]
y_err_list.append(y_errs)
v = np.array(dipole_s_list)*0
orbit_responses_src = ColumnDataSource(data=dict(
dipole=dipole_list,
dipole_s=dipole_s_list,
text_y=text_y_list,
x_resp=x_resp_list,
x_resp_model=x_resp_model_list,
x_err=x_err_list,
y_resp=y_resp_list,
y_resp_model=y_resp_model_list,
y_err=y_err_list,
))
#print(orbit_responses_src.data)
selected_dipole_src = ColumnDataSource(data=dict(
dipole_index = [0],
dipole_s = [float('nan')],
text_y = [float('nan')],
))
s = np.array(s_BPMs)
v = s*0
selected_orbit_response_src = ColumnDataSource(data=dict(
bpm_s = list(s),
bpm_name = list(BPMs),
x_response = x_slopes,
x_response_err = x_errs,
y_response = y_slopes,
y_response_err = y_errs,
))
s = np.array(df_twiss.S.values)
v = s*0
selected_orbit_response_model_src = ColumnDataSource(data=dict(
s = list(s),
x_response = x_slopes_model,
y_response = y_slopes_model,
))
p = figure(tools="xbox_zoom,pan,reset", toolbar_location="above", y_axis_type=None, #x_axis_type=None,
logo="grey", plot_width=800, plot_height=100, active_drag="xbox_zoom")
p.grid.grid_line_color = None
p.axis.axis_line_color = None
#p.axis.major_tick_line_color = None
#p.yaxis.major_tick_line_color = None
p.x_range.start=-15; p.x_range.end = max(s)+15
p.y_range.start=-4.5; p.y_range.end = +4.5
p.rect("dipole_s", 'text_y', width=18, height=3, source=selected_dipole_src,
line_color="black", line_alpha=1.0, fill_alpha=0.4)
r = p.rect("dipole_s", 'text_y', width=18, height=3, source=orbit_responses_src,
#hover_color='red', hover_alpha=0.4,
line_color="black", line_alpha=0.2, fill_alpha=0.2)
#p.line([p.x_range.start, p.x_range.end], 0, line_color="black", line_alpha=0.3)
labels = LabelSet(x="dipole_s", y='text_y', text="dipole", x_offset=0, y_offset=0,
text_font_size="9pt", text_color="black",
source=orbit_responses_src, text_align='center', text_baseline='middle')
p.add_layout(labels)
hover_callback = CustomJS(args={
'src': orbit_responses_src,
'dipole_src': selected_dipole_src,
's_src': selected_orbit_response_src,
'm_src': selected_orbit_response_model_src,
},
code='''
var indices = cb_data.index['1d'].indices;
if(indices.length > 0){
i = indices[0];
if(i != dipole_src.data.dipole_index[0]){
dipole_src.data.dipole_index[0] = i;
data = src.data;
dipole = data.dipole[i];
s = data.dipole_s[i];
y = data.text_y[i];
dipole_src.data.dipole_s[0] = s;
dipole_src.data.text_y[0] = y;
s_src.data.x_response = src.data.x_resp[i];
s_src.data.x_response_err = src.data.x_err[i];
s_src.data.y_response = src.data.y_resp[i];
s_src.data.y_response_err = src.data.y_err[i];
m_src.data.x_response = src.data.x_resp_model[i];
m_src.data.y_response = src.data.y_resp_model[i];
//console.log(dipole);
s_src.trigger('change');
m_src.trigger('change');
dipole_src.trigger('change');
}
}
''')
p.add_tools(HoverTool(tooltips=None, callback=hover_callback, renderers=[r]))
p_dipoles = p
p = figure(tools="save,box_zoom,pan,reset,hover", toolbar_location="above", x_range=p_layout.x_range,
logo="grey", plot_width=800, plot_height=400, active_drag="box_zoom")
p.x_range.start=-15; p.x_range.end = max(s)+15
p.y_range.start=-60; p.y_range.end = +60
#p.xaxis.axis_label='s (m)'
p.yaxis.axis_label=ylabel
p.xaxis.major_label_text_alpha = 0
p.xaxis.major_label_text_font_size = '1pt'
p.line('s', 'x_response', source=selected_orbit_response_model_src, color='red',
line_width=2, line_alpha=0.5, legend='dx/dkick')
c1 = p.circle('bpm_s', 'x_response', source=selected_orbit_response_src, color='red',
legend='dx/dkick')
r1 = p.rect("bpm_s", 'x_response', width=1.5, height='x_response_err', source=selected_orbit_response_src,
line_color=None, fill_color='red', fill_alpha=0.3)
p.line('s', 'y_response', source=selected_orbit_response_model_src, color='blue',
line_width=2, line_alpha=0.5, legend='dy/dkick')
c2 = p.circle('bpm_s', 'y_response', source=selected_orbit_response_src, color='blue',
legend='dy/dkick')
r2 = p.rect("bpm_s", 'y_response', width=1.5, height='y_response_err', source=selected_orbit_response_src,
line_color=None, fill_color='blue', fill_alpha=0.3)
tips = [
('BPM', '@bpm_name'),
]
hover = p.select_one(HoverTool)
hover.tooltips = tips
hover.renderers=[c1,c2,r1,r2]
#p.add_tools(HoverTool(tooltips=tips, renderers=[c1,c2]))
p.legend.background_fill_alpha = 0.4
lay=gridplot([
[p_dipoles],
[p],
[p_layout]
], toolbar_location='right')
return lay
lay = plot_responses(xDipoles, 'BPM response to x-kick (mm/mrad)')
show(lay)
lay = plot_responses(yDipoles, 'BPM response to y-kick (mm/mrad)')
show(lay)