This notebook also can be downloaded or executed at MS Azure cloud: https://notebooks.azure.com/library/fnal (press 'Clone and run' button).
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 the function to read MAD-X twiss output
from urllib.request import urlopen
from io import StringIO
import shlex
def read_twiss(url):
txt = urlopen(url).read().decode()
# make stream out of string (it's easier to read then):
DATA=StringIO(txt)
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)
return df, twiss_pars
df_twiss, twiss_pars = read_twiss('https://apetrenko.blob.core.windows.net/fnal-booster-model/twiss.outx')
print(twiss_pars.TITLE + ' ' + twiss_pars.DATE)
# full list of parameters from the Twiss table
print(list(twiss_pars.index.values))
df_twiss.head()
# beta-functions:
p = figure(tools="save,box_zoom,pan,reset", toolbar_location="above",
logo="grey", plot_width=800, plot_height=300, active_drag="box_zoom")
p.line(df_twiss.S, df_twiss.BETX, color='red', line_width=2, line_alpha=0.5, legend='Beta x')
p.line(df_twiss.S, df_twiss.BETY, color='blue', line_width=2, line_alpha=0.5, legend='Beta y')
p.legend.background_fill_alpha = 0.5
p.y_range.start=0; p.y_range.end = df_twiss[['BETX','BETY']].values.max() + 3
p.xaxis.axis_label='s (m)'
p.yaxis.axis_label='Beta x,y (m)'
show(p)
# plotting layout of elements in twiss.outx file:
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)
# beta-functions:
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='Beta x')
p.line(df_twiss.S, df_twiss.BETY, color='blue', line_width=2, line_alpha=0.5, legend='Beta y')
p.legend.background_fill_alpha = 0.5
p.y_range.start=0; p.y_range.end = 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='Beta x,y (m)'
lay=gridplot([
[p],
[p_layout]
], toolbar_location='right')
show(lay)
# dispersion:
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=2, line_alpha=0.5, legend='DISPX')
p.line(df_twiss.S, df_twiss.DISPY, color='blue', line_width=2, line_alpha=0.5, legend='DISPY')
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.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)
To calculate the model orbit responses in the general case of coupled optics we need the transport matrix. The transport matrix is calculated from tracking results of particle with shifted initial coordinates.
MAD-X file generates 5 twiss.outx-like files with trajectories corresponding to different initial particle coordinates: r_dX.outx, r_dY.outx, r_dPX.outx, r_dPY.outx, r_dPT.outx.
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('https://apetrenko.blob.core.windows.net/fnal-booster-model/r_dX.outx')
df_dY, pars = read_twiss('https://apetrenko.blob.core.windows.net/fnal-booster-model/r_dY.outx')
df_dPX, pars = read_twiss('https://apetrenko.blob.core.windows.net/fnal-booster-model/r_dPX.outx')
df_dPY, pars = read_twiss('https://apetrenko.blob.core.windows.net/fnal-booster-model/r_dPY.outx')
df_dPT, pars = read_twiss('https://apetrenko.blob.core.windows.net/fnal-booster-model/r_dPT.outx')
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
df_twiss.tail()
Now let's calculate the betatron tunes from this transport matrix and compare them to the twiss_pars.Q1,Q2 values. 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)) )
The following function finds the orbit response of a coasting beam, without the effect of horizontal feedback (it will be included later):
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
dipole = 'HL11'
xresp, yresp = coasting_xy_response_to(dipole)
df_twiss['dx_d' + dipole + 'coast'] = xresp
df_twiss['dy_d' + dipole + 'coast'] = yresp
df_twiss.head()
# Plot the response to see if everything is fine:
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)
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
Let's compare the obtained Dx and Dy to the twiss.outx DISPX DISPY:
# dispersion:
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)
The RPOS feedback loop adjusts the energy of the beam so that it's location at the RPOS BPM does not change.
df_twiss[df_twiss['NAME']=='RPOS']
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'] * x_RPOS_resp / Dx_RPOS
y_term = - df_twiss['Dy'] * x_RPOS_resp / Dx_RPOS
return x_term, y_term
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)
# saved at the end of the ORM_viewer.ipynb:
url = 'https://apetrenko.blob.core.windows.net/fnal-booster-orm/2016-11-22_ORM_nb.txt'
df_ORM = pd.read_csv(url, delim_whitespace=True)
df_ORM.head()
Now we can calculate all model orbit responses
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()
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')]
According to J. DiMarco et al. Test Results of the AC Field Measurements of Fermilab Booster Corrector Magnets // 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
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)
It looks like the kick vs angle coefficient is underestimated significantly (by ~30%) (30% higher current is required to get the required kick). Let's try to rescale this coefficient:
A_per_mrad = 1.3 * 8.55*np.sqrt(GAMMA*GAMMA - 1)
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)
The agreement in both planes is decent. It looks like a good starting point for the LOCO fit. However it's important to find the source of the 30% corrector or BPM calibration error!