fix: line fit formula, b1+b2 calculation;add: data process

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.get_kb(bias, deltaE)
-        # bind.draw_fit_line("result/FIT-LINE/" + bind.name + ".png")
+        bind.draw_fit_line("result/FIT-LINE/" + bind.name + ".png")
-        # bind.draw_cluster("result/GMM/" + bind.name + ".png")
+        bind.draw_cluster("result/GMM/" + bind.name + ".png")
-        # bind.draw_peak("result/PEAK/" + 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
             )
         )

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()

qdx/Bind.py

@@ -81,16 +81,56 @@ class Bind(object):
         self.y[1] = data[:, 1]
         self.px[1] = data[:, 2]
 
-    def get_energy(self):
+    def get_line(self):
-        """Get energy (in channel) by fit Gaussian Model with ldata + rdata"""
+        """Fit data with $x = -\\frac{k_2}{k_1}y + \\frac{E-b_1-b_2}{k_1}$."""
-        eng = self.x[0] + self.y[0]
+        model = Linear1D()
-        fitted_model = fit_hist_gaussian(eng)
+
+        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]
+        eng = self.k1 * self.x[1] + self.k2 * self.y[1]
-        fitted_model = fit_hist_gaussian(eng)
+        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)
+        E = np.mean(self.k1 * self.fx + self.k2 * self.fy + self.b)
-        reg = fit_line(model, np.array(list(zip(self.fx, self.fy)), dtype=object), self.fz)
+        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]

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)
 

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)
+        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_E]
+        return [d_L, d_C, d_k1, d_k2, d_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):
+def readFileData(file, n=6, m=8, minT=800, maxT=4000):
-    """Read data from root file
+    """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()
+        x = data["adc{:d}ch{:d}".format(na, nc)].array(library="np")
-        y = data["adc{:d}ch{:d}".format(na, nc + m)].array()
+        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])