In [3]:
from __future__ import print_function, division
from datetime import datetime, timedelta
from pytimber import LoggingDB
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from itertools import cycle
ldb = LoggingDB()
In [4]:
# BSRT
sigma = ['LHC.BSRT.5L4.B2:FIT_SIGMA_V', 'LHC.BSRT.5L4.B2:FIT_SIGMA_H']
lsf = ['LHC.BSRT.5L4.B2:LSF_V', 'LHC.BSRT.5L4.B2:LSF_H']
pos = ['LHC.BSRT.5L4.B2:FIT_POSITION_V', 'LHC.BSRT.5L4.B2:FIT_POSITION_H']
beta = ['LHC.BSRT.5L4.B2:BETA_V', 'LHC.BSRT.5L4.B2:BETA_H'] # This needs timestamps from the hole fill
# BPM
bpm_data = ['LHC.BOFSU:POSITIONS_V', 'LHC.BOFSU:POSITIONS_H']
bpm_nvar = ['LHC.BOFSU:BPM_NAMES_V', 'LHC.BOFSU:BPM_NAMES_H']
bpm_n = ['.5L4.B2', '.6L4.B2', '.7L4.B2']
bpm_d = [135.338, 172.338, 268.942]
# CTRL - BGV
fillnum = 6371
bgv_beta = {'h': 225.76, 'v': 158.93} # \pm .43, \pm .13
bgv_pos = 220 # make sure position signs etc are OK w.r.t BPMs
In [5]:
fill = ldb.getLHCFillData(fillnum, unixtime=False)
tf1, tf2 = fill['startTime'], fill['endTime']
print('Fill %s\n\tStart: %s\tEnd :%s' % (fill["fillNumber"], tf1, tf2))
for item in fill['beamModes']:
print("Mode: %s\n\tStart: %s\tEnd :%s" % (item["mode"], item['startTime'], item['endTime']))
In [6]:
# Set the times that I want
mode = "STABLE"
for item in fill['beamModes']:
if item["mode"] == mode:
# t1, t2 = item['startTime']+timedelta(minutes=10), item['endTime']
t1, t2 = item['startTime'], item['endTime']
# t1, t2 = datetime(2017, 11, 2, 19, 5, 0), datetime(2017, 11, 2, 19, 20, 0)
In [7]:
bsrt = {}
In [9]:
# For the beta I need the fill timestamps
# one measurement at the beginning, one at the end (hopefully)
d_beta = ldb.get(beta, tf1, tf2, unixtime=False)
minFmt = mdates.DateFormatter('%Y-%m-%d %H:%M')
f = plt.figure()
ax = f.add_subplot(111)
for k, d in d_beta.items():
print("%s\t: x: %s\ty: %s" % (k, str(d[0].size), str(d[1].shape)))
x = d[0]
y = d[1]
# get the last beta measurement before the beggining of the daq
try:
ind_min = np.argmin(np.absolute(filter(lambda x: x.days < 0, np.subtract(x, t1))))
bsrt["b_%s" % k[-1].lower()] = round(y[ind_min], 1)
print('\tIndex = %d\tValue = %.1f' % (ind_min, y[ind_min]))
except ValueError:
# most probably lack of measurements
bsrt["b_%s" % k[-1].lower()] = round(y[0], 1)
ax.plot(x, y, lw=3, label=k)
ax.set_ylabel('Value [m]', fontsize=12)
ax.legend(loc=3, numpoints=1, framealpha=.7)
ax.grid(which='both', lw=2)
ax.set_ylim((100, 400))
ax.set_title('[BSRT] Beta')
ax.xaxis.set_major_formatter(minFmt)
f.autofmt_xdate()
plt.show()
print(bsrt)
In [17]:
d_lsf = ldb.get(lsf, tf1, tf2, unixtime=False)
for k, d in d_lsf.iteritems():
print("%s\t: x: %s\ty: %s" % (k, str(d[0].size), str(d[1].shape)))
x = d[0]
y = d[1]
# get the last beta measurement before the beggining of the daq
try:
ind_min = np.argmin(np.absolute(filter(lambda x: x.days < 0, np.subtract(x, t1))))
bsrt["c_%s" % k[-1].lower()] = y[ind_min]
print('\tIndex = %d\tValue = %.3f' % (ind_min, y[ind_min]))
except ValueError:
# most probably lack of measurements
bsrt["c_%s" % k[-1].lower()] = round(y[0], 1)
print(bsrt)
In [18]:
norm = 6900/0.938 # normalization factor for emmitance
window = 60
def fill_array_with_nan(data):
lens = np.array([len(i) for i in data])
# Mask of valid places in each row
mask = np.arange(lens.max()) < lens[:,None]
# Setup output array and put elements from data into masked positions
ret = np.empty(mask.shape, dtype=np.float)
ret[:] = np.nan
ret[mask] = np.concatenate(data)
return set(lens), ret
def mavg(a, n=30) :
ret = np.cumsum(a, dtype=float)
ret[n:] = ret[n:] - ret[:-n]
return ret[n - 1:] / n
d_sigma = ldb.get(sigma, t1, t2, unixtime=False)
minFmt = mdates.DateFormatter('%Y-%m-%d %H:%M:%S')
ma = cycle(('+','1', 'x', '2', '3', '4', '|', '_'))
f = plt.figure(figsize=(8, 12))
ax = f.add_subplot(311)
ax2 = f.add_subplot(312)
ax3 = f.add_subplot(313)
for k, d in d_sigma.iteritems():
print("%s\t: x: %s\ty: %s" % (k, str(d[0].size), str(d[1].shape)))
s = k[-1].lower()
x = d[0]
y = d[1]
if len(y.shape) == 1:
lens, y = fill_array_with_nan(d[1])
print("\tLens: %s\ty: %s" % (lens, str(y.shape)))
y = np.nanmean(y, axis=1)
mu = y.mean()
print('\tMean = %.3f\tSigma = %.3f' % (mu, y.std()))
bsrt["s_m%s" % s] = mu
bsrt["sc_m%s" % s] = np.sqrt(mu**2 - bsrt["c_%s" % s]**2)
y2 = np.sqrt(y**2 - bsrt["c_%s" % s]**2)
y3 = norm*(y2**2/bsrt['b_h'])
ax.plot(x, y, ls='None', marker=next(ma), label=k)
ax.plot(x[window-1:], mavg(y, n=window), ls='--', c='k')
ax2.plot(x, y2, ls='None', marker=next(ma), label=k)
ax2.plot(x[window-1:], mavg(y2, n=window), ls='--', c='k')
ax3.plot(x, y3, ls='None', marker=next(ma), label=k)
ax3.plot(x[window-1:], mavg(y3, n=window), ls='--', c='k')
for a in [ax, ax2, ax3]:
a.set_ylabel('Value [mm2]', fontsize=12)
a.legend(loc=1, numpoints=1, framealpha=.7)
a.grid(which='both', lw=2)
ax.set_title('[BSRT] Sigma')
ax2.set_title('[BSRT] Sigma Corrected')
ax3.set_title('Normalized emmitance')
ax3.xaxis.set_major_formatter(minFmt)
f.autofmt_xdate()
plt.show()
In [19]:
# e = norm*(sigma**2 / beta(x))
norm = 6900/0.938 # normalization factor for emmitance
print('\n*** Mean Beam Size @BSRT\n\tmX=%.3f\tmY=%.3f\nmXY=%.3f'
% (bsrt['s_mh'], bsrt['s_mv'], np.mean([bsrt['s_mh'], bsrt['s_mv']])))
e1x = norm*(bsrt['s_mh']**2/bsrt['b_h'])
e1y = norm*(bsrt['s_mv']**2/bsrt['b_v'])
print('\teX=%.3f\teY=%.3f' % (e1x, e1y))
print('\n*** Mean Beam Size @BSRT - Corrected (LSF)\n\tmX=%.3f\tmY=%.3f\nmXY=%.3f'
% (bsrt['sc_mh'], bsrt['sc_mv'], np.mean([bsrt['sc_mh'], bsrt['sc_mv']])))
e2x = norm*(bsrt['sc_mh']**2/bsrt['b_h'])
e2y = norm*(bsrt['sc_mv']**2/bsrt['b_v'])
print('\teX=%.3f\teY=%.3f' % (e2x, e2y))
# Extrapolating @BGV
print('\n*** Extrapolating @BGV')
mu1x = bsrt['s_mh']*np.sqrt(bgv_beta['h'] / bsrt['b_h'])
mu1y = bsrt['s_mv']*np.sqrt(bgv_beta['v'] / bsrt['b_v'])
print("\tmX = %.3f\tmY = %.3f" % (mu1x, mu1y))
print('mXY = %.3f' % (np.mean([mu1x, mu1y])))
em1x = norm*mu1x**2/bgv_beta['h']
em1y = norm*mu1y**2/bgv_beta['v']
print('\teX=%.3f\teY=%.3f' % (em1x, em1y))
print('\n*** Extrapolating @BGV - Corrected (LSF)')
mu2x = bsrt['sc_mh']*np.sqrt(bgv_beta['h'] / bsrt['b_h'])
mu2y = bsrt['sc_mv']*np.sqrt(bgv_beta['v'] / bsrt['b_v'])
print("\tmX = %.3f\tmY = %.3f" % (mu2x, mu2y))
print('mXY = %.3f' % (np.mean([mu2x, mu2y])))
print('mXY = %.3f' % (np.sqrt((mu2x**2+mu2y**2)/2.0)))
em2x = norm*mu2x**2/bgv_beta['h']
em2y = norm*mu2y**2/bgv_beta['v']
print('\teX=%.3f\teY=%.3f' % (em2x, em2y))
In [10]:
d_pos = ldb.get(pos, t1, t2)
minFmt = mdates.DateFormatter('%Y-%m-%d %H:%M:%S')
ma = cycle(('+','1', 'x', '2', '3', '4'))
f = plt.figure()
ax = f.add_subplot(111)
for k, d in d_pos.iteritems():
print("%s\t: x: %s\ty: %s" % (k, str(d[0].size), str(d[1].shape)))
x = [datetime.utcfromtimestamp(ts) for ts in d[0]]
y = d[1]
if len(y.shape) == 1:
lens, y = fill_array_with_nan(d[1])
print("\tLens: %s\ty: %s" % (lens, str(y.shape)))
y = np.nanmean(y, axis=1)
print('\tMean = %.3f\tSigma = %.3f' % (y.mean(), y.std()))
bsrt["s_m%s" % k[-1].lower()] = y.mean()
bsrt["s_s%s" % k[-1].lower()] = y.std()
ax.plot(x, y, ls='None', marker=next(ma), label=k)
ax.set_ylabel('Value [mm]', fontsize=12)
ax.legend(loc=5, numpoints=1, framealpha=.7)
ax.grid(which='both', lw=2)
ax.set_title('[BSRT] Position')
ax.xaxis.set_major_formatter(minFmt)
f.autofmt_xdate()
plt.show()
In [11]:
bpm = {}
In [12]:
# Retrieve data from LogDB
n_bpm = ldb.get(bpm_nvar, tf1, tf2)
d_bpm = ldb.get(bpm_data, t1, t2)
for k, d in n_bpm.iteritems():
print("%s\t: x: %s\ty: %s" % (k, str(d[0].size), str(d[1].shape)))
for k, d in d_bpm.iteritems():
print("%s\t: x: %s\ty: %s" % (k, str(d[0].size), str(d[1].shape)))
In [ ]:
indexes = {'ih': [], 'iv': []}
tmp = []
for k, d in n_bpm.iteritems():
print('%s' % k)
for ind, name in enumerate(d[1][0]):
for t in bpm_n:
if t in name:
indexes['i%s' % k[-1].lower()].append((int(ind), name))
print('\t%d\t%s - %s' % (ind, k, name))
print('\n')
f = plt.figure()
for cnt, (k, d) in enumerate(d_bpm.iteritems()):
ma = cycle(('+','1', 'x', '2', '3', '4'))
ax = f.add_subplot(int('21%s' % (cnt+1)))
print("%s\t: x: %s\ty: %s" % (k, str(d[0].size), str(d[1].shape)))
x = [datetime.utcfromtimestamp(ts) for ts in d[0]]
s = k[-1].lower()
for t in indexes['i%s' % s]:
i, n = t
y = d[1][:,i]
mu, si = y.mean(), y.std()
print('\t%s\t: Mean = %.3f\tSigma = %.3f' % (n, mu, si))
try:
bpm[n][s] = (mu, si)
except KeyError:
bpm[n] = {}
bpm[n][s] = (mu, si)
ax.plot(x, y, ls='None', marker=next(ma), label=n)
ax.set_ylabel('Value [um]', fontsize=12)
ax.legend(loc=1, numpoints=1, framealpha=.7, bbox_to_anchor=(1.44, 1))
ax.grid(which='both', lw=2)
ax.set_title('[BPM] %s' % k)
ax.xaxis.set_major_formatter(minFmt)
f.autofmt_xdate()
plt.show()
In [ ]:
# Extrapolating @BGV
# BGV is after 5L4, 6L4 @-220m
for k, d in zip(bpm_n, bpm_d):
for key in bpm.keys():
if k in key:
print("%s @ %s\t(%.1f, %.1f)" % (k, -d, bpm[key]['h'][0], bpm[key]['v'][0]))
In [ ]:
p1, p2, p3 = 'BPMYB.5L4.B2', 'BPMYA.6L4.B2', 'BPM.7L4.B2'
print('Extrapolating 5L4-6L4')
slope_h = (bpm[p2]['h'][0] - bpm[p1]['h'][0]) / (bpm_d[1] - bpm_d[0])
slope_v = (bpm[p2]['v'][0] - bpm[p1]['v'][0]) / (bpm_d[1] - bpm_d[0])
pos_h = np.rint(bpm[p2]['h'][0] + slope_h*(bgv_pos-bpm_d[1]))
pos_v = np.rint(bpm[p2]['v'][0] + slope_v*(bgv_pos-bpm_d[1]))
print('\t@BGV: (%d, %d) um' % (pos_h, pos_v))
print('Interpolating 6L4-7L4')
slope_h = (bpm[p3]['h'][0] - bpm[p2]['h'][0]) / (bpm_d[2] - bpm_d[1])
slope_v = (bpm[p3]['v'][0] - bpm[p2]['v'][0]) / (bpm_d[2] - bpm_d[1])
pos_h = np.rint(bpm[p2]['h'][0] + slope_h*(bgv_pos-bpm_d[1]))
pos_v = np.rint(bpm[p2]['v'][0] + slope_v*(bgv_pos-bpm_d[1]))
print('\t@BGV: (%d, %d) um' % (pos_h, pos_v))
In [ ]: