add: Faraday cylinder normalize;fix: bind formula;change: cal and pro class
This commit is contained in:
parent
8cbcc431f1
commit
98b652fec7
16
main.py
16
main.py
|
@ -1,21 +1,19 @@
|
|||
import numpy as np
|
||||
from qdx import Process
|
||||
|
||||
p1 = Process("coef1.csv", n=5)
|
||||
p2 = Process("coef2.csv", n=5)
|
||||
p1 = Process()
|
||||
p2 = Process()
|
||||
|
||||
# Mg25 + d -> Mg26* + p
|
||||
angles = [i for i in np.arange(8, 64, 4)] + [10]
|
||||
for angle in angles:
|
||||
task = "Mg25(dp)Mg26/angle-{:d}.csv".format(angle)
|
||||
task = "result/Mg25(dp)Mg26/angle-{:d}.csv".format(angle)
|
||||
print("Processing Task :", task)
|
||||
|
||||
p1.read_data(task)
|
||||
p1("coef1.csv", task, n=5)
|
||||
p1.energy_filter(12, 30)
|
||||
p1.draw_filter("Mg25(dp)Mg26/figure/filter/angle-{:d}.png".format(angle))
|
||||
p1.draw_result("Mg25(dp)Mg26/figure/result/angle-{:d}.png".format(angle))
|
||||
p1.draw_result("result/Mg25(dp)Mg26/figure/angle-{:d}.png".format(angle))
|
||||
|
||||
p2.read_data(task)
|
||||
p2("coef2.csv", task, n=5)
|
||||
p2.energy_filter(12, 30)
|
||||
p2.draw_filter("Mg25(dp)Mg26/figure/filter/angle-{:d}-o.png")
|
||||
p2.draw_result("Mg25(dp)Mg26/figure/result/angle-{:d}-o.png".format(angle))
|
||||
p2.draw_result("result/Mg25(dp)Mg26/figure/angle-{:d}-o.png".format(angle))
|
||||
|
|
12
qdx/Bind.py
12
qdx/Bind.py
|
@ -10,10 +10,10 @@ 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: $x = -\\frac{k_2}{k_1}y + \\frac{E-b_1-b_2}{k_1}$.
|
||||
Plotting the data in a 2D plane for linear regression, we can get $\\frac{k_2}{k_1}$ and $C_i$
|
||||
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_1 = \\frac{E_1 - E_2}{C_1 - C_2}$.
|
||||
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$.
|
||||
|
@ -30,8 +30,8 @@ class Bind(object):
|
|||
$b = b_1 + b_2$
|
||||
|
||||
K1/K2/C1/C2: float
|
||||
$K_i = -\\frac{k_2}{k_1}$
|
||||
$C_i = \\frac{E_i-b}{k_1}$
|
||||
$K_i = -\\frac{k_1}{k_2}$
|
||||
$C_i = \\frac{E_i-b}{k_2}$
|
||||
"""
|
||||
|
||||
def __init__(self, n, m):
|
||||
|
@ -82,7 +82,7 @@ class Bind(object):
|
|||
self.px[1] = data[:, 2]
|
||||
|
||||
def get_line(self):
|
||||
"""Fit data with $x = -\\frac{k_2}{k_1}y + \\frac{E-b_1-b_2}{k_1}$."""
|
||||
"""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])
|
||||
|
|
|
@ -9,74 +9,69 @@ from .utils import readBlockData, get_hist
|
|||
|
||||
|
||||
class Calibration(object):
|
||||
"""Calibrate the detector according to the calibration data
|
||||
"""Calibrate the detector according to the calibration data"""
|
||||
|
||||
Parameters
|
||||
----------
|
||||
bias_b : float
|
||||
bias $b = b_1 + b_2$
|
||||
delta_e : float
|
||||
delta energy between two beams
|
||||
n : int, optional
|
||||
number of blocks, default 6
|
||||
m : int, optional
|
||||
number of binds, default 8
|
||||
"""
|
||||
def __init__(self):
|
||||
pass
|
||||
|
||||
def __init__(self, bias_b, delta_e, n=6, m=8):
|
||||
self.bias = bias_b
|
||||
self.delta_e = delta_e
|
||||
self.n = n
|
||||
self.m = m
|
||||
|
||||
self.binds = [[Bind(i, j) for j in range(m)] for i in range(n)]
|
||||
|
||||
def __call__(self, file1, file2):
|
||||
def __call__(self, file1, file2, bias_b, delta_e, n=6, m=8):
|
||||
"""Calibration
|
||||
|
||||
Parameters
|
||||
----------
|
||||
file1/2 : str
|
||||
data file path of energy 1/2
|
||||
bias_b : float
|
||||
bias $b = b_1 + b_2$, in MeV
|
||||
delta_e : float
|
||||
delta energy between two beams, in MeV
|
||||
n : int, optional
|
||||
number of blocks, default 6
|
||||
m : int, optional
|
||||
number of binds, default 8
|
||||
"""
|
||||
# Initialization
|
||||
self.n, self.m = n, m
|
||||
self.binds = [[Bind(i, j) for j in range(m)] for i in range(n)]
|
||||
|
||||
# Read Data
|
||||
file_list = csv.reader(open(file1, "r"))
|
||||
pbar = tqdm(desc="Read Data E1", total=len(open(file1, "r").readlines()))
|
||||
for row in file_list:
|
||||
pn = int(row[1])
|
||||
ldata, rdata = readBlockData(row[0], pn, self.m)
|
||||
for i in range(self.m):
|
||||
ldata, rdata = readBlockData(row[0], int(row[5]), pn, m)
|
||||
for i in range(m):
|
||||
bind = self.binds[pn][i]
|
||||
bind.add_data(0, ldata[i], rdata[i], float(row[2]))
|
||||
bind.add_data(0, ldata[i], rdata[i], float(row[4]))
|
||||
pbar.update(1)
|
||||
pbar.close()
|
||||
|
||||
file_list = csv.reader(open(file2, "r"))
|
||||
pbar = tqdm(desc="Read Data E1", total=len(open(file2, "r").readlines()))
|
||||
pbar = tqdm(desc="Read Data E2", total=len(open(file2, "r").readlines()))
|
||||
for row in file_list:
|
||||
pn = int(row[1])
|
||||
ldata, rdata = readBlockData(row[0], pn, self.m)
|
||||
for i in range(self.m):
|
||||
ldata, rdata = readBlockData(row[0], int(row[5]), pn, m)
|
||||
for i in range(m):
|
||||
bind = self.binds[pn][i]
|
||||
bind.add_data(1, ldata[i], rdata[i], float(row[2]))
|
||||
bind.add_data(1, ldata[i], rdata[i], float(row[4]))
|
||||
pbar.update(1)
|
||||
pbar.close()
|
||||
|
||||
# Data preprocessing
|
||||
pbar = tqdm(desc="Bind Process", total=self.n * self.m)
|
||||
for i in range(self.n):
|
||||
for j in range(self.m):
|
||||
pbar = tqdm(desc="Bind Process", total=n * m)
|
||||
for i in range(n):
|
||||
for j in range(m):
|
||||
bind: Bind = self.binds[i][j]
|
||||
bind.slash()
|
||||
bind.get_line()
|
||||
bind.get_kb(self.bias, self.delta_e)
|
||||
bind.get_kb(bias_b, delta_e)
|
||||
pbar.update(1)
|
||||
pbar.close()
|
||||
|
||||
# Fit
|
||||
pbar = tqdm(desc="Bind Fit", total=self.n * self.m)
|
||||
for i in range(self.n):
|
||||
for j in range(self.m):
|
||||
pbar = tqdm(desc="Bind Fit", total=n * m)
|
||||
for i in range(n):
|
||||
for j in range(m):
|
||||
bind: Bind = self.binds[i][j]
|
||||
bind.get_peak_center()
|
||||
bind.fit_px()
|
||||
|
@ -89,8 +84,9 @@ class Calibration(object):
|
|||
Parameters
|
||||
----------
|
||||
path : str
|
||||
save folder, there must be `Fit-line`, `GMM`, and `PEAK` 3 subfolders.
|
||||
save folder, there must be `FIT-LINE`, `GMM`, and `PEAK` 3 subfolders.
|
||||
"""
|
||||
pbar = tqdm(desc="Draw Fit Figure", total=self.n * self.m)
|
||||
for i in range(self.n):
|
||||
for j in range(self.m):
|
||||
bind: Bind = self.binds[i][j]
|
||||
|
@ -98,6 +94,9 @@ class Calibration(object):
|
|||
bind.draw_cluster(os.path.join(path, "GMM", bind.name + ".png"))
|
||||
bind.draw_peak(os.path.join(path, "PEAK", bind.name + ".png"))
|
||||
|
||||
pbar.update(1)
|
||||
pbar.close()
|
||||
|
||||
def draw_check(self, path="Check.png"):
|
||||
"""Draw check figure
|
||||
|
||||
|
@ -106,13 +105,15 @@ class Calibration(object):
|
|||
path : str, optional
|
||||
save path
|
||||
"""
|
||||
pbar = tqdm(desc="Draw Figure Check", total=self.n * self.m)
|
||||
fig = plt.figure(figsize=(24, 15), dpi=200)
|
||||
ax1 = fig.add_subplot(2, 1, 1)
|
||||
ax2 = fig.add_subplot(2, 1, 2)
|
||||
pbar = tqdm(desc="Draw Check Figure", total=self.n * self.m)
|
||||
fig = plt.figure(figsize=(24, 24), dpi=200)
|
||||
ax1 = fig.add_subplot(3, 1, 1)
|
||||
ax2 = fig.add_subplot(3, 1, 2)
|
||||
ax3 = fig.add_subplot(3, 1, 3)
|
||||
peaks = np.array([])
|
||||
|
||||
for i in range(self.n):
|
||||
engs = np.array([])
|
||||
for j in range(self.m):
|
||||
bind = self.binds[i][j]
|
||||
peaks = np.hstack((np.unique(bind.px[0]), peaks))
|
||||
|
@ -120,23 +121,34 @@ class Calibration(object):
|
|||
eng = bind.predict_energy(bind.x[0], bind.y[0])
|
||||
pX = bind.predict_px(bind.x[0], bind.y[0])
|
||||
count, center = get_hist(pX, delta=0.5)
|
||||
ax1.scatter(pX, eng, s=0.1, color="k")
|
||||
ax2.scatter(center, count + 3000 * (7 - j), s=0.5, color="k")
|
||||
engs = np.hstack((engs, eng))
|
||||
|
||||
ax2.scatter(pX, eng, s=0.1, color="k")
|
||||
ax3.scatter(center, count + 3000 * (7 - j), s=0.5, color="k")
|
||||
|
||||
pbar.update(1)
|
||||
count, center = get_hist(engs, delta=0.01)
|
||||
ax1.plot(center, count)
|
||||
|
||||
peaks = np.unique(peaks)
|
||||
for x in peaks:
|
||||
ax2.vlines(x, 0, 3000 * self.m, color="gray", linestyles="dashed")
|
||||
ax3.vlines(x, 0, 3000 * self.m, color="gray", linestyles="dashed")
|
||||
for j in range(self.m):
|
||||
ax2.hlines(
|
||||
2500 * j,
|
||||
ax3.hlines(
|
||||
3000 * j,
|
||||
(np.min(peaks) // 50) * 50,
|
||||
(np.min(peaks) // 50 + 1) * 50,
|
||||
(np.max(peaks) // 50 + 1) * 50,
|
||||
color="r",
|
||||
linestyles="dashdot",
|
||||
)
|
||||
|
||||
ax1.set_xlabel("Energy (MeV)")
|
||||
ax1.set_ylabel("Count per bin")
|
||||
ax2.set_xlabel("x (mm)")
|
||||
ax2.set_ylabel("Energy (MeV)")
|
||||
ax3.set_xlabel("x (mm)")
|
||||
ax3.set_ylabel("Count per bin")
|
||||
|
||||
fig.savefig(path, facecolor="w", transparent=False)
|
||||
plt.close()
|
||||
pbar.close()
|
||||
|
|
|
@ -10,46 +10,40 @@ from .utils import readFileData, get_hist
|
|||
|
||||
|
||||
class Process(object):
|
||||
"""Process the experimental data according to the calibration results.
|
||||
"""Process the experimental data according to the calibration results."""
|
||||
|
||||
Parameters
|
||||
----------
|
||||
path : str
|
||||
def __init__(self) -> None:
|
||||
pass
|
||||
|
||||
def __call__(self, coef, task, n=6, m=8):
|
||||
"""Read Process Data
|
||||
coef : str
|
||||
coefficient file
|
||||
task : str
|
||||
task file
|
||||
n : int, optional
|
||||
number of blocks, default 6
|
||||
m : int, optional
|
||||
number of binds, default 8
|
||||
"""
|
||||
|
||||
def __init__(self, path, n=6, m=8) -> None:
|
||||
self.n = n
|
||||
self.m = m
|
||||
|
||||
# Initialization
|
||||
self.n, self.m = n, m
|
||||
self.binds = [[Bind(i, j) for j in range(m)] for i in range(n)]
|
||||
|
||||
# Read Calibration Data
|
||||
pbar = tqdm(desc="Bind Initialization", total=n * m)
|
||||
data = list(csv.reader(open(path, "r")))
|
||||
data = list(csv.reader(open(coef, "r")))
|
||||
data = np.array(data, dtype=np.float64)
|
||||
for i in range(n):
|
||||
for j in range(m):
|
||||
bind = self.binds[i][j]
|
||||
bind(data[j + i * m][2:])
|
||||
|
||||
pbar.update(1)
|
||||
pbar.close()
|
||||
|
||||
def read_data(self, file):
|
||||
"""Read Process Data
|
||||
|
||||
file : str
|
||||
task file path
|
||||
"""
|
||||
# Read Data
|
||||
n, m = self.n, self.m
|
||||
total = len(open(file, "r").readlines()) * self.n * self.m
|
||||
file_list = csv.reader(open(file, "r"))
|
||||
total = len(open(task, "r").readlines()) * n * m
|
||||
file_list = csv.reader(open(task, "r"))
|
||||
|
||||
self.pX = np.array([])
|
||||
self.eng = np.array([])
|
||||
|
@ -104,31 +98,6 @@ class Process(object):
|
|||
self.eng_n = y
|
||||
self.reg = reg
|
||||
|
||||
def draw_filter(self, path="filter.png"):
|
||||
"""Draw the energy filter result
|
||||
|
||||
Parameters
|
||||
----------
|
||||
path : str, optional
|
||||
save path
|
||||
"""
|
||||
px_min = (np.min(self.pX_n) // 50) * 50
|
||||
px_max = (np.max(self.pX_n) // 50 + 1) * 50
|
||||
px_x = np.linspace(px_min, px_max, int(px_max - px_min))
|
||||
|
||||
fig = plt.figure(figsize=(20, 8), dpi=200)
|
||||
ax = fig.add_subplot(1, 1, 1)
|
||||
|
||||
ax.scatter(self.pX, self.eng, s=0.01, color="black")
|
||||
ax.scatter(self.pX_n, self.eng_n, s=0.01, color="orange")
|
||||
ax.plot(px_x, self.reg(px_x))
|
||||
|
||||
ax.set_xlabel("x (mm)")
|
||||
ax.set_ylabel("Energy (MeV)")
|
||||
|
||||
fig.savefig(path)
|
||||
plt.close()
|
||||
|
||||
def draw_result(self, path="result.png"):
|
||||
"""Draw the processing result
|
||||
|
||||
|
@ -137,23 +106,26 @@ class Process(object):
|
|||
path : str, optional
|
||||
save path
|
||||
"""
|
||||
fig = plt.figure(figsize=(20, 8), dpi=200)
|
||||
fig = plt.figure(figsize=(24, 12), dpi=200)
|
||||
ax1 = fig.add_subplot(2, 1, 1)
|
||||
ax2 = fig.add_subplot(2, 1, 2)
|
||||
|
||||
count, center = get_hist(self.pX_n, delta=0.1)
|
||||
ax1.scatter(self.pX_n, self.eng_n, s=0.05, color="k")
|
||||
ax1.scatter(self.pX, self.eng, s=0.01, color="black")
|
||||
ax1.scatter(self.pX_n, self.eng_n, s=0.01, color="orange")
|
||||
ax2.step(center, count, where="post", color="k")
|
||||
|
||||
px_min = (np.min(self.pX_n) // 50) * 50
|
||||
px_max = (np.max(self.pX_n) // 50 + 1) * 50
|
||||
px_x = np.linspace(px_min, px_max, int(px_max - px_min))
|
||||
ax1.plot(px_x, self.reg(px_x))
|
||||
|
||||
ax1.set_xticks(np.arange(px_min, px_max, 50))
|
||||
ax2.set_xticks(np.arange(px_min, px_max, 50))
|
||||
|
||||
ax1.set_xlabel("x (mm)")
|
||||
ax1.set_ylabel("Energy (MeV)")
|
||||
ax2.set_xlabel("x (mm)")
|
||||
ax2.set_ylabel("Count per bin")
|
||||
|
||||
fig.savefig(path)
|
||||
fig.savefig(path, facecolor="w", transparent=False)
|
||||
plt.close()
|
||||
|
|
16
qdx/utils.py
16
qdx/utils.py
|
@ -4,13 +4,15 @@ import matplotlib.pyplot as plt
|
|||
from sklearn.mixture import GaussianMixture
|
||||
|
||||
|
||||
def readFileData(file, n=6, m=8, minT=800, maxT=4000):
|
||||
def readFileData(file, count, n=6, m=8, minT=800, maxT=4000):
|
||||
"""Read whole data from root file
|
||||
|
||||
Parameters
|
||||
----------
|
||||
file : str
|
||||
root file path
|
||||
count : int
|
||||
count that normalized by counts of Faraday cylinder
|
||||
n : int, optional
|
||||
number of blocks, default 6
|
||||
m : int, optional
|
||||
|
@ -26,8 +28,8 @@ def readFileData(file, n=6, m=8, minT=800, maxT=4000):
|
|||
for j in range(m):
|
||||
na = i // 2
|
||||
nc = j + 2 * m * (i % 2)
|
||||
x = data["adc{:d}ch{:d}".format(na, nc)].array(library="np")
|
||||
y = data["adc{:d}ch{:d}".format(na, nc + m)].array(library="np")
|
||||
x = data["adc{:d}ch{:d}".format(na, nc)].array(library="np")[:count]
|
||||
y = data["adc{:d}ch{:d}".format(na, nc + m)].array(library="np")[:count]
|
||||
idx = np.where((x + y >= minT) & (x + y <= maxT))[0]
|
||||
ldata.append(x[idx])
|
||||
rdata.append(y[idx])
|
||||
|
@ -35,13 +37,15 @@ def readFileData(file, n=6, m=8, minT=800, maxT=4000):
|
|||
return ldata, rdata
|
||||
|
||||
|
||||
def readBlockData(file, n, m=8, minT=800, maxT=4000):
|
||||
def readBlockData(file, count, n, m=8, minT=800, maxT=4000):
|
||||
"""Read block data from root file
|
||||
|
||||
Parameters
|
||||
----------
|
||||
file : str
|
||||
root file path
|
||||
count : int
|
||||
count that normalized by counts of Faraday cylinder
|
||||
n : int
|
||||
No.n block
|
||||
m : int, optional
|
||||
|
@ -56,8 +60,8 @@ def readBlockData(file, n, m=8, minT=800, maxT=4000):
|
|||
for j in range(m):
|
||||
na = n // 2
|
||||
nc = j + 2 * m * (n % 2)
|
||||
x = data["adc{:d}ch{:d}".format(na, nc)].array(library="np")
|
||||
y = data["adc{:d}ch{:d}".format(na, nc + m)].array(library="np")
|
||||
x = data["adc{:d}ch{:d}".format(na, nc)].array(library="np")[:count]
|
||||
y = data["adc{:d}ch{:d}".format(na, nc + m)].array(library="np")[:count]
|
||||
idx = np.where((x + y >= minT) & (x + y <= maxT))[0]
|
||||
ldata.append(x[idx])
|
||||
rdata.append(y[idx])
|
||||
|
|
Loading…
Reference in New Issue