In [1]:
cd ~/Workspace/gm2/studies/field-shimming/ring-geometry
In [2]:
%matplotlib inline
from IPython.display import display, Markdown
import numpy as np
from scipy.optimize import curve_fit
from scipy.interpolate import interp1d
import matplotlib as mpl
import matplotlib.pyplot as plt
import ROOT as r
import seaborn as sns
# Set the figure size to span the notebook.
mpl.rcParams['figure.figsize'] = (18, 6)
mpl.style.use('fivethirtyeight')
mpl.rcParams['axes.grid'] = False
colors = sns.hls_palette(10, l=0.3, s=0.8)
def finish_plot(ylabel=r'z [$\mu$ m]'):
plt.xlabel(r'$\theta$ [deg]')
plt.xlim(0, 360)
plt.ylabel(ylabel)
yoke_label = ord('A')
for i in xrange(0, 37):
x = (i * 10.0 - 15.0 - 0.1012) % 360.0
if i % 3 == 0:
plt.axvline(x, linestyle='-', color='k', alpha=0.2)
else:
plt.axvline(x, linestyle='--', color='k', alpha=0.2)
if i is not 0:
if i < 10:
s = str(i) + ' '
else:
s = str(i)
plt.figtext(0.052 + i * (0.868 / 36.0), 0.86, s, color='k', alpha=0.4)
if (i % 3) == 1:
plt.figtext(0.052 + i * (0.868 / 36.0), 0.82, chr(yoke_label), color='k', alpha=0.4)
yoke_label += 1
plt.show()
def tilt_plane(phi, amp=832.173, baseline=1.059, phase=171.41):
return amp * (np.cos((phi - phase) * np.pi / 180.0) + baseline)
def sinusoid(phi, amp, baseline, phase):
return amp * np.cos((phi - phase) * np.pi / 180.0) + baseline
def pole2(x, c1, c2, a0, a1, a2):
v = (a0 + c1 - a1 * 5.0 + a2 * 5.0**2) * (x < -5.0)
v += (a0 + c2 + a1 * 5.0 + a2 * 5.0**2) * (x > 5.0)
v += (a0 + a1 * x + a2 * x**2) * (x >= -5.0) * (x <= 5.0)
return v
def pole3(x, c1, c2, a0, a1, a2, a3):
v = (a0 + c1 - a1 * 5.0 + a2 * 5.0**2 - a3 * 5.0**3) * (x < -5.0)
v += (a0 + c2 + a1 * 5.0 + a2 * 5.0**2 + a3 * 5.0**3) * (x > 5.0)
v += (a0 + a1 * x + a2 * x**2 + a3 * x**3) * (x >= -5.0) * (x <= 5.0)
return v
def pole4(x, c1, c2, a0, a1, a2, a3, a4):
return a0 + a1 * x + a2 * x**2 + a3 * x**3 + a4 * x**4
def pole5(x, c1, c2, a0, a1, a2, a3, a4, a5):
v = (a0 + c1 - a1 * 5.0 + a2 * 5.0**2 - a3 * 5.0**3 + a4 * 5.0**4 - a5 * 5.0**5) * (x < -5.0)
v += (a0 + c2 + a1 * 5.0 + a2 * 5.0**2 + a3 * 5.0**3 + a4 * 5.0**4 + a5 * 5.0**5) * (x > 5.0)
v += (a0 + a1 * x + a2 * x**2 + a3 * x**3 + a4 * x**4 + a5 * x**5) * (x >= -5.0) * (x <= 5.0)
return v
Global Bottom Pole Adjustment Plan¶
With the ring leveling complete, we are about to start a program of adjusting all the poles. The plan needs to make heavy use of our tilt measurements and laser tracker, since those devices are calibrated to gravity.
In [3]:
display(Markdown('## The Data'))
# Open the laser data from the last full scan.
f = r.TFile('data-in/full_scan_025.root')
t = f.Get('t')
z_laser = np.empty(t.GetEntries())
phi_laser = np.empty(t.GetEntries())
for i in xrange(t.GetEntries()):
t.GetEntry(i)
z_laser[i] = t.z_2 * 1e6 # Convert inches to microns
phi_laser[i] = (t.phi_2 + 0.368) % 360
# Remove points with bad data.
indices = np.where(z_laser > 0.5)
z_laser = z_laser[indices]
phi_laser = phi_laser[indices]
# Subtract the height offset to put things in the same space as our reconstruction.
z_laser -= z_laser.mean()
z_laser -= z_laser.min() + 110
# Construct weights based on azimuthal spacing
w_laser = np.empty(phi_laser.shape)
t_phi = (phi_laser - 0.368) % 360
w_laser[1:-1] = 0.5 * (t_phi[2:] - 2 * t_phi[1:-1] + t_phi[:-2])
w_laser[0] = 0.5 * (t_phi[1] - 2 * t_phi[0] + t_phi[-1] - 360.0)
w_laser[-1] = 0.5 * (t_phi[0] + 360.0 - 2 * t_phi[-1] + t_phi[-2])
w_laser = np.abs(w_laser)
plt.plot(w_laser)
plt.show()
display(Markdown('### Laser Height vs Azimuth'))
plt.scatter(phi_laser, z_laser, color=colors[0], alpha=0.5)
finish_plot()
# Open the data for pole boundaries (elog 432) and azimuthal tilts (elog 430)
d_pb = np.genfromtxt('data-in/tilt_pole_boundaries_02_22_2016.csv', skip_header=True)
d_tilt_cal = np.genfromtxt('data-in/rad_tilt_cal_03_13_2016.csv', skip_header=True)
d_tilt = np.genfromtxt('data-in/pole_tilts_03_13_2016.csv', skip_header=True)
d_tilt_azi = (d_tilt[:, 2] - d_tilt_cal[:, 2].mean()) * 1.29
# Take the average of the two steps.
dz = (2 * d_pb[:, 4] - d_pb[:, 2] - d_pb[:, 6]) * 0.5
dz *= 1.29 * 2.25 * 0.0254 # Convert the value to microns
drad = (d_pb[:, 5] - d_pb[:, 3]) * 1.29 * 11.75 * 0.0254
dz += 0.5 * drad
dz_outer = dz - 0.5 * drad
# plt.plot(np.linspace(5, 355, 36), dz)
# plt.plot(np.linspace(5, 355, 36), dz - 0.5 * drad)
# plt.plot(np.linspace(5, 355, 36), dz + 0.5 * drad)
# finish_plot()
dzds = np.array(d_tilt_azi)
# Refill the z values and shift into [0, 360]
z = np.empty(dz.shape[0] * 100)
phi = np.linspace(0, 360, z.shape[0])
for i in xrange(z.shape[0]):
# To set between [0, 360] instead of [-5, 355]
idx = (i - 50) % z.shape[0]
z[idx] = dz[:i/100].sum()
z[idx] += dzds[i/100] * 0.012 * (i % 100)
z[idx] += dzds[:i/100].sum() * 1.2
# Divide the different by the number of weight of pieces in between z_i and z_f
dzds0 = (z[-51] - z[-50] - dz[-1]) / (1.2 * 35.99)
# Add the offset to the slopes
dzds -= dzds0
# Refill the z values and shift into [0, 360]
z = np.empty(dz.shape[0] * 100)
phi = np.linspace(0, 360, z.shape[0])
for i in xrange(z.shape[0]):
# To set between [0, 360] instead of [-5, 355]
idx = (i - 50) % z.shape[0]
z[idx] = dz[:i/100].sum()
z[idx] += dzds[i/100] * 0.012 * (i % 100)
z[idx] += dzds[:i/100].sum() * 1.2
display(Markdown('### Tilt Model (With Closure)'))
plt.scatter(phi_laser, z_laser, color=colors[2], alpha=0.5)
plt.plot(phi, z, c=colors[0])
finish_plot()
Dealing with the ring tilt plane and laser tilt plane¶
One can observe from differencing the simple tilt model and laser data, that there are residual tilt planes still involved in this story. One tilt plane is the ring with respect to gravity which we can get a proxy for by comparing the tilt model with a flat plane. The second tilt plane is the laser tracker with respect to the ring tilt plane (or gravity, but we are using the ring as our basis). This one we tease out by looking for one omega between the laser data and the corrected tilt sensor model. I ignore yoke D in both of these calculations, since we know it to live above the plane by our choice.
In [4]:
display(Markdown('### Ring Tilt Plane'))
# Pick out all points besides yoke D
idx1 = np.where((phi - 75.0) % 360 > 30.0)
f_z = z[idx1]
f_phi = phi[idx1]
p, cov = curve_fit(sinusoid, f_phi, f_z)
print "(amp, baseline, phase) => (%.2f, %.2f, %.2f)" % (p[0], p[1], p[2] % 360.0)
plt.plot(phi, z, color=colors[0], alpha=0.3)
plt.plot(phi, z - sinusoid(phi, p[0], p[1], p[2]), color=colors[0])
finish_plot()
z_corr = np.array(z - sinusoid(phi, p[0], p[1], p[2]))
display(Markdown('### Laser Height vs Tilt Model (Ring Tilt Plane Removed)'))
f = interp1d(phi, z_corr)
plt.scatter(phi_laser, z_laser, color=colors[0])
plt.scatter(phi_laser, f(phi_laser), color=colors[2])
finish_plot()
display(Markdown('### Laser Tilt Plane (Difference of Laser Height and Tilt Model)'))
# Put the data onto an evenly spaced grid
f_z_laser = interp1d(phi_laser, z_laser - f(phi_laser))
phi_laser_grid = np.linspace(0.5, 359.5, 5001)
z_laser_grid = f_z_laser(phi_laser_grid)
p, cov = curve_fit(sinusoid, phi_laser_grid, z_laser_grid)
print "(amp, baseline, phase) => (%.2f, %.2f, %.2f)" % (p[0], p[1], p[2] % 360.0)
plt.plot(phi_laser_grid, sinusoid(phi_laser_grid, p[0], p[1], p[2]), color=colors[2])
plt.plot(phi_laser_grid, z_laser_grid, color=colors[0])
finish_plot()
plt.scatter(phi_laser, z_laser - sinusoid(phi_laser, p[0], p[1], p[2]), color=colors[0])
plt.scatter(phi, z_corr, color=colors[2])
finish_plot()
z_laser_corr = np.array(z_laser - sinusoid(phi_laser, p[0], p[1], p[2]))
Fitting for Curvature of Poles¶
We need to estimate how to fix the curvature of each pole. About the best we can do here is use the laser data.
In [5]:
# Try fitting each pole to a parabola
fit_params = np.empty([36, 8])
fit_tilt = np.empty(36)
bounds = (np.array([-200, -200, -250.0, -50, -4, -0.4, -0.16, -0.032]), np.array([200, 200, 250.0, 50, 0, 0.4, 0.0, 0.032]))
for i in xrange(36):
f_phi = (phi_laser - i * 10 + 10) % 360
f_z = np.array(z_laser_corr)
indices = np.where(f_phi > 0.0)
f_phi = f_phi[indices]
f_z = f_z[indices]
f_w = w_laser[indices]
indices = np.where(f_phi < 20.0)
f_phi = f_phi[indices]
f_z = f_z[indices]
f_w = f_w[indices]
indices = f_phi.argsort()
f_phi = f_phi[indices]
f_z = f_z[indices]
f_phi -= 10.0
if f_phi.shape[0] == 0:
continue
fit_params[i, :] = np.array([0, 0, 0, 0, 0, 0, 0, 0])
else:
bounds[0][0] = -dz[i-1]
bounds[1][0] = -dz[i-1] + 0.001
bounds[0][1] = dz[i]
bounds[1][1] = dz[i] + 0.001
ff = interp1d(f_phi, f_z)
f_phi = np.linspace(-6.0, 6.0, 501)
f_z = ff(f_phi)
sigma = np.ones(f_phi.shape[0])
sigma[np.where(f_phi < -5.0)] *= 0.01
sigma[np.where(f_phi > 5.0)] *= 0.01
p, pcor = curve_fit(pole5, f_phi, f_z, bounds=bounds, sigma=sigma, max_nfev=10000)
plt.scatter(f_phi, f_z, color=colors[i%10], alpha=0.8)
fit_params[i, :] = p
display(Markdown('### Laser Data For All Poles'))
plt.show()
for i in xrange(36):
f_phi = np.linspace(10.0 * i - 4.99, 10.0 * i + 4.99, 100)
p_phi = np.linspace(-4.99, 4.99, 100)
par = fit_params[i, :]
plt.plot(f_phi, pole5(p_phi, par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7]), color=colors[0])
fit_tilt[i] = pole5(1.2214, par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7])
fit_tilt[i] -= pole5(-1.2214, par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7])
fit_tilt[i] /= 0.29845
display(Markdown('### Quintic Fits to Each Pole Curvature'))
plt.scatter(phi_laser, z_laser_corr, color=colors[2], alpha=0.5)
finish_plot()
Radial Pole Tilts¶
The last piece of the puzzle here is inflecting the correct radial pole tilt on the bottom pole pieces. We have data from field off measurements, and we began the process of taking field on measurements today, 3/15/2016. Comparing the two and knowing that we want to shoot for 45 $\mu rad$ of radial closure with the field on, we can make a plan for the radial adjustment in our bottom pole movement program.
In [6]:
# Load the data for round two.
d1_tilt_rad_cal = np.genfromtxt('data-in/rad_tilt_cal_1_03_07_2016.csv')
d1_tilt_rad_cal_2 = np.genfromtxt('data-in/rad_tilt_cal_2_03_07_2016.csv')
d1_tilt = np.genfromtxt('data-in/pole_tilts_03_07_2016.csv', skip_header=True)
# Set the tilts from day 4 where we measured the radial tilt well, and the azimuthal symbiotically.
d1_tilt_rad = (d1_tilt[:, 1] - d1_tilt_rad_cal[:, 2].mean()) * 1.29
d1_tilt_rad[:9] = (d1_tilt[:9, 1] - d1_tilt_rad_cal[:, 2].mean()) * 1.29
d1_tilt_rad[9:19] = (d1_tilt[9:19, 1] - d1_tilt_rad_cal_2[:, 2].mean()) * 1.29
d1_tilt_rad[19:] = (d1_tilt[19:, 1] - d1_tilt_rad_cal[:, 2].mean()) * 1.29
d2_tilt_rad_cal = np.genfromtxt('data-in/rad_tilt_cal_03_16_2016.csv', skip_header=True)
d2_tilt = np.genfromtxt('data-in/pole_tilts_03_16_2016.csv', skip_header=True)
#d2_tilt_rad = (d2_tilt[:, 2] - 31053.5) * 1.29
d2_tilt_rad = (d2_tilt[:, 2] - d2_tilt_rad_cal[:, 2].mean()) * 1.29
plt.plot(np.linspace(0, 350, 36), d1_tilt_rad, color=colors[0])
plt.plot(np.linspace(0, 350, 36), d2_tilt_rad, color=colors[2])
display(Markdown('### Measured Radial Tilts with Field On/Off'))
finish_plot('radial tilt [$\mu rad$]')
display(Markdown('### Radial Tilt Differences with Field On/Off'))
plt.hist(d2_tilt_rad - d1_tilt_rad, bins=20, color=colors[0])
plt.show()
Pole Adjustments¶
We can use the steps from tilt sensor measurements, the curvature from most of the polynomial fits, and the radial tilts from measurements field on/off. With all of this info together, we can create a pole movement plan which might still need some spot checking.
In [8]:
# Set target height based on Yoke C
z0 = 55.0
drad = 45.0 - d2_tilt_rad
display(Markdown('### Pole Adjustment Height Goal'))
plt.scatter(phi_laser, z_laser_corr, color=colors[2], alpha=0.75)
for i in xrange(-4, 10):
plt.axhline(z0 + i * 25, color='gray', linestyle='--', alpha=0.25)
plt.axhline(z0 + 0.0, color='gray', alpha=0.5)
finish_plot()
display(Markdown('### Pole Adjustment Radial Tilt Goal'))
plt.plot(np.linspace(0, 350, 36), 45.0 - d2_tilt_rad + d1_tilt_rad, color='gray', alpha=0.8)
plt.plot(np.linspace(0, 350, 36), d1_tilt_rad, color=colors[0])
finish_plot('radial tilt [$\mu rad$]')
# Calculate pole foot adjustments
# Convenience variables
moves = np.zeros([36, 10])
mean_pole_z = np.zeros(36)
output = ['| Radius | A | B | C | D | E |']
output.append('|-------|------|------|------|------|------|')
output.append('| Inner | %.2f | %.2f | %.2f | %.2f | %.2f |')
output.append('| Outer | %.2f | %.2f | %.2f | %.2f | %.2f |')
for i in xrange(36):
# Calculate the height displacement
idx = np.where((phi_laser - (10.0 * i - 5.0)) % 360 < 10.0)
t_z = z_laser_corr[idx]
t_phi = phi_laser[idx]
t_dz = z0 - t_z.mean()
moves[i, :] = t_dz
if i not in [19, 20]:
par = fit_params[i, :]
for j in xrange(5):
a = -4.788 + 2.394 * j
moves[i, j] -= pole5(a, par[0], par[1], 0.0, par[3], par[4], par[5], par[6], par[7])
moves[i, j + 5] -= pole5(a, par[0], par[1], 0.0, par[3], par[4], par[5], par[6], par[7])
else:
for j in xrange(5):
f = interp1d(phi, z)
t_dz0 = (f(t_phi) - t_z).mean()
a = 10.0 * i - 4.788 + 2.394 * j
moves[i, j] -= f(a) - t_dz0
moves[i, j + 5] -= f(a) - t_dz0
# Take care of the radial tilt as an average
for j in xrange(5):
moves[i, j + 5] += (d2_tilt_rad[i] - 45.0) * 0.51
# Now take care bad pole steps
t_moves = np.array(moves)
for i in xrange(36):
t_dz_inner = -t_moves[i, 0] + t_moves[i - 1, 4]
t_dz_outer = -t_moves[i, 5] + t_moves[i - 1, 9]
a3_inner = 0.5 * (t_dz_inner - dz[i-1]) / 4.878**3
a3_outer = 0.5 * (t_dz_outer - dz_outer[i-1]) / 4.878**3
for j in xrange(0, 2):
a = -4.788 + 2.394 * j
moves[i, j] += a3_inner * np.abs(a**3)
moves[i, j + 5] += a3_outer * np.abs(a**3)
t_dz_inner = -t_moves[(i + 1) % 36, 0] + t_moves[i, 4]
t_dz_outer = -t_moves[(i + 1) % 36, 5] + t_moves[i, 9]
a3_inner = 0.5 * (t_dz_inner - dz[i]) / 4.878**3
a3_outer = 0.5 * (t_dz_outer - dz_outer[i]) / 4.878**3
for j in xrange(3, 5):
a = -4.788 + 2.394 * j
moves[i, j] -= a3_inner * np.abs(a**3)
moves[i, j + 5] -= a3_outer * np.abs(a**3)
# print "% .2f % .2f % .2f % .2f % .2f" % tuple(moves[i, :5] - t_moves[i, :5])
# print "% .2f % .2f % .2f % .2f % .2f" % tuple(moves[i, 5:] - t_moves[i, 5:])
display(Markdown('### Pole Adjustment Plan'))
plt.scatter(phi_laser, z_laser_corr, color='#bdbdbd', alpha=1.0)
for i in xrange(36):
t_phi = np.linspace(10.0 * i - 4.788, 10.0 * i + 4.788, 5)
plt.scatter(t_phi, z0 - moves[i, :5], alpha=0.9, s=50.0, color=colors[i%10])
plt.scatter(t_phi, z0 - moves[i, 5:], alpha=0.9, s=50.0, marker='^', color=colors[i%10])
for i in xrange(-4, 10):
plt.axhline(z0 + i * 25, color='gray', linestyle='--', alpha=0.25)
plt.axhline(z0 + 0.0, color='gray', alpha=0.5)
finish_plot()
display(Markdown('### Perfect Pole Adjustment Results'))
final_z = np.array(z_laser_corr)
for i in xrange(phi_laser.shape[0]):
n = int(((phi_laser[i] + 5.0) % 360) / 10.0)
par = fit_params[n, :]
a = (phi_laser[i] + 5.0) % 10.0 - 5.0
final_z[i] -= pole5(a, par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7]) - z0
plt.scatter(phi_laser, final_z, color=colors[2], alpha=0.8)
plt.axhline(z0 + 0.0, color='gray', alpha=0.5)
finish_plot()
display(Markdown('### Pole Steps'))
plt.plot(np.linspace(5, 355, 36), dz, color=colors[0], alpha=0.2)
plt.plot(np.linspace(5, 355, 36), dz_outer, color=colors[2], alpha=0.2)
plan_dz_inner = np.empty(36)
plan_dz_outer = np.empty(36)
plan_dz_inner[:35] = moves[1:, 0] - moves[:-1, 4]
plan_dz_inner[-1] = moves[0, 0] - moves[-1, 4]
plan_dz_outer[:35] = moves[1:, 5] - moves[:-1, 9]
plan_dz_outer[-1] = moves[0, 5] - moves[-1, 9]
plt.plot(np.linspace(5, 355, 36), -plan_dz_inner, color=colors[0], alpha=0.8)
plt.plot(np.linspace(5, 355, 36), -plan_dz_outer, color=colors[2], alpha=0.8)
plt.axhline(50, color='gray', alpha=0.5)
plt.axhline(25, color='gray', alpha=0.5)
plt.axhline(0, color='gray', alpha=0.5)
plt.axhline(-25, color='gray', alpha=0.5)
plt.axhline(-50, color='gray', alpha=0.5)
finish_plot()
display(Markdown('### Inner Feet Shim Changes'))
plt.hist(moves[:, :5].flatten() / 25.4, bins=40, color=colors[0])
plt.show()
display(Markdown('### Outer Feet Shim Changes'))
plt.hist(moves[:, 5:].flatten() / 25.4, bins=40, color=colors[0])
plt.show()
# Print the values.
for i in xrange(36):
display(Markdown('### Pole %i Movements' % (i + 1)))
display(Markdown('\n'.join(output) % tuple(moves[i, :] / 25.4)))
print d2_tilt_rad.mean()
print (moves[:, :5].flatten().mean() - moves[:, 5:].flatten().mean()) / 0.51
In [ ]: