feat: data process script

This commit is contained in:
liuyihui 2023-06-01 15:15:10 +08:00
parent c6fab20c6c
commit 0397ef71f8
20 changed files with 653 additions and 22 deletions

1
.gitignore vendored
View File

@ -15,3 +15,4 @@ test/
# Data
*.root
*.data
data/

View File

@ -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

View File

@ -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));
// }
// }
}

View File

@ -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

BIN
script/draw Executable file

Binary file not shown.

76
script/draw.cpp Normal file
View File

@ -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;
}

127
script/draw.py Normal file
View File

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

32
script/eff.py Normal file
View File

@ -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")

80
script/eff.txt Normal file
View File

@ -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

24
script/eff2.py Normal file
View File

@ -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")

15
script/eff2.txt Normal file
View File

@ -0,0 +1,15 @@
5171
5605
5750
5789
5649
5669
5720
5642
5492
5353
5165
5030
4957
4831
4642

24
script/eff3.py Normal file
View File

@ -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")

21
script/ej-200.txt Normal file
View File

@ -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

15
script/inner.txt Normal file
View File

@ -0,0 +1,15 @@
2295
2811
3105
3346
3515
3543
3887
4070
4253
4331
3835
3960
4153
4236
4213

16
script/main.py Normal file
View File

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

15
script/outer.txt Normal file
View File

@ -0,0 +1,15 @@
2876
2794
2645
2443
2134
2126
1833
1572
1239
1022
1330
1070
804
595
429

27
script/roi.py Normal file
View File

@ -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")

80
script/roi.txt Normal file
View File

@ -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

28
script/roi2.py Normal file
View File

@ -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")

15
script/roi2.txt Normal file
View File

@ -0,0 +1,15 @@
3798
4083
4283
4214
4161
4194
4125
4138
3990
3903
3745
3638
3624
3532
3410