Comparing ORM vs Model

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.

Input required python modules

The required libraries are:

  • numpy scientific computing library
  • pandas data analysis library
  • bokeh plotting library
In [1]:
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)
Loading BokehJS ...

Read in the model data

Define functions to read in the twiss.out file generated by MADX

In [2]:
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

In [3]:
dir="/Users/cytan/expt/booster/13Mar2017/check/";
In [4]:
df_twiss, twiss_pars = read_twiss(dir+"twiss.outx");

Display the first few lines of df

In [5]:
df_twiss.head()
Out[5]:
NAME S L TILT ANGLE X PX Y PY PT ... K1SL K2SL BETX BETY ALFX ALFY MUX MUY DISPX DISPY
0 BOOSTER$START 0.000 0.000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 33.316890 5.239851 0.074862 -0.033595 0.000000 0.000000 3.243458 0.005476
1 START_1 0.000 0.000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 33.316890 5.239851 0.074862 -0.033595 0.000000 0.000000 3.243458 0.005476
2 SA 0.088 0.176 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 33.303948 5.247243 0.072206 -0.050408 0.000420 0.002671 3.243285 0.005586
3 HS24 0.188 0.024 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 33.289808 5.259235 0.069188 -0.069514 0.000898 0.005701 3.243089 0.005712
4 VS24 0.212 0.024 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 33.286505 5.262682 0.068464 -0.074099 0.001013 0.006427 3.243042 0.005742

5 rows × 22 columns

and the twiss parameters

In [6]:
twiss_pars.index.values
Out[6]:
array(['NAME', 'TYPE', 'SEQUENCE', 'PARTICLE', 'MASS', 'CHARGE', 'ENERGY',
       'PC', 'GAMMA', 'KBUNCH', 'BCURRENT', 'SIGE', 'SIGT', 'NPART', 'EX',
       'EY', 'ET', 'LENGTH', 'ALFA', 'ORBIT5', 'GAMMATR', 'Q1', 'Q2',
       'DQ1', 'DQ2', 'DXMAX', 'DYMAX', 'XCOMAX', 'YCOMAX', 'BETXMAX',
       'BETYMAX', 'XCORMS', 'YCORMS', 'DXRMS', 'DYRMS', 'DELTAP',
       'SYNCH_1', 'SYNCH_2', 'SYNCH_3', 'SYNCH_4', 'SYNCH_5', 'TITLE',
       'ORIGIN', 'DATE', 'TIME'], dtype=object)

Plot out $\beta$ and dispersion functions with MADX elements

First create the MADX elements graphics

In [7]:
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

In [8]:
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);

Creating transport matrices from MADX results

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:

In [9]:
dX = 1e-4; dPX = 0.5e-5; dY = 1e-4; dPY = 1e-5; dPT = 0.5e-4; # from the booster.madx
In [10]:
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

In [11]:
df_dPT.head()
Out[11]:
NAME S L X PX Y PY
0 BOOSTER$START 0.000 0.000 0.0 0.0 0.0 0.0
1 START_1 0.000 0.000 0.0 0.0 0.0 0.0
2 SA 0.088 0.176 0.0 0.0 0.0 0.0
3 HS24 0.188 0.024 0.0 0.0 0.0 0.0
4 VS24 0.212 0.024 0.0 0.0 0.0 0.0

Calculate the matrix elements

In [12]:
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

In [13]:
df_twiss.tail()
Out[13]:
NAME S L TILT ANGLE X PX Y PY PT ... R33 R43 R14 R24 R34 R44 R16 R26 R36 R46
959 MINS 470.464106 0.500000 0.0 0.000000 0.0 0.0 0.0 0.0 0.0 ... -0.010044 0.130596 0.17954 0.03008 -7.75907 1.39163 3.320480 0.428825 0.026777 -0.006990
960 FMAGD24 472.158912 2.889612 0.0 0.070742 0.0 0.0 0.0 0.0 0.0 ... 0.216152 0.139876 0.21929 0.01393 -5.77779 0.88435 3.869969 0.174617 0.016143 -0.005407
961 SC 473.903718 0.600000 0.0 0.000000 0.0 0.0 0.0 0.0 0.0 ... 0.483942 0.164842 0.22543 -0.00377 -4.66420 0.47624 3.892914 -0.099026 0.007778 -0.004431
962 MARK24 474.203718 0.000000 0.0 0.000000 0.0 0.0 0.0 0.0 0.0 ... 0.533395 0.164842 0.22430 -0.00377 -4.52133 0.47624 3.863208 -0.099026 0.006449 -0.004431
963 BOOSTER$END 474.203718 0.000000 0.0 0.000000 0.0 0.0 0.0 0.0 0.0 ... 0.533395 0.164842 0.22430 -0.00377 -4.52133 0.47624 3.863208 -0.099026 0.006449 -0.004431

5 rows × 45 columns

In [14]:
df_twiss.tail(1)
Out[14]:
NAME S L TILT ANGLE X PX Y PY PT ... R33 R43 R14 R24 R34 R44 R16 R26 R36 R46
963 BOOSTER$END 474.203718 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.533395 0.164842 0.2243 -0.00377 -4.52133 0.47624 3.863208 -0.099026 0.006449 -0.004431

1 rows × 45 columns

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.

In [15]:
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)) )
1-turn transport matrix M =
[[ -2.03832000e-01  -3.30282000e+01  -7.49000000e-04   2.24300000e-01]
 [  2.99220000e-02  -5.47000000e-02  -2.65000000e-03  -3.77000000e-03]
 [  4.54000000e-04  -4.71100000e-01   5.33395000e-01  -4.52133000e+00]
 [  1.25700000e-03   1.85200000e-02   1.64842000e-01   4.76240000e-01]]
Coupled Q1    = 0.7292009
twiss.outx Q1 = 6.7292030
Coupled Q2    = 0.8344151
twiss.outx Q2 = 6.8344117

Calculating the orbit response from the transport matrices

First create the functions that finds the orbit response of coasting beam without RPOS feedback

In [16]:
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

Dispersion function

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}. $$
In [17]:
[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
In [18]:
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

Check

In [19]:
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)

RPOS

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.

df_twiss[df_twiss["NAME"]=="RPOS"]

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

In [20]:
rfactor = 1.0;
In [21]:
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

Check without RPOS for HL11

In [22]:
dipole = "HL02";
xresp, yresp = coasting_xy_response_to(dipole);
df_twiss['dx_d' + dipole + 'coast'] = xresp;
df_twiss['dy_d' + dipole + 'coast'] = yresp;
In [23]:
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)

Check with RPOS for HL11

In [24]:
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;
In [25]:
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)

Measured ORM

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

In [26]:
dir_m = "/Users/cytan/expt/booster/13Mar2017/"
In [27]:
df_ORM = pd.read_csv(dir_m+"13Mar2017_ORM_nb.txt", delim_whitespace=True )
In [28]:
df_ORM.head()
Out[28]:
Dipole Dipole_base_name s HS24slope HS24intercept HS24std_err HL01slope HL01intercept HL01std_err HS01slope ... VS22std_err VL23slope VL23intercept VL23std_err VS23slope VS23intercept VS23std_err VL24slope VL24intercept VL24std_err
47 HS24 S24 0.260 -1.90040 -0.07407 0.237230 0.10009 5.0342 0.055319 1.33840 ... 0.032003 0.095871 -1.6124 0.029904 -0.004359 0.18922 0.024111 0.160460 0.42975 0.058518
0 HL01 L01 11.629 0.35163 -0.38526 0.069987 -0.30850 4.9474 0.067034 -0.15692 ... 0.058416 0.033833 -1.6763 0.052726 0.056027 0.13534 0.048521 0.009270 0.41792 0.046888
24 HS01 S01 20.018 1.46350 -0.80726 0.074073 -0.57039 4.5570 0.045247 -2.30870 ... 0.018754 0.000057 -1.7111 0.038081 0.065246 0.12505 0.014674 0.121380 0.47432 0.038765
1 HL02 L02 32.334 0.74475 -0.59164 0.078350 0.32993 4.7472 0.057168 0.29747 ... 0.028056 -0.055304 -1.6893 0.041744 -0.040239 0.11852 0.027224 -0.030561 0.50915 0.049816
25 HS02 S02 39.777 1.42580 -0.34374 0.096052 1.25530 5.0003 0.068613 2.07970 ... 0.018218 -0.013665 -1.7111 0.035313 -0.063617 0.15333 0.012205 0.011926 0.40189 0.066220

5 rows × 309 columns

Now we can calculate the model orbit response

In [29]:
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.')
Working on VL24...
Done.
In [30]:
df_twiss.head()
Out[30]:
NAME S L TILT ANGLE X PX Y PY PT ... dx_dVL23 dy_dVL23 dx_dVS23coast dy_dVS23coast dx_dVS23 dy_dVS23 dx_dVL24coast dy_dVL24coast dx_dVL24 dy_dVL24
0 BOOSTER$START 0.000 0.000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.065001 10.345412 0.245388 3.461599 0.338741 3.461780 0.192193 -2.740644 0.663973 -2.739731
1 START_1 0.000 0.000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.065001 10.345412 0.245388 3.461599 0.338741 3.461780 0.192193 -2.740644 0.663973 -2.739731
2 SA 0.088 0.176 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.064059 10.371126 0.245211 3.534199 0.338560 3.534383 0.192696 -2.576797 0.664455 -2.575868
3 HS24 0.188 0.024 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.062989 10.400346 0.245010 3.616700 0.338354 3.616887 0.193267 -2.390607 0.665003 -2.389660
4 VS24 0.212 0.024 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.062732 10.407359 0.244962 3.636500 0.338305 3.636688 0.193404 -2.345922 0.665134 -2.344970

5 rows × 431 columns

Corrector kick vs current

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}}. $$

In [31]:
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
MAD-X Twiss file GAMMA = 1.44 => proton kinetic energy = 415 MeV

Plot Measured ORM

Setup for the plot

In [32]:
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')]
In [33]:
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

Measured ORM to horizontal kicks

In [34]:
lay = plot_responses(xDipoles, 'BPM response to x-kick (mm/mrad)')
show(lay)

Measured ORM to vertical kicks

In [35]:
lay = plot_responses(yDipoles, 'BPM response to y-kick (mm/mrad)')
show(lay)
In [ ]: