diff --git a/.gitignore b/.gitignore
index d58425e..be48023 100644
--- a/.gitignore
+++ b/.gitignore
@@ -15,3 +15,4 @@ test/
 # Data
 *.root
 *.data
+data/
diff --git a/init_vis.mac b/init_vis.mac
index 2ed5342..b7f97a6 100644
--- a/init_vis.mac
+++ b/init_vis.mac
@@ -1,7 +1,9 @@
 /control/verbose 1
-/run/verbose 2
-/event/verbose 0
-/tracking/verbose 0
+/run/verbose 1
+/event/verbose 1
+/tracking/verbose 1
+
+/Gene/PrimaryGA/SetBeamEnergy 1 MeV
 
 # Initialize kernel
 /run/initialize
diff --git a/integral.cpp b/integral.cpp
index 15380fc..5008b91 100644
--- a/integral.cpp
+++ b/integral.cpp
@@ -1,39 +1,72 @@
 #include <TFile.h>
+#include <TString.h>
 #include <TTree.h>
 #include <stdio.h>
 
-void integral(TString file) {
+int integral(TString file) {
     double e[24];
     int id[24], N[1];
-    int sum1, sum2;
+    int sum1, sum2, sum3;
 
     TFile f(file);
-    if (f.IsZombie()) {
-        printf("cannot open root file\n");
-        return;
-    }
     TTree* t = (TTree*)f.Get("SimData1");
-    if (!t) {
-        printf("cannot find tree!\n");
-        return;
-    }
     t->SetBranchAddress("detE", e);
     t->SetBranchAddress("detId", id);
     t->SetBranchAddress("Ndets", N);
     int n = t->GetEntries();
     //	printf("n=%d\n", n);
     int i, j;
-    sum1 = sum2 = 0;
+    sum1 = sum2 = sum3 = 0;
     for (i = 0; i < n; i++) {
         t->GetEntry(i);
         for (j = 0; j < N[0]; j++) {
-            if (e[j] > 0.18 && id[j] < 12) {
-                sum1++;
-            }
-            if (e[j] > 0.18 && id[j] > 11) {
-                sum2++;
-            }
+            if (e[j] > 0.18 && id[j] < 12) sum1++;
+            if (e[j] > 0.18 && id[j] > 11) sum2++;
+            if (e[j] >= 0.64 && e[j] <= 0.8) sum3++;
         }
     }
-    printf("%d,%d,%d,%f\n", sum1 + sum2, sum1, sum2, 1.0 * sum1 / sum2);
+    // printf("%d, %d, %d, %d, %f\n", n, sum1 + sum2, sum1, sum2, 1.0 * sum1 / sum2);
+    // return sum1 + sum2;
+    return sum3;
+}
+
+void process() {
+    double or1[5] = {10, 11, 12, 13, 14};
+    double or2[5] = {12.14, 13.35, 14.57, 15.79, 17};
+    double or3[5] = {14.28, 15.71, 17.14, 18.57, 20};
+    // int energy[8] = {600, 650, 700, 750, 800, 850, 900, 950};
+
+    for (int i = 0; i < 5; i++) {
+        TString file =
+            TString::Format("data/7/%d/Simulation_RunNo_650kev_70mm_%dmm.root", int(or1[i]), int(or1[i] * 10));
+        printf("%d\n", integral(file));
+    }
+
+    for (int i = 0; i < 5; i++) {
+        TString file =
+            TString::Format("data/8.5/%.2f/Simulation_RunNo_650kev_85mm_%dmm.root", or2[i], int(or2[i] * 10));
+        printf("%d\n", integral(file));
+    }
+
+    for (int i = 0; i < 5; i++) {
+        TString file =
+            TString::Format("data/10/%.2f/Simulation_RunNo_650kev_100mm_%dmm.root", or3[i], int(or3[i] * 10));
+        printf("%d\n", integral(file));
+    }
+
+    // for (int i = 0; i < 5; i++) {
+    //     for (int j = 0; j < 8; j++) {
+    //         TString file = TString::Format("data/7/%d/Simulation_RunNo_%dkev_70mm_%dmm.root", int(or1[i]), energy[j],
+    //                                        int(or1[i] * 10));
+    //         printf("%d\n", integral(file));
+    //     }
+    // }
+
+    // for (int i = 0; i < 5; i++) {
+    //     for (int j = 0; j < 8; j++) {
+    //         TString file = TString::Format("data/8.5/%.2f/Simulation_RunNo_%dkev_85mm_%dmm.root", or2[i], energy[j],
+    //                                        int(or2[i] * 10));
+    //         printf("%d\n", integral(file));
+    //     }
+    // }
 }
diff --git a/run.mac b/run.mac
index f238e54..dc6d816 100644
--- a/run.mac
+++ b/run.mac
@@ -6,4 +6,4 @@
 /run/initialize
 
 /gun/particle neutron
-/control/loop loop.mac energy 0.1 5 0.1
+/control/loop loop.mac energy 0.5 1 0.05
diff --git a/script/draw b/script/draw
new file mode 100755
index 0000000..7b999d4
Binary files /dev/null and b/script/draw differ
diff --git a/script/draw.cpp b/script/draw.cpp
new file mode 100644
index 0000000..b853e36
--- /dev/null
+++ b/script/draw.cpp
@@ -0,0 +1,76 @@
+#include "math.h"
+#include "stdlib.h"
+
+#include "iostream"
+
+const double pi = 3.1415926535;
+
+double mass[4];
+double vbeta, vgamma;
+double t3_cm;
+double t4_cm;
+double p3_cm;
+double p4_cm;
+
+void SetBeamEnergy(double energy) {
+    double umass = 931.494;
+    double massHe4 = 4.00260325413 * umass;
+    double massNe22 = 21.991385114 * umass;
+    double massMg25 = 24.98583697 * umass;
+    double massn = 939.56563;
+
+    mass[0] = massHe4;
+    mass[1] = massNe22;
+    mass[2] = massMg25;
+    mass[3] = massn;
+
+    double eb = energy + mass[0];
+    double pb2 = energy * energy + 2.0 * energy * mass[0];
+    double pb = std::sqrt(pb2);
+    double esum = eb + mass[1];
+    vbeta = pb / esum;
+    vgamma = 1.0 / std::sqrt(1.0 - vbeta * vbeta);
+    double e_cm2 = esum * esum - pb2;
+    double e_cm = std::sqrt(e_cm2);
+    double t_cm = e_cm - mass[2] - mass[3];
+    if (t_cm < 0.0) {
+        std::cout << "No solution!" << std::endl;
+        return;
+    }
+    double t_cm2 = t_cm * t_cm;
+    t3_cm = (t_cm2 + 2.0 * mass[3] * t_cm) / e_cm / 2.0;
+    t4_cm = (t_cm2 + 2.0 * mass[2] * t_cm) / e_cm / 2.0;
+    double p3_cm2 = t3_cm * t3_cm + 2.0 * t3_cm * mass[2];
+    double p4_cm2 = t4_cm * t4_cm + 2.0 * t4_cm * mass[3];
+    p3_cm = std::sqrt(p3_cm2);
+    p4_cm = std::sqrt(p4_cm2);
+}
+
+void GeneratePrimaries(double thetacmr) {
+    double thetalabr;
+
+    double tg_thetalabr =
+        p4_cm * std::sin(thetacmr) /
+        (vgamma * (p4_cm * std::cos(thetacmr) + vbeta * std::sqrt(p4_cm * p4_cm + mass[3] * mass[3])));
+    if (tg_thetalabr >= 1.0e6)
+        thetalabr = pi / 2.0;
+    else
+        thetalabr = std::atan(tg_thetalabr);
+    if (thetalabr < 0.0) thetalabr = pi + thetalabr;
+
+    double p4_cmx = p4_cm * std::sin(thetacmr);
+    double p4_cmz = p4_cm * std::cos(thetacmr);
+    double p4_labx = p4_cmx;
+    double p4_labz = vgamma * (p4_cmz + vbeta * (t4_cm + mass[3]));
+    double p4_lab = std::sqrt(p4_labx * p4_labx + p4_labz * p4_labz);
+    double e4_lab = sqrt(p4_lab * p4_lab + mass[3] * mass[3]) - mass[3];
+
+    // std::cout << thetalabr << ", " << e4_lab << std::endl;
+    std::cout << e4_lab << std::endl;
+}
+
+int main() {
+    SetBeamEnergy(1);
+    for (int i = 0; i < 100; i++) GeneratePrimaries(i * pi / 100.0);
+    return 0;
+}
diff --git a/script/draw.py b/script/draw.py
new file mode 100644
index 0000000..7526c8e
--- /dev/null
+++ b/script/draw.py
@@ -0,0 +1,127 @@
+import numpy as np
+from matplotlib import pyplot as plt
+
+mu = 931.4940954  # MeV/c^2
+
+
+class Particle(object):
+    """Define a Particle
+
+    Parameters
+    ----------
+    A : int
+        mass number
+    m : float
+        mass
+    """
+
+    def __init__(self, A, m) -> None:
+        self.A = A
+        self.m = m
+
+
+class Reaction(object):
+    """A Nuclear Reaction A(a,b)B
+    target nuclei (A) is assumed to be quiescent
+
+    Parameters
+    ----------
+    a/A/b/B : Particle
+        Particle a, A, b, B
+    """
+
+    def __init__(self, a, A, b, B) -> None:
+        self.a: Particle = a
+        self.A: Particle = A
+        self.b: Particle = b
+        self.B: Particle = B
+
+    @property
+    def Q(self):
+        """Returns Q value of the reaction"""
+        return (self.A.m + self.a.m - self.B.m - self.b.m) * mu
+
+    @property
+    def threshold(self):
+        """Returns threshold energy of the reaction"""
+        if self.Q < 0:
+            return -(self.A.m + self.a.m) / self.A.m * self.Q
+        return 0
+
+    def eb(self, theta, ke):
+        """Returns energy of b in the lab frame
+
+        Parameters
+        ----------
+        theta : float
+            angle between a and b in the lab frame in rad
+        ke : float
+            kinetic energy of a in MeV
+        """
+        i1 = np.sqrt(self.a.A * self.b.A * ke) / (self.B.A + self.b.A) * np.cos(theta)
+        i2 = (self.B.A - self.a.A) / (self.B.A + self.b.A) + (self.a.A * self.b.A) / ((self.B.A + self.b.A) ** 2) * (np.cos(theta) ** 2)
+        i3 = self.B.A * self.Q / (self.B.A + self.b.A)
+        return (i1 + np.sqrt(i2 * ke + i3)) ** 2
+
+    def eB(self, theta, ke):
+        """Returns energy of B in the lab frame
+
+        Parameters
+        ----------
+        theta : float
+            angle between a and b in the lab frame in rad
+        ke : float
+            kinetic energy of a in MeV
+        """
+        eb = self.eb(theta, ke)
+        i1 = self.a.A / self.B.A * ke + self.b.A / self.B.A * eb
+        i2 = 2 * np.sqrt(self.a.A * self.b.A * ke * eb) / self.B.A
+        return i1 - i2 * np.cos(theta)
+
+
+neutron = Particle(1, 1.008665)
+proton = Particle(1, 1.007825)
+tritium = Particle(3, 3.016049)
+he_3 = Particle(3, 3.016029)
+he_4 = Particle(4, 4.002603)
+ne_22 = Particle(22, 21.991385)
+mg_25 = Particle(25, 24.985837)
+
+theta = np.linspace(0, np.pi, 100)
+
+r1 = Reaction(he_4, ne_22, neutron, mg_25)
+res = r1.eb(theta, 1)
+
+plt.plot(theta, res)
+plt.title(r"$^{22}$Ne($\alpha$,n)$^{25}$Mg, $T_\alpha$=0.8 MeV")
+plt.ylabel(r"$T_n$ (MeV)")
+plt.xlabel(r"$\theta$ (rad)")
+plt.xticks([0, np.pi / 4, np.pi / 2, np.pi * 3 / 4, np.pi], [r"$0$", r"$\pi/4$", r"$\pi/2$", r"$3\pi/4$", r"$\pi$"])
+plt.xlim(0, np.pi)
+plt.show()
+
+exit(0)
+
+energy = np.linspace(np.min(res), np.max(res), 100)
+
+r2 = Reaction(he_3, neutron, proton, tritium)
+
+# res = r2.eb(theta, 0.0114)
+# res2 = r2.eB(theta, 0.0114)
+
+res = r2.eb(np.pi / 4, energy)
+res2 = r2.eB(np.pi / 4, energy)
+
+fig, ax = plt.subplots()
+
+ax.plot(energy, res, color="C0")
+ax.set_title(r"$^3$He(n,p)$^3$H, $\theta$=$\pi/4$")
+ax.set_xlabel(r"$T_n$ (MeV)")
+ax.set_ylabel(r"$T_p$ (MeV)")
+
+ax1 = ax.twinx()
+ax1.plot(energy, res2, color="C1")
+ax1.set_ylabel(r"$T_t$ (MeV)")
+
+plt.tight_layout()
+plt.show()
diff --git a/script/eff.py b/script/eff.py
new file mode 100644
index 0000000..edff3da
--- /dev/null
+++ b/script/eff.py
@@ -0,0 +1,32 @@
+import numpy as np
+from matplotlib import pyplot as plt
+
+or1 = [10, 11, 12, 13, 14]
+or2 = [12.14, 13.35, 14.57, 15.79, 17]
+energy = [600, 650, 700, 750, 800, 850, 900, 950]
+markers = ["^", "+", "s", "D", "."]
+
+data = np.loadtxt("test/eff.txt")
+data = data.reshape(10, 8)
+data = data / 10000
+
+plt.figure(dpi=400)
+
+ax = plt.subplot(1, 1, 1)
+for k in range(5):
+    ax.scatter(energy, data[k], label="{:d}cm".format(or1[k]), marker=markers[k])
+ax.legend(loc="lower left")
+ax.set_xlabel("Energy [keV]", fontsize=15)
+ax.set_ylabel("Efficiency", fontsize=15)
+
+# ax = plt.subplot(2, 1, 2)
+# for k in range(5):
+#     ax.scatter(energy, data[k + 5], label="{:.2f}cm".format(or2[k]), marker=markers[k])
+# ax.legend(loc="lower left")
+# ax.set_xlabel("Energy [keV]", fontsize=15)
+# ax.set_ylabel("Efficiency", fontsize=15)
+
+plt.tight_layout()
+plt.subplots_adjust()
+# plt.show()
+plt.savefig("/home/fox/Documents/Thesis/figures/chap3/Fig3.png")
diff --git a/script/eff.txt b/script/eff.txt
new file mode 100644
index 0000000..9fa32b8
--- /dev/null
+++ b/script/eff.txt
@@ -0,0 +1,80 @@
+5389
+5171
+4988
+4774
+4722
+4599
+4533
+4401
+5726
+5605
+5383
+5255
+5089
+5022
+4959
+4824
+5880
+5750
+5613
+5496
+5341
+5233
+5190
+5094
+5885
+5789
+5669
+5515
+5480
+5345
+5317
+5121
+5751
+5649
+5455
+5535
+5436
+5266
+5141
+5253
+5766
+5669
+5512
+5378
+5405
+5283
+5127
+5085
+5819
+5720
+5632
+5544
+5472
+5485
+5456
+5162
+5650
+5642
+5507
+5456
+5437
+5366
+5368
+5361
+5567
+5492
+5435
+5347
+5350
+5248
+5275
+5182
+5412
+5353
+5310
+5181
+5219
+5121
+5061
+5055
\ No newline at end of file
diff --git a/script/eff2.py b/script/eff2.py
new file mode 100644
index 0000000..498a1a7
--- /dev/null
+++ b/script/eff2.py
@@ -0,0 +1,24 @@
+import numpy as np
+from matplotlib import pyplot as plt
+
+or1 = [10, 11, 12, 13, 14]
+markers = ["^", "+", "s", "D", "."]
+
+data = np.loadtxt("test/eff2.txt")
+data = data.reshape(3, 5)
+data = data / 10000
+
+fig = plt.figure(dpi=400)
+
+ax = plt.subplot(1, 1, 1)
+for k in range(5):
+    ax.scatter([7, 8.5, 10], data[:, k], label=r"$R_o/R_i$ = {:.2f}".format(or1[k] / 7), marker=markers[k])
+ax.legend(loc="lower left")
+ax.set_xlabel(r"$R_i$ [cm]", fontsize=15)
+ax.set_ylabel("Efficiency", fontsize=15)
+ax.set_xticks([7, 8.5, 10])
+
+plt.tight_layout()
+plt.subplots_adjust()
+# plt.show()
+plt.savefig("/home/fox/Documents/Thesis/figures/chap3/Fig5.png")
diff --git a/script/eff2.txt b/script/eff2.txt
new file mode 100644
index 0000000..7909d6d
--- /dev/null
+++ b/script/eff2.txt
@@ -0,0 +1,15 @@
+5171
+5605
+5750
+5789
+5649
+5669
+5720
+5642
+5492
+5353
+5165
+5030
+4957
+4831
+4642
\ No newline at end of file
diff --git a/script/eff3.py b/script/eff3.py
new file mode 100644
index 0000000..2d4dfdc
--- /dev/null
+++ b/script/eff3.py
@@ -0,0 +1,24 @@
+import numpy as np
+from matplotlib import pyplot as plt
+
+or1 = [10, 11, 12, 13, 14]
+markers = ["^", "+", "s", "D", "."]
+
+data = np.loadtxt("test/inner.txt")
+data = data.reshape(3, 5)
+data = data / 10000
+
+fig = plt.figure(dpi=400)
+
+ax = plt.subplot(1, 1, 1)
+for k in range(5):
+    ax.scatter([7, 8.5, 10], data[:, k], label=r"$R_o/R_i$ = {:.2f}".format(or1[k] / 7), marker=markers[k])
+ax.legend(loc="lower right")
+ax.set_xlabel(r"$R_i$ [cm]", fontsize=15)
+ax.set_ylabel("Efficiency", fontsize=15)
+ax.set_xticks([7, 8.5, 10])
+
+plt.tight_layout()
+plt.subplots_adjust()
+# plt.show()
+plt.savefig("/home/fox/Documents/Thesis/figures/chap3/Fig6(a).png")
diff --git a/script/ej-200.txt b/script/ej-200.txt
new file mode 100644
index 0000000..c5188ae
--- /dev/null
+++ b/script/ej-200.txt
@@ -0,0 +1,21 @@
+395.00000, 0.00000
+400.00000, 0.00926
+403.19005, 0.03801
+408.51147, 0.19876
+410.83959, 0.30444
+412.75197, 0.40119
+418.82171, 0.79712
+421.48241, 0.93703
+422.81277, 0.97870
+425.00000, 1.00000
+427.63530, 0.98168
+433.37245, 0.89833
+440.10736, 0.74651
+446.75913, 0.60064
+460.06266, 0.43394
+465.71666, 0.37440
+472.45158, 0.26426
+480.10111, 0.17197
+490.16190, 0.10053
+499.97326, 0.05885
+515.00000, 0.00000
\ No newline at end of file
diff --git a/script/inner.txt b/script/inner.txt
new file mode 100644
index 0000000..df7a5aa
--- /dev/null
+++ b/script/inner.txt
@@ -0,0 +1,15 @@
+2295
+2811
+3105
+3346
+3515
+3543
+3887
+4070
+4253
+4331
+3835
+3960
+4153
+4236
+4213
\ No newline at end of file
diff --git a/script/main.py b/script/main.py
new file mode 100644
index 0000000..e3eb57a
--- /dev/null
+++ b/script/main.py
@@ -0,0 +1,16 @@
+import numpy as np
+
+data = np.loadtxt("ej-200.txt", delimiter=",")
+
+c = 299792458
+h = 6.62606957 * 1e-34
+hc = h * c / (1.6 * 1e-19)
+
+L = len(data)
+for k in range(L):
+    print(str(round(hc * 1e9 / data[L - k - 1][0], 3)) + " * eV", end=", ")
+print()
+
+for k in range(L):
+    print("{:.5f}".format(round(data[L - k - 1][1], 5)), end=", ")
+print()
diff --git a/script/outer.txt b/script/outer.txt
new file mode 100644
index 0000000..e05c2ca
--- /dev/null
+++ b/script/outer.txt
@@ -0,0 +1,15 @@
+2876
+2794
+2645
+2443
+2134
+2126
+1833
+1572
+1239
+1022
+1330
+1070
+804
+595
+429
\ No newline at end of file
diff --git a/script/roi.py b/script/roi.py
new file mode 100644
index 0000000..8e7619c
--- /dev/null
+++ b/script/roi.py
@@ -0,0 +1,27 @@
+import numpy as np
+from matplotlib import pyplot as plt
+
+or1 = [10, 11, 12, 13, 14]
+or2 = [12.14, 13.35, 14.57, 15.79, 17]
+energy = [600, 650, 700, 750, 800, 850, 900, 950]
+markers = ["^", "+", "s", "D", "."]
+
+data = np.loadtxt("test/eff.txt")
+data = data.reshape(10, 8)
+
+roi = np.loadtxt("test/roi.txt")
+roi = roi.reshape(10, 8)
+roi = 1 - roi / data
+
+plt.figure(dpi=400)
+
+for k in range(5):
+    plt.scatter(energy, roi[k], label="{:d}cm".format(or1[k]), marker=markers[k])
+plt.xlabel("Energy [keV]", fontsize=15)
+plt.ylabel(r"Energy Leakage Rate $\epsilon$", fontsize=15)
+plt.legend(loc="lower left")
+
+plt.tight_layout()
+plt.subplots_adjust()
+# plt.show()
+plt.savefig("/home/fox/Documents/Thesis/figures/chap3/Fig7(a).png")
diff --git a/script/roi.txt b/script/roi.txt
new file mode 100644
index 0000000..0812ed9
--- /dev/null
+++ b/script/roi.txt
@@ -0,0 +1,80 @@
+3973
+3798
+3656
+3505
+3442
+3391
+3363
+3277
+4219
+4083
+3952
+3875
+3694
+3634
+3680
+3532
+4233
+4283
+4186
+4057
+3921
+3818
+3758
+3745
+4342
+4214
+4178
+4024
+3992
+3954
+3888
+3721
+4202
+4161
+4027
+4012
+3971
+3855
+3760
+3819
+4252
+4194
+4073
+3957
+3945
+3873
+3761
+3638
+4296
+4125
+4112
+4117
+4037
+3986
+3951
+3748
+4144
+4138
+4041
+4001
+3965
+3924
+3951
+3880
+4067
+3990
+3938
+3902
+3933
+3812
+3806
+3772
+3922
+3903
+3913
+3755
+3813
+3772
+3612
+3724
\ No newline at end of file
diff --git a/script/roi2.py b/script/roi2.py
new file mode 100644
index 0000000..f86db10
--- /dev/null
+++ b/script/roi2.py
@@ -0,0 +1,28 @@
+import numpy as np
+from matplotlib import pyplot as plt
+
+or1 = [10, 11, 12, 13, 14]
+or2 = [12.14, 13.35, 14.57, 15.79, 17]
+energy = [600, 650, 700, 750, 800, 850, 900, 950]
+markers = ["^", "+", "s", "D", "."]
+
+data = np.loadtxt("test/eff2.txt")
+data = data.reshape(3, 5)
+
+roi = np.loadtxt("test/roi2.txt")
+roi = roi.reshape(3, 5)
+roi = 1 - roi / data
+
+plt.figure(dpi=400)
+
+for k in range(5):
+    plt.scatter([7, 8.5, 10], roi[:, k], label=r"$R_o/R_i$ = {:.2f}".format(or1[k] / 7), marker=markers[k])
+plt.xlabel(r"$R_i$ [cm]", fontsize=15)
+plt.xticks([7, 8.5, 10])
+plt.ylabel(r"Energy Leakage Rate $\epsilon$", fontsize=15)
+plt.legend(loc="lower right")
+
+plt.tight_layout()
+plt.subplots_adjust()
+# plt.show()
+plt.savefig("/home/fox/Documents/Thesis/figures/chap3/Fig7(b).png")
diff --git a/script/roi2.txt b/script/roi2.txt
new file mode 100644
index 0000000..d466fe3
--- /dev/null
+++ b/script/roi2.txt
@@ -0,0 +1,15 @@
+3798
+4083
+4283
+4214
+4161
+4194
+4125
+4138
+3990
+3903
+3745
+3638
+3624
+3532
+3410
\ No newline at end of file