297 lines
9.3 KiB
Python
297 lines
9.3 KiB
Python
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
|
|
from .utils import get_hist, GMM_clip
|
|
from .model import Linear1D, FixedSlopeLine, pXLine
|
|
from .fit import fit_line, fit_hist_gaussian
|
|
|
|
|
|
class Bind(object):
|
|
"""Bind of a Block, using data under two energy conditions for fitting.
|
|
The energy of an ADC device can be expressed as: $E = kx + b$.
|
|
So the complete energy is represented by $E(x,y) = k_1x + b_1 + k_2y + b_2$.
|
|
This formula can be converted to: $y = -\\frac{k_1}{k_2}x + \\frac{E-b_1-b_2}{k_2}$.
|
|
Plotting the data in a 2D plane for linear regression, we can get $\\frac{k_1}{k_2}$ and $C_i$
|
|
|
|
Using data of $E_1$ and $E_2$, now $k_2 = \\frac{E_1 - E_2}{C_1 - C_2}$.
|
|
|
|
The incident position produces different responses in the left and right detectors, which can be expressed as:
|
|
$Px = \\frac{k_2y - k_1x}{E}L+CL$.
|
|
|
|
Attributes
|
|
----------
|
|
L/C : float
|
|
parameters in $Px = \\frac{k_2y - k_1x}{E}L+CL$.
|
|
|
|
k1, k2: float
|
|
parameters in $Px = \\frac{k_2y - k_1x}{E}L+CL$.
|
|
|
|
b : float
|
|
$b = b_1 + b_2$
|
|
|
|
K1/K2/C1/C2: float
|
|
$K_i = -\\frac{k_1}{k_2}$
|
|
$C_i = \\frac{E_i-b}{k_2}$
|
|
"""
|
|
|
|
def __init__(self, n, m):
|
|
self.n, self.m = n, m
|
|
self.px = [None, None]
|
|
self.x, self.y = [None, None], [None, None]
|
|
|
|
self.L, self.C = 40, 0
|
|
self.k1, self.k2, self.b = 0, 0, 0
|
|
self.K1, self.C1, self.K2, self.C2 = 0, 0, 0, 0
|
|
|
|
def add_data(self, k, x, y, px):
|
|
"""k-th energy data
|
|
|
|
Parameters
|
|
----------
|
|
k : int
|
|
k-th energy
|
|
x : array
|
|
left data
|
|
y : array
|
|
right data
|
|
px : int
|
|
px (set value)
|
|
"""
|
|
self.x[k] = np.concatenate((self.x[k], x)) if self.x[k] is not None else x
|
|
self.y[k] = np.concatenate((self.y[k], y)) if self.y[k] is not None else y
|
|
self.px[k] = (
|
|
np.concatenate((self.px[k], np.full(len(x), px)))
|
|
if self.px[k] is not None
|
|
else np.full(len(x), px)
|
|
)
|
|
|
|
def clip(self):
|
|
"""Using Gaussian Mixture Method (GMM) to decompose the data into noise and available data"""
|
|
data = GMM_clip(
|
|
np.array(list(zip(self.x[0], self.y[0], self.px[0])), dtype=object)
|
|
)
|
|
self.x[0] = data[:, 0]
|
|
self.y[0] = data[:, 1]
|
|
self.px[0] = data[:, 2]
|
|
|
|
data = GMM_clip(
|
|
np.array(list(zip(self.x[1], self.y[1], self.px[1])), dtype=object)
|
|
)
|
|
self.x[1] = data[:, 0]
|
|
self.y[1] = data[:, 1]
|
|
self.px[1] = data[:, 2]
|
|
|
|
def get_line(self):
|
|
"""Fit data with $y = -\\frac{k_1}{k_2}x + \\frac{E-b_1-b_2}{k_2}$."""
|
|
model = Linear1D()
|
|
|
|
reg = fit_line(model, self.x[0], self.y[0])
|
|
self.K1 = reg.slope.value
|
|
self.C1 = reg.intercept.value
|
|
|
|
reg = fit_line(model, self.x[1], self.y[1])
|
|
self.K2 = reg.slope.value
|
|
self.C2 = reg.intercept.value
|
|
|
|
def get_kb(self, bias, deltaE=4):
|
|
"""Get $k_2$, $k_2$, $b$ from $K_i$, $C_i$
|
|
Set slope equal average of $K_i$.
|
|
|
|
Parameters
|
|
----------
|
|
bias : float
|
|
bias $b = b_1 + b_2$
|
|
deltaE : float, optional
|
|
delta energy between two beams
|
|
"""
|
|
K = (self.K1 + self.K2) / 2
|
|
self.K1 = self.K2 = K
|
|
|
|
model = FixedSlopeLine(slope=K, intercept=self.C1)
|
|
fitted_model = fit_line(model, self.x[0], self.y[0])
|
|
self.C1 = fitted_model.intercept.value
|
|
|
|
model = FixedSlopeLine(slope=K, intercept=self.C2)
|
|
fitted_model = fit_line(model, self.x[1], self.y[1])
|
|
self.C2 = fitted_model.intercept.value
|
|
|
|
self.k2 = deltaE / abs(self.C1 - self.C2)
|
|
self.k1 = -self.k2 * K
|
|
|
|
eng = self.k1 * self.x[0] + self.k2 * self.y[0]
|
|
fitted_model = fit_hist_gaussian(eng, step=0.01)
|
|
self.E1 = fitted_model.mean.value
|
|
|
|
eng = self.k1 * self.x[1] + self.k2 * self.y[1]
|
|
fitted_model = fit_hist_gaussian(eng, step=0.01)
|
|
self.E2 = fitted_model.mean.value
|
|
|
|
self.b = bias - self.C1 * self.k2
|
|
|
|
self.E1 += self.b
|
|
self.E2 += self.b
|
|
|
|
def get_peak_center(self):
|
|
"""Get peak center (in channel) using Gaussian Model"""
|
|
self.fx, self.fy, self.fz = [], [], []
|
|
|
|
peaks = np.unique(self.px[0])
|
|
for px in peaks:
|
|
idx = np.where(self.px[0] == px)[0]
|
|
x, y = self.x[0][idx], self.y[0][idx]
|
|
|
|
if len(idx) < 400:
|
|
continue
|
|
|
|
fitted_model = fit_hist_gaussian(x)
|
|
self.fx.append(fitted_model.mean.value)
|
|
|
|
fitted_model = fit_hist_gaussian(y)
|
|
self.fy.append(fitted_model.mean.value)
|
|
|
|
self.fz.append(px)
|
|
|
|
self.fx = np.array(self.fx)
|
|
self.fy = np.array(self.fy)
|
|
self.fz = np.array(self.fz)
|
|
|
|
def fit_px(self):
|
|
"""Fit using $Px = \\frac{k_2y - k_1x}{E}L+CL$."""
|
|
model = Linear1D()
|
|
reg = fit_line(model, self.fx, self.fz)
|
|
|
|
E = np.mean(self.k1 * self.fx + self.k2 * self.fy + self.b)
|
|
L = reg.slope.value * E / (self.k2 * self.K1 - self.k1)
|
|
C = (reg.intercept.value - self.k2 * self.C1 * L / E) / L
|
|
|
|
model = pXLine(L=L, C=C, k1=self.k1, k2=self.k2, b=self.b)
|
|
reg = fit_line(
|
|
model, np.array(list(zip(self.fx, self.fy)), dtype=object), self.fz
|
|
)
|
|
self.L = reg.L.value
|
|
self.C = reg.C.value
|
|
|
|
def predict_energy(self, x, y):
|
|
"""Use $E(x,y) = k_1x + b_1 + k_2y + b_2$. to calculate energy.
|
|
|
|
Parameters
|
|
----------
|
|
x/y : array
|
|
data
|
|
"""
|
|
return self.k1 * x + self.k2 * y + self.b
|
|
|
|
def predict_px(self, x, y):
|
|
"""Use $Px = \\frac{k_2y - k_1x}{E}L+CL$ to calculate pX.
|
|
|
|
Parameters
|
|
----------
|
|
x/y : array
|
|
data
|
|
"""
|
|
eng = self.predict_energy(x, y)
|
|
return (self.k2 * y - self.k1 * x) / eng * self.L + self.C * self.L
|
|
|
|
def draw_fit_line(self, title):
|
|
fig = plt.figure(figsize=(8, 8))
|
|
ax = fig.add_subplot(1, 1, 1)
|
|
ax.scatter(self.x[0], self.y[0], s=0.1, c="black", label=r"$E_1$")
|
|
ax.scatter(self.x[1], self.y[1], s=0.1, c="dimgray", label=r"$E_2$")
|
|
ax.plot(
|
|
self.x[0],
|
|
self.K1 * self.x[0] + self.C1,
|
|
c="red",
|
|
label=r"$x_2={:.4f}x_1+{:.4f},\ R^2={:.5f}$".format(
|
|
self.K1, self.C1, self.RSquare1
|
|
),
|
|
)
|
|
ax.plot(
|
|
self.x[1],
|
|
self.K2 * self.x[1] + self.C2,
|
|
c="orangered",
|
|
label=r"$x_2={:.4f}x_1+{:.4f},\ R^2={:.5f}$".format(
|
|
self.K2, self.C2, self.RSquare2
|
|
),
|
|
)
|
|
|
|
plt.title("Y - X Line")
|
|
plt.xlabel("Output - Left (X) / (channel)")
|
|
plt.ylabel("Output - Right (Y) / (channel)")
|
|
|
|
ax.legend()
|
|
fig.savefig(title, facecolor="w", transparent=False, bbox_inches="tight")
|
|
plt.close()
|
|
|
|
def draw_peak(self, title):
|
|
fig = plt.figure(figsize=(8, 8))
|
|
|
|
peaks = np.unique(self.px[0])
|
|
n, k = len(peaks), 1
|
|
for px in peaks:
|
|
ax = fig.add_subplot(n, 1, k)
|
|
ax.set_title("pX = {:.2f}".format(px))
|
|
|
|
idx = np.where(self.px[0] == px)[0]
|
|
x, y = self.x[0][idx], self.y[0][idx]
|
|
|
|
count1, center1 = get_hist(x)
|
|
count2, center2 = get_hist(y)
|
|
|
|
ax.scatter(center1, count1, s=0.5)
|
|
ax.scatter(center2, count2, s=0.5)
|
|
|
|
ax.set_xlabel("Channel")
|
|
ax.set_ylabel("Counts per Channel")
|
|
|
|
k += 1
|
|
|
|
fig.suptitle("Counts Curve")
|
|
plt.tight_layout()
|
|
fig.savefig(title, facecolor="w", transparent=False)
|
|
plt.close()
|
|
|
|
def draw_cluster(self, title):
|
|
fig = plt.figure(figsize=(8, 8))
|
|
ax = fig.add_subplot(1, 1, 1)
|
|
|
|
peaks = np.unique(self.px[0])
|
|
for px in peaks:
|
|
idx = np.where(self.px[0] == px)[0]
|
|
x, y = self.x[0][idx], self.y[0][idx]
|
|
ax.scatter(x, y, s=0.1, label="pX={:.0f}".format(px))
|
|
|
|
plt.title("Y - X Line")
|
|
plt.xlabel("Output - Left (X) / (channel)")
|
|
plt.ylabel("Output - Right (Y) / (channel)")
|
|
|
|
plt.legend()
|
|
plt.tight_layout()
|
|
fig.savefig(title, facecolor="w", transparent=False)
|
|
plt.close()
|
|
|
|
@property
|
|
def name(self):
|
|
return "{:d}-{:d}-{:.2f}-{:.2f}".format(self.n, self.m, self.E1, self.E2)
|
|
|
|
def _r_square(self, y, yp):
|
|
mean = np.mean(y)
|
|
SST = np.sum((y - mean) ** 2)
|
|
SSE = np.sum((y - yp) ** 2)
|
|
return 1 - SSE / SST
|
|
|
|
@property
|
|
def RSquare1(self):
|
|
return self._r_square(self.y[0], self.K1 * self.x[0] + self.C1)
|
|
|
|
@property
|
|
def RSquare2(self):
|
|
return self._r_square(self.y[1], self.K2 * self.x[1] + self.C2)
|
|
|
|
def __call__(self, data):
|
|
"""Data is read to complete initialization"""
|
|
self.k1 = data[0]
|
|
self.k2 = data[1]
|
|
self.b = data[2]
|
|
self.L = data[3]
|
|
self.C = data[4]
|