From 98b652fec7109626bdba3a896874a70910e8d323 Mon Sep 17 00:00:00 2001 From: YiHui Liu Date: Wed, 27 Jul 2022 11:33:47 +0800 Subject: [PATCH] add: Faraday cylinder normalize;fix: bind formula;change: cal and pro class --- main.py | 16 +++---- qdx/Bind.py | 12 +++--- qdx/calibration.py | 104 +++++++++++++++++++++++++-------------------- qdx/process.py | 80 ++++++++++++---------------------- qdx/utils.py | 16 ++++--- 5 files changed, 107 insertions(+), 121 deletions(-) diff --git a/main.py b/main.py index 051a131..6453847 100644 --- a/main.py +++ b/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)) diff --git a/qdx/Bind.py b/qdx/Bind.py index c7b070a..82d0832 100644 --- a/qdx/Bind.py +++ b/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]) diff --git a/qdx/calibration.py b/qdx/calibration.py index e9977a8..52a034d 100644 --- a/qdx/calibration.py +++ b/qdx/calibration.py @@ -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() diff --git a/qdx/process.py b/qdx/process.py index 7ebc0d7..ccd2af0 100644 --- a/qdx/process.py +++ b/qdx/process.py @@ -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 - coefficient 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 + 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 + """ + # 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() diff --git a/qdx/utils.py b/qdx/utils.py index b6e929f..e0c2703 100644 --- a/qdx/utils.py +++ b/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])