diff --git a/calibration.py b/calibration.py index 00d8f77..db81d70 100644 --- a/calibration.py +++ b/calibration.py @@ -1,12 +1,14 @@ import csv import numpy as np from qdx import Bind -from qdx.utils import readData, get_hist +from qdx.utils import readBlockData, get_hist from tqdm import tqdm from matplotlib import pyplot as plt # Initialization n, m = 5, 8 +bias = 12.97 +deltaE = 4 binds = [[Bind(i, j) for j in range(m)] for i in range(n)] # Read Data @@ -14,7 +16,7 @@ file_list = csv.reader(open("./config1.csv", "r")) pbar = tqdm(desc="Read Data E1", total=len(open("./config1.csv", "r").readlines())) for row in file_list: pn = int(row[1]) - ldata, rdata = readData("2016Q3D/root/raw/201609Q3D" + row[0] + ".root", pn, m) + ldata, rdata = readBlockData(row[0], pn, m) for i in range(m): bind = binds[pn][i] bind.add_data(0, ldata[i], rdata[i], 585 - float(row[2])) @@ -25,7 +27,7 @@ file_list = csv.reader(open("./config2.csv", "r")) pbar = tqdm(desc="Read Data E2", total=len(open("./config2.csv", "r").readlines())) for row in file_list: pn = int(row[1]) - ldata, rdata = readData("2016Q3D/root/raw/201609Q3D" + row[0] + ".root", pn, m) + ldata, rdata = readBlockData(row[0], pn, m) for i in range(m): bind = binds[pn][i] bind.add_data(1, ldata[i], rdata[i], 585 - float(row[2])) @@ -38,12 +40,11 @@ for i in range(n): for j in range(m): bind: Bind = binds[i][j] bind.slash() - bind.get_energy() bind.get_line() - bind.get_kb() - # bind.draw_fit_line("result/FIT-LINE/" + bind.name + ".png") - # bind.draw_cluster("result/GMM/" + bind.name + ".png") - # bind.draw_peak("result/PEAK/" + bind.name + ".png") + bind.get_kb(bias, deltaE) + bind.draw_fit_line("result/FIT-LINE/" + bind.name + ".png") + bind.draw_cluster("result/GMM/" + bind.name + ".png") + bind.draw_peak("result/PEAK/" + bind.name + ".png") pbar.update(1) pbar.close() @@ -59,7 +60,7 @@ pbar.close() # Draw check figure pbar = tqdm(desc="Figure Check", total=n * m) -fig = plt.figure(figsize=(16, 16), dpi=200) +fig = plt.figure(figsize=(16, 10), dpi=200) ax1 = fig.add_subplot(2, 1, 1) ax2 = fig.add_subplot(2, 1, 2) peaks = np.array([]) @@ -72,7 +73,6 @@ for i in range(n): 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 + 2500 * (7 - j), s=0.5, color="k") @@ -82,19 +82,19 @@ peaks = np.unique(peaks) for x in peaks: ax2.vlines(x, 0, 20000, color="gray", linestyles="dashed") for j in range(m): - ax2.hlines(2500 * j, 0, 650, color="r", linestyles="dashdot") + ax2.hlines(2500 * j, -50, 600, color="r", linestyles="dashdot") fig.savefig("./result/Check.png", facecolor="w", transparent=False) plt.close() pbar.close() # Save coefficient to file -f = open("./coef.txt", "w") +f = open("./coef.csv", "w") for i in range(n): for j in range(m): bind = binds[i][j] f.writelines( - "{:d} {:d} {:.9f} {:.9f} {:.9f} {:.9f} {:.9f}\n".format( + "{:d},{:d},{:.9f},{:.9f},{:.9f},{:.9f},{:.9f}\n".format( i, j, bind.k1, bind.k2, bind.b, bind.L, bind.C ) ) diff --git a/process.py b/process.py index e69de29..eb2fd24 100644 --- a/process.py +++ b/process.py @@ -0,0 +1,69 @@ +import csv +import numpy as np +from qdx import Bind +from qdx.utils import readFileData, get_hist +from tqdm import tqdm +from matplotlib import pyplot as plt + +# Initialization +n, m = 5, 8 +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("coef1.csv", "r"))) +data = np.array(data, dtype=np.float64) +for i in range(n): + for j in range(m): + bind = binds[i][j] + bind(data[j + i * m][2:]) + + pbar.update(1) +pbar.close() + +# Read Data +total = len(open("task3.csv", "r").readlines()) * n * m +file_list = csv.reader(open("task3.csv", "r")) + +pX = np.array([]) +eng = np.array([]) + +pX_full = np.array([]) +eng_full = np.array([]) + +pbar = tqdm(desc="Task - Mg25(d,p)Mg26*", total=total) +for row in file_list: + ldata, rdata = readFileData(row[0], n, m) + for i in range(n): + for j in range(m): + bind = binds[i][j] + + x = bind.predict_px(ldata[j + i * m], rdata[j + i * m]) + float(row[1]) + e = bind.predict_energy(ldata[j + i * m], rdata[j + i * m]) + edge_l = 5 + 130 * i + float(row[1]) - 35 + edge_r = edge_l + 65 + idx = np.where((x >= edge_l) & (x <= edge_r))[0] + + pX = np.hstack((pX, x[idx])) + eng = np.hstack((eng, e[idx])) + + pX_full = np.hstack((pX_full, x)) + eng_full = np.hstack((eng_full, e)) + + pbar.update(1) +pbar.close() + +# Draw +fig = plt.figure(figsize=(20, 8), dpi=200) +ax1 = fig.add_subplot(2, 1, 1) +ax2 = fig.add_subplot(2, 1, 2) + +count, center = get_hist(pX, delta=0.1) +ax1.scatter(pX, eng, s=0.05, color="k") +py, = ax2.step(center, count, where="post", color="k") + +ax1.set_xticks(np.arange((np.min(pX) // 50) * 50, (np.max(pX) // 50 + 1) * 50, 50)) +ax2.set_xticks(np.arange((np.min(pX) // 50) * 50, (np.max(pX) // 50 + 1) * 50, 50)) + +fig.savefig("Task3.png") +plt.close() diff --git a/qdx/Bind.py b/qdx/Bind.py index 045426e..c7b070a 100644 --- a/qdx/Bind.py +++ b/qdx/Bind.py @@ -81,16 +81,56 @@ class Bind(object): self.y[1] = data[:, 1] self.px[1] = data[:, 2] - def get_energy(self): - """Get energy (in channel) by fit Gaussian Model with ldata + rdata""" - eng = self.x[0] + self.y[0] - fitted_model = fit_hist_gaussian(eng) + def get_line(self): + """Fit data with $x = -\\frac{k_2}{k_1}y + \\frac{E-b_1-b_2}{k_1}$.""" + 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, delta=0.01) self.E1 = fitted_model.mean.value - eng = self.x[1] + self.y[1] - fitted_model = fit_hist_gaussian(eng) + eng = self.k1 * self.x[1] + self.k2 * self.y[1] + fitted_model = fit_hist_gaussian(eng, delta=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 = [], [], [] @@ -115,52 +155,25 @@ class Bind(object): self.fy = np.array(self.fy) self.fz = np.array(self.fz) - def get_line(self): - """Fit data with $x = -\\frac{k_2}{k_1}y + \\frac{E-b_1-b_2}{k_1}$.""" - model = Linear1D() - - self.reg1 = fit_line(model, self.x[0], self.y[0]) - self.K1 = self.reg1.slope.value - self.C1 = self.reg1.intercept.value - - self.reg2 = fit_line(model, self.x[1], self.y[1]) - self.K2 = self.reg2.slope.value - self.C2 = self.reg2.intercept.value - - def get_kb(self): - """Get $k_2$, $k_2$, $b$ from $K_i$, $C_i$ - Set slope equal average of $K_i$. - """ - 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.k1 = (self.E1 - self.E2) / (self.C1 - self.C2) - self.k2 = -(self.K1 + self.K2) * self.k1 / 2 - self.b = (self.E1 - self.C1 * self.k1 + self.E2 - self.C2 * self.k1) / 2 - 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) - L = reg.slope.value * self.E1 / (self.k2 * self.K1 - self.k1) - C = (reg.intercept.value - self.k2 * self.C1 * L / self.E1) / L - model = pXLine(L=L, C=C, k1=self.k1, k2=self.k2, E=self.E1) - reg = fit_line(model, np.array(list(zip(self.fx, self.fy)), dtype=object), 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 @@ -176,7 +189,8 @@ class Bind(object): x/y : array data """ - return (self.k2 * y - self.k1 * x) / self.E1 * self.L + self.C * self.L + 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)) @@ -185,7 +199,7 @@ class Bind(object): ax.scatter(self.x[1], self.y[1], s=0.1, c="dimgray", label=r"$E_2$") ax.plot( self.x[0], - self.reg1(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 @@ -193,7 +207,7 @@ class Bind(object): ) ax.plot( self.x[1], - self.reg2(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 @@ -244,7 +258,7 @@ class Bind(object): @property def name(self): - return "{:d}-{:d}-{:.1f}-{:.1f}".format(self.n, self.m, self.E1, self.E2) + return "{:d}-{:d}-{:.2f}-{:.2f}".format(self.n, self.m, self.E1, self.E2) def _r_square(self, y, yp): mean = np.mean(y) @@ -254,8 +268,16 @@ class Bind(object): @property def RSquare1(self): - return self._r_square(self.y[0], self.reg1(self.x[0])) + 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.reg2(self.x[1])) + 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] diff --git a/qdx/fit.py b/qdx/fit.py index 4002bd6..186d2b3 100644 --- a/qdx/fit.py +++ b/qdx/fit.py @@ -31,7 +31,7 @@ def fit_line(model, x, y): return fitted_model -def fit_hist_gaussian(x): +def fit_hist_gaussian(x, delta=1): """ Gaussian fitting is performed on the histogram @@ -39,6 +39,8 @@ def fit_hist_gaussian(x): ---------- x : array data point + delta : int, optional + Minimum bin width. The bin width is an integer multiple of delta. Returns ------- @@ -47,7 +49,7 @@ def fit_hist_gaussian(x): """ fitter = fitting.LMLSQFitter() - count, center = get_hist(x) + count, center = get_hist(x, delta=delta) model = models.Gaussian1D(amplitude=count.max(), mean=x.mean(), stddev=x.std()) fitted_model = fitter(model, center, count) diff --git a/qdx/model.py b/qdx/model.py index fc5b860..ad467b7 100644 --- a/qdx/model.py +++ b/qdx/model.py @@ -70,15 +70,17 @@ class pXLine(Fittable2DModel): linear = True L, C = Parameter(default=40), Parameter(default=0) - k1, k2, E = Parameter(), Parameter(), Parameter() + k1, k2, b = Parameter(), Parameter(), Parameter() @staticmethod - def evaluate(x, y, L, C, k1, k2, E): + def evaluate(x, y, L, C, k1, k2, b): + E = k1 * x + k2 * y + b return (k2 * y - k1 * x) / E * L + C * L @staticmethod - def fit_deriv(x, y, L, C, k1, k2, E): + def fit_deriv(x, y, L, C, k1, k2, b): + E = k1 * x + k2 * y + b d_L = (k2 * y - k1 * x) / E + C d_C = np.full(x.shape, L) - d_k1, d_k2, d_E = np.zeros_like(x), np.zeros_like(x), np.zeros_like(x) - return [d_L, d_C, d_k1, d_k2, d_E] + d_k1, d_k2, d_b = np.zeros_like(x), np.zeros_like(x), np.zeros_like(x) + return [d_L, d_C, d_k1, d_k2, d_b] diff --git a/qdx/utils.py b/qdx/utils.py index 622e5e3..b6e929f 100644 --- a/qdx/utils.py +++ b/qdx/utils.py @@ -4,8 +4,39 @@ import matplotlib.pyplot as plt from sklearn.mixture import GaussianMixture -def readData(file, n, m=8, minT=800, maxT=4000): - """Read data from root file +def readFileData(file, n=6, m=8, minT=800, maxT=4000): + """Read whole data from root file + + Parameters + ---------- + file : str + root file path + n : int, optional + number of blocks, default 6 + m : int, optional + number of binds, default 8 + minT/maxT : int, optional + Filtering data, the sum of the left and right sides needs to be in the interval [minT, maxT] + min / max threshold + """ + data = uproot.open(file)["Tree1"] + + ldata, rdata = [], [] + for i in range(n): + 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") + idx = np.where((x + y >= minT) & (x + y <= maxT))[0] + ldata.append(x[idx]) + rdata.append(y[idx]) + + return ldata, rdata + + +def readBlockData(file, n, m=8, minT=800, maxT=4000): + """Read block data from root file Parameters ---------- @@ -25,8 +56,8 @@ def readData(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() - y = data["adc{:d}ch{:d}".format(na, nc + m)].array() + 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") idx = np.where((x + y >= minT) & (x + y <= maxT))[0] ldata.append(x[idx]) rdata.append(y[idx])