add: simulation result, Particle Definition;

change: model format;
add: template class for PrimaryGeneratorAction;
add: random algorithm in PrimaryGeneratorAction;
This commit is contained in:
liuyihui 2022-05-16 18:08:47 +08:00
parent 3b2b20209d
commit f09c1d4407
16 changed files with 358 additions and 96 deletions

View File

@ -5,9 +5,9 @@ project(DESCSS)
option(WITH_GEANT4_UIVIS "Build example with Geant4 UI and Vis drivers" ON)
if(WITH_GEANT4_UIVIS)
find_package(Geant4 REQUIRED ui_all vis_all)
find_package(Geant4 REQUIRED ui_all vis_all)
else()
find_package(Geant4 REQUIRED)
find_package(Geant4 REQUIRED)
endif()
include(${Geant4_USE_FILE})
@ -20,7 +20,10 @@ add_executable(DESCSS main.cpp ${sources} ${headers})
target_link_libraries(DESCSS ${Geant4_LIBRARIES})
set(DESCSS_SCRIPTS
debug.mac
electron.mac
vis.mac
assets/model.txt
)
foreach(_script ${DESCSS_SCRIPTS})

View File

@ -105,7 +105,12 @@
"wcsv_ntuple": "cpp",
"bufobj": "cpp",
"aidas": "cpp",
"element": "cpp"
"element": "cpp",
"colorf": "cpp",
"cid": "cpp",
"mem": "cpp",
"sout": "cpp",
"typedefs": "cpp"
}
}
}

View File

@ -108,7 +108,7 @@ $$
两者同时进行拟合与比较,同时考虑简洁性,除`H`以外,其余元素均采用双指数模型,`H`采用简化模型。
#### 俘获辐射带
1. 电子[^5]
1. 电子[^5]
在轨道较高($L>2.5$)的内带中,电子能谱中有$\sim64\%$符合指数分布,$\sim4\%$符合幂律分布。
但在本例中,轨道较低($L\sim1.06$),因此使用指数模型$j = j_0e^{-\frac{E}{E_0}}$和幂律模型$j=j_0E^{-\alpha}$进行拟合,随后使用`Origin Lab`对两者进行`BIC test`,最终选用指数模型$j_0=1.333680\times10^6,\ E_0=0.0824$。
@ -124,7 +124,7 @@ $$
<div align=center><img src="docs/ExFit-trapped-electron.png" style="max-width: 80%;"><br /></div>
2. 质子
2. 质子
与电子保持一致,幂律模型修改为$j=j_0(E-\beta)^{-\alpha}$,进行`BIC test`,最终选用幂律模型$j_0=559.76377,\ \alpha=2.40873,\ \beta=-1.53424$。
<div align=center>
@ -154,34 +154,34 @@ $$
|:--------:|:--------------------------------------------------------------------:|:-------:|:-----------------------:|:-----------------------:|:-----------------------:|:-------:|:-----------:|
| $e^{-}$ | $j=j_0e^{-\frac{E}{E_0}}$ | 0.98107 | 0.1 | 10 | 9.23242E+11 | 923242 | 999999.8096 |
| `p` | $j=j_0(E-\beta)^{-\alpha}$ | 0.99883 | 50 | 1000 | 42853994.15 | 100000 | 428.5399415 |
| `H` | $\frac{C_1\beta^{\alpha-1}}{pc}\left(\frac{pc}{pc+C_2}\right)^{C_3}$ | 0.98915 | 200 | 100000 | 12854310.92 | 100000 | 128.5431092 |
| `He` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98936 | 200 | 100000 | 1718828.12 | 100000 | 17.1882812 |
| `Li` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99755 | 200 | 100000 | 12982.36284 | 100000 | 0.129823628 |
| `Be` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99708 | 200 | 100000 | 6687.159271 | 100000 | 0.066871593 |
| `B` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99569 | 200 | 100000 | 17093.65985 | 100000 | 0.170936599 |
| `C` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98936 | 200 | 100000 | 52245.07239 | 100000 | 0.522450724 |
| `N` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98892 | 200 | 100000 | 13564.66838 | 100000 | 0.135646684 |
| `O` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98936 | 200 | 100000 | 48807.93662 | 100000 | 0.488079366 |
| `F` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99323 | 200 | 100000 | 1122.681148 | 100000 | 0.011226811 |
| `Ne` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99014 | 200 | 100000 | 8065.236798 | 100000 | 0.080652368 |
| `Na` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99263 | 200 | 100000 | 1865.503552 | 100000 | 0.018655036 |
| `Mg` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98972 | 200 | 100000 | 10928.24228 | 100000 | 0.109282423 |
| `Al` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99217 | 200 | 100000 | 1938.943418 | 100000 | 0.019389434 |
| `Si` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98873 | 200 | 100000 | 8262.870777 | 100000 | 0.082628708 |
| `P` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99183 | 200 | 100000 | 413.5749377 | 100000 | 0.004135749 |
| `S` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98873 | 200 | 100000 | 1659.707 | 100000 | 0.01659707 |
| `Cl` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99087 | 200 | 100000 | 371.7198798 | 100000 | 0.003717199 |
| `Ar` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99432 | 200 | 100000 | 751.2576633 | 100000 | 0.007512577 |
| `K` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98963 | 200 | 100000 | 466.3178791 | 100000 | 0.004663179 |
| `Ca` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98935 | 200 | 100000 | 1161.134206 | 100000 | 0.011611342 |
| `Sc` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99251 | 200 | 100000 | 231.3722059 | 100000 | 0.002313722 |
| `Ti` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.9934 | 200 | 100000 | 828.6605823 | 100000 | 0.008286606 |
| `V` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99424 | 200 | 100000 | 403.8859746 | 100000 | 0.00403886 |
| `Cr` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99313 | 200 | 100000 | 783.4909673 | 100000 | 0.00783491 |
| `Mn` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99385 | 200 | 100000 | 570.5548534 | 100000 | 0.005705549 |
| `Fe` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99351 | 200 | 100000 | 6024.476681 | 100000 | 0.060244767 |
| `Co` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99426 | 200 | 100000 | 20.89166947 | 100000 | 0.000208917 |
| `Ni` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99231 | 200 | 100000 | 292.8131737 | 100000 | 0.002928132 |
| `H` | $\frac{C_1\beta^{\alpha-1}}{pc}\left(\frac{pc}{pc+C_2}\right)^{C_3}$ | 0.98915 | 220 | 100000 | 12854310.92 | 100000 | 128.5431092 |
| `He` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98936 | 220 | 100000 | 1720188.128 | 100000 | 17.20188128 |
| `Li` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99755 | 220 | 100000 | 12988.97006 | 100000 | 0.129889701 |
| `Be` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99708 | 220 | 100000 | 6691.176848 | 100000 | 0.066911768 |
| `B` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99569 | 220 | 100000 | 17105.83698 | 100000 | 0.17105837 |
| `C` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98936 | 220 | 100000 | 52286.55444 | 100000 | 0.522865544 |
| `N` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98892 | 220 | 100000 | 13576.5542 | 100000 | 0.135765542 |
| `O` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98936 | 220 | 100000 | 48846.63508 | 100000 | 0.488466351 |
| `F` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99323 | 220 | 100000 | 1123.492483 | 100000 | 0.011234925 |
| `Ne` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99014 | 220 | 100000 | 8071.585382 | 100000 | 0.080715854 |
| `Na` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99263 | 220 | 100000 | 1866.8702 | 100000 | 0.018668702 |
| `Mg` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98972 | 220 | 100000 | 10936.4597 | 100000 | 0.109364597 |
| `Al` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99217 | 220 | 100000 | 1940.381408 | 100000 | 0.019403814 |
| `Si` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98873 | 220 | 100000 | 8269.257173 | 100000 | 0.082692572 |
| `P` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99183 | 220 | 100000 | 413.8897994 | 100000 | 0.004138898 |
| `S` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98873 | 220 | 100000 | 1660.990128 | 100000 | 0.016609901 |
| `Cl` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99087 | 220 | 100000 | 372.0271181 | 100000 | 0.003720271 |
| `Ar` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99432 | 220 | 100000 | 751.7977189 | 100000 | 0.007517977 |
| `K` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98963 | 220 | 100000 | 466.7118907 | 100000 | 0.004667119 |
| `Ca` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.98935 | 220 | 100000 | 1161.900661 | 100000 | 0.011619007 |
| `Sc` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99251 | 220 | 100000 | 231.5541271 | 100000 | 0.002315541 |
| `Ti` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.9934 | 220 | 100000 | 829.2810844 | 100000 | 0.008292811 |
| `V` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99424 | 220 | 100000 | 404.1724078 | 100000 | 0.004041724 |
| `Cr` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99313 | 220 | 100000 | 784.0902938 | 100000 | 0.007840903 |
| `Mn` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99385 | 220 | 100000 | 570.9732265 | 100000 | 0.005709732 |
| `Fe` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99351 | 220 | 100000 | 6028.059963 | 100000 | 0.0602806 |
| `Co` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99426 | 220 | 100000 | 20.90356212 | 100000 | 0.000209036 |
| `Ni` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 0.99231 | 220 | 100000 | 292.9961042 | 100000 | 0.002929961 |
</div>
@ -212,17 +212,9 @@ $$
2. `Acceptance-Rejection Method``ARM`
对于一个在$[a, b]$的分布,设其概率密度函数(PDF)为$F(x)$,首先生成一个均匀分布随机数$X \sim U(x_{min}, x_{max})$,随后独立生成另一个均匀分布随机数$Y \sim U(y_{min}, y_{max})$,如果$Y \le F(X)$,则保留$X$,保留下来的$X$满足$F(x)$的分布。
`ITM`也被称为反演法,关键是需要获取累积分布函数的逆函数,效率较高;`ARM`本质上是一种模拟算法,效率低,但是适应性更广,特别是对于一些复杂分布函数。`GCR`质子和俘获辐射粒子的模型有三类,其中`GCR`质子的函数较为简单,其余两类较为复杂,因此分别采用`ITM`和`ARM`。
`ITM`也被称为反演法,关键是需要获取累积分布函数的逆函数,效率较高;`ARM`本质上是一种模拟算法,效率低,但是适应性更广,特别是对于一些复杂分布函数。
<div align=center>
| $F(E)$ | $\int F(E)\mathrm{d}E$ | $C(x)=\int_a^xF(u)\mathrm{d}u$ | $C^{-1}(y)$ |
|:--------------------------:|:------------------------------------------:|:----------------------------------------------------------------------------:|:---------------------------------------------------------------------------------------------------:|
| $j=j_0(E-\beta)^{-\alpha}$ | $\frac{j_0(E-\beta)^{1-\alpha}}{1-\alpha}$ | $\frac{j_0}{1-\alpha}\left((x-\beta)^{1-\alpha}-(a-\beta)^{1-\alpha}\right)$ | $\sqrt[1-\alpha]{\frac{(y+C)(1-\alpha)}{j_0}}+\beta,\ C=\frac{j_0(a-\beta)^{1-\alpha}}{1-\alpha}$ |
</div>
1. 位置生成
1. 位置生成
而发射位置要求在球面上均匀分布。由于球面的面积元$\mathrm{d}\Omega=\sin\theta\mathrm{d}\phi\mathrm{d}\theta$是$\theta$的函数,如果位置向量$\vec{r}(\rho, \phi, \theta)$的$\phi$和$\theta$是在$\phi\in[0, 2\pi)$和$\theta\in[0, \pi]$上均匀分布,选取的点将会在两极较为密集(即),因此可以采用`ITM`方法生成$\theta$$C^{-1}(z)=\cos^{-1}(z)$。
```c
double u = rand();
@ -234,7 +226,7 @@ $$
double z = rho * cos(theta);
```
2. 发射方向
2. 发射方向
在`Geant4`中,假定粒子的位置矢量为$\vec{r}$,则以$\vec{r}$反方向作为$z$轴正方向建立右手坐标系,粒子的指向为:
$$
P_x = -\sin\theta\cos\phi \\
@ -243,6 +235,21 @@ P_z = -\cos\theta
$$
为了保证发射方向是朝向球面内,因此需要限制$\theta$的最大角度为$90^{\circ}$。
3. 能量
* `GCR`质子的函数较为简单,采用`ITM`;(下表中$j_0$已经归一化)
<div align=center>
| $F(E)$ | $\int F(E)\mathrm{d}E$ | $C(x)=\int_a^xF(u)\mathrm{d}u$ | $C^{-1}(y)$ |
|:--------------------------:|:------------------------------------------:|:----------------------------------------------------------------------------:|:---------------------------------------------------------------------------------------------------:|
| $j=j_0(E-\beta)^{-\alpha}$ | $\frac{j_0(E-\beta)^{1-\alpha}}{1-\alpha}$ | $\frac{j_0}{1-\alpha}\left((x-\beta)^{1-\alpha}-(a-\beta)^{1-\alpha}\right)$ | $\sqrt[1-\alpha]{\frac{(y+C)(1-\alpha)}{j_0}}+\beta,\ C=\frac{j_0(a-\beta)^{1-\alpha}}{1-\alpha}$ |
</div>
<div align=center><img src="docs/simu-tp.png" style="max-width: 50%;"></div>
* 俘获质子的函数较为复杂,使用`ARM`生成`100000`个随机数需要`12`分钟,效率满足要求;
<div align=center><img src="docs/simu-gcrp.png" style="max-width: 50%;"></div>
* 其余俘获辐射粒子同样使用`ARM`,生成`100000`个随机数需要`4.5`分钟。
<div align=center><img src="docs/simu-gcri.png" style="max-width: 50%;"></div>
## 空间站建模
1. 尺寸与分区[^6]
<div align=center><img src="docs/size.webp" style="max-width: 50%;"></div>

View File

@ -1,37 +1,29 @@
# GCR
125.89223, 44.00673, -0.000460574, 112356
0.01815, 0.000297726, 0.00184, 0.49271
0.000184794, 0.000389118, 0.0021, 0.49406
0.0000953693, 0.000386546, 0.002, 0.48222
0.000245835, 0.000383147, 0.00182, 0.45495
0.000551454, 0.000297607, 0.00184, 0.49310
0.000177397, 0.000347376, 0.00161, 0.43464
0.000515269, 0.000297659, 0.00184, 0.49294
0.0000121025, 0.000307953, 0.00208, 0.53260
0.0000855093, 0.000299579, 0.00188, 0.50022
0.0000200236, 0.000306074, 0.00204, 0.52582
0.000107752, 0.000283677, 0.00197, 0.52149
0.0000207575, 0.000304806, 0.00201, 0.52089
0.0000809795, 0.000281074, 0.00191, 0.51232
0.00000442291, 0.000303868, 0.00198, 0.51717
0.0000162654, 0.000281067, 0.00191, 0.51235
0.00000463438, 0.00033962, 0.00178, 0.46921
0.00000956898, 0.000351681, 0.00199, 0.49986
0.00000578307, 0.000336363, 0.00172, 0.45907
0.0000109853, 0.000266092, 0.00154, 0.42213
0.00000291072, 0.000344777, 0.00187, 0.48297
0.0000104814, 0.000348052, 0.00193, 0.49105
0.0000051388, 0.000351427, 0.00199, 0.49888
0.0000098941, 0.000346998, 0.00191, 0.48874
0.00000723814, 0.000349723, 0.00196, 0.49543
0.0000592702, 0.000279954, 0.00175, 0.45382
0.00000020753, 0.000283367, 0.0018, 0.46005
0.00000284117, 0.000275123, 0.00168, 0.44440
# Trapped Electron
# j0 E0
1333680 0.0824
# Trapped Proton
# j0 alpha beta
559.76377 2.40873 1.53424
559.76377, -1.53424, 2.40873, 1.5156509
125.89223, 44.00673, -0.000460574, 112356, 361.781666
0.01815, 0.000297726, 0.00184, 0.49271, 48.4143048
0.000184794, 0.000389118, 0.0021, 0.49406, 0.3655716
0.0000953693, 0.000386546, 0.002, 0.48222, 0.1883217
0.000245835, 0.000383147, 0.00182, 0.45495, 0.4814399
0.000551454, 0.000297607, 0.00184, 0.49310, 1.4715932
0.000177397, 0.000347376, 0.00161, 0.43464, 0.382109
0.000515269, 0.000297659, 0.00184, 0.49294, 1.3747775
0.0000121025, 0.000307953, 0.00208, 0.53260, 0.0316204
0.0000855093, 0.000299579, 0.00188, 0.50022, 0.2271729
0.0000200236, 0.000306074, 0.00204, 0.52582, 0.0525426
0.000107752, 0.000283677, 0.00197, 0.52149, 0.3078042
0.0000207575, 0.000304806, 0.00201, 0.52089, 0.0546116
0.0000809795, 0.000281074, 0.00191, 0.51232, 0.2327364
0.00000442291, 0.000303868, 0.00198, 0.51717, 0.0116488
0.0000162654, 0.000281067, 0.00191, 0.51235, 0.0467482
0.00000463438, 0.00033962, 0.00178, 0.46921, 0.0104706
0.00000956898, 0.000351681, 0.00199, 0.49986, 0.0211592
0.00000578307, 0.000336363, 0.00172, 0.45907, 0.0131355
0.0000109853, 0.000266092, 0.00154, 0.42213, 0.0327014
0.00000291072, 0.000344777, 0.00187, 0.48297, 0.006517
0.0000104814, 0.000348052, 0.00193, 0.49105, 0.0233399
0.0000051388, 0.000351427, 0.00199, 0.49888, 0.0113753
0.0000098941, 0.000346998, 0.00191, 0.48874, 0.022068
0.00000723814, 0.000349723, 0.00196, 0.49543, 0.0160699
0.0000592702, 0.000279954, 0.00175, 0.45382, 0.1696584
0.00000020753, 0.000283367, 0.0018, 0.46005, 0.0005883
0.00000284117, 0.000275123, 0.00168, 0.44440, 0.0082463

View File

@ -1 +0,0 @@
22210

BIN
docs/simu-gcri.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 23 KiB

BIN
docs/simu-gcrp.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 27 KiB

BIN
docs/simu-tp.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 22 KiB

View File

@ -3,12 +3,15 @@
#include "G4VUserActionInitialization.hh"
extern G4String particleType;
class ActionInitialization : public G4VUserActionInitialization {
public:
ActionInitialization();
~ActionInitialization() override;
void Build() const override;
void BuildForMaster() const override;
void Build() const;
};
#endif

69
include/Particle.h Normal file
View File

@ -0,0 +1,69 @@
#ifndef DESCSS_Particle_h
#define DESCSS_Particle_h
#include "G4IonTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTable.hh"
#include "globals.hh"
class Particle {
public:
Particle(G4String name, G4int Z, G4int A) {
if (Z == 1)
ion = G4ParticleTable::GetParticleTable()->FindParticle("proton");
else if (Z == 2)
ion = G4ParticleTable::GetParticleTable()->FindParticle("alpha");
else
ion = G4IonTable::GetIonTable()->GetIon(Z, A, 0);
this->name = name;
this->Z = Z, this->A = A;
};
void setC1(G4double v) { C1 = v; };
void setC2(G4double v) { C2 = v; };
void setC3(G4double v) { C3 = v; };
void setC4(G4double v) { C4 = v; };
void setSum(G4double v) { sum = v; };
virtual G4double pdf(G4double E);
protected:
G4ParticleDefinition* ion;
G4int Z, A;
G4String name;
G4double sum;
G4double C1, C2, C3, C4;
};
class TProton : public Particle {
public:
TProton();
G4double pdf(G4double E) {
G4double C = C1 * pow(50 - C2, 1 - C3) / sum / (1 - C3);
return pow(sum * (C + y) * (1 - C3) / C1, 1 / (1 - C3)) + C2;
}
};
class GCRProton : public Particle {
public:
GCRProton();
G4double pdf(G4double E) {
G4double pc = getPc(E);
G4double beta = getBeta(E);
return C1 * pow(beta, C2 - 1) * pow(pc / (pc + C3), C4) / pc / sum;
}
private:
G4double M = 931.5;
G4double getPc(G4double E) { return sqrt(2 * M * E + E * E); }
G4double getBeta(G4double E) {
G4double pc = getPc(E);
return pc / sqrt(M * M + pc * pc);
}
};
class GCRIon : public Particle {
public:
GCRIon();
G4double pdf(G4double E) { return C1 * exp(-C2 * E) * (1 - exp(-C3 * E + C4)) / sum; }
};
#endif

View File

@ -1,24 +1,27 @@
#ifndef DESCSS_PrimaryGeneratorAction_h
#define DESCSS_PrimaryGeneratorAction_h
#include "G4ParticleGun.hh"
#include "G4GeneralParticleSource.hh"
#include "G4ParticleGun.hh"
#include "G4VUserPrimaryGeneratorAction.hh"
#include "globals.hh"
class G4GeneralParticleSource;
class G4ParticleGun;
class G4Event;
class G4Box;
template <typename T>
class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction {
public:
PrimaryGeneratorAction();
~PrimaryGeneratorAction();
virtual void GeneratePrimaries(G4Event*);
const G4GeneralParticleSource* GetParticleGun() const { return fParticleGun; }
const T* GetParticleGun() const { return fParticleGun; }
private:
G4GeneralParticleSource* fParticleGun;
G4String particleType;
T* fParticleGun = nullptr;
};
#endif

View File

@ -9,9 +9,15 @@
#include "G4VisExecutive.hh"
#include "QBBC.hh"
extern G4String particleType;
int main(int argc, char** argv) {
G4UIExecutive* ui = nullptr;
if (argc == 1) ui = new G4UIExecutive(argc, argv);
if (argc >= 3)
particleType = argv[2];
else
particleType = "TE";
G4MTRunManager* runManager = new G4MTRunManager;
G4VisExecutive* visManager = new G4VisExecutive;

View File

@ -1,8 +1,27 @@
#include "ActionInitialization.h"
#include "PrimaryGeneratorAction.h"
#include <time.h>
ActionInitialization::ActionInitialization() {}
#include "Randomize.hh"
G4String particleType;
class G4GeneralParticleSource;
class G4ParticleGun;
ActionInitialization::ActionInitialization() {
G4long seed = time(NULL);
CLHEP::HepRandom::setTheEngine(new CLHEP::RanecuEngine());
CLHEP::HepRandom::setTheSeed(seed);
}
ActionInitialization::~ActionInitialization() {}
void ActionInitialization::Build() const { SetUserAction(new PrimaryGeneratorAction); }
void ActionInitialization::BuildForMaster() const {}
void ActionInitialization::Build() const {
if (particleType == "TE")
SetUserAction(new PrimaryGeneratorAction<G4GeneralParticleSource>);
else
SetUserAction(new PrimaryGeneratorAction<G4ParticleGun>);
}

View File

@ -17,8 +17,6 @@ DetectorConstruction::DetectorConstruction() {}
DetectorConstruction::~DetectorConstruction() {}
std::ofstream dataFile("F:/Project/Geant4/DESCSS/debug.txt");
static void ConstructSectionSphere(G4LogicalVolume* fMotherLogical, G4double zBias) {
G4double pRadius = 2.8 * m / 2;
G4double pRadius2 = 1.2 * m;
@ -255,7 +253,5 @@ G4VPhysicalVolume* DetectorConstruction::Construct() {
// 大柱段
ConstructSectionBig(logicWorld, 0.815 * m);
dataFile << logicWorld->GetMass() / kg << std::endl;
return physicsWorld;
}

View File

@ -1,9 +1,86 @@
#include "PrimaryGeneratorAction.h"
PrimaryGeneratorAction::PrimaryGeneratorAction() {
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "Randomize.hh"
PrimaryGeneratorAction<G4GeneralParticleSource>::PrimaryGeneratorAction() {
fParticleGun = new G4GeneralParticleSource();
}
PrimaryGeneratorAction::~PrimaryGeneratorAction() { delete fParticleGun; }
PrimaryGeneratorAction<G4ParticleGun>::PrimaryGeneratorAction() {
fParticleGun = new G4ParticleGun();
}
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* e) { fParticleGun->GeneratePrimaryVertex(e); }
template <typename T>
PrimaryGeneratorAction<T>::~PrimaryGeneratorAction() {
delete fParticleGun;
}
std::string line;
G4double M = 931.5;
G4double C1, C2, C3, C4, sum;
std::ifstream modelFile("F:/Project/Geant4/DESCSS/assets/model.txt");
G4double randomPhi() {
G4double v = G4UniformRand();
return 2 * pi * v;
}
G4double randomTheta() {
G4double u = G4UniformRand();
return std::acos(2 * u - 1);
}
G4ThreeVector randomPos(G4double rho) {
G4double phi = randomPhi();
G4double theta = randomTheta();
G4double x = rho * std::sin(theta) * std::cos(phi);
G4double y = rho * std::sin(theta) * std::sin(phi);
G4double z = rho * std::cos(theta);
return G4ThreeVector(x, y, z);
}
G4ThreeVector randomDir() {
G4double phi = randomPhi();
G4double theta = randomTheta();
G4double x = -std::sin(theta) * std::cos(phi);
G4double y = -std::sin(theta) * std::sin(phi);
G4double z = -std::cos(theta);
return G4ThreeVector(x, y, z);
}
G4double ITM(G4double Emin, G4double Emax, G4double (*f)(G4double)) {
G4double x = G4UniformRand();
return (*f)(x);
}
G4double ARM(G4double Emin, G4double Emax, G4double (*f)(G4double)) {
G4double x, y;
while (true) {
x = G4UniformRand() * (Emax - Emin) + Emin;
y = G4UniformRand();
if (y <= (*f)(x)) break;
}
return x;
}
void PrimaryGeneratorAction<G4GeneralParticleSource>::GeneratePrimaries(G4Event* e) {
fParticleGun->GeneratePrimaryVertex(e);
}
void PrimaryGeneratorAction<G4ParticleGun>::GeneratePrimaries(G4Event* e) {
G4double E = 0;
G4double R = 15. * m;
G4ThreeVector pos = G4ThreeVector(); // randomPos(R);
G4ThreeVector dir = G4ThreeVector(); // randomDir();
fParticleGun->SetParticleEnergy(E);
fParticleGun->SetParticlePosition(pos);
fParticleGun->SetParticleMomentumDirection(dir);
fParticleGun->GeneratePrimaryVertex(e);
}

83
utils/check.py Normal file
View File

@ -0,0 +1,83 @@
import math
import numpy as np
from matplotlib import pyplot as plt
def trap_proton():
def pdf(E, C1=559.76377, C2=-1.53424, C3=2.40873, sum=1.51565):
return C1 * pow(E - C2, -C3) / sum
pdf = np.vectorize(pdf)
def pdf_int(E, C1=559.76377, C2=-1.53424, C3=2.40873, sum=1.5156509):
return C1 * math.pow(E - C2, 1 - C3) / (1 - C3) / sum
def pdf_sum(E):
return pdf_int(E) - pdf_int(50)
def pdf_inv(y, C1=559.76377, C2=-1.53424, C3=2.40873, sum=1.5156509):
return pow(sum * (pdf_int(50) + y) * (1 - C3) / C1, 1 / (1 - C3)) + C2
with open('debug.txt', 'r') as f:
data = f.readlines()
data = np.array([float(s) for s in data[0].split(', ')[:-1]])
x = np.linspace(50, 1000, 1000)
plt.plot(x, pdf(x), label="f(x)")
plt.hist(data, bins=np.linspace(50, 1000, 500), density=True, histtype='step', label="Simulation")
plt.legend()
plt.ylabel('Frequency')
plt.xlabel('Energy/(MeV)')
plt.savefig('docs/simu-tp.png', bbox_inches='tight')
def gcr_proton():
M = 931.5
def get_pc(E):
return math.sqrt(2 * M * E + E * E)
def get_beta(E):
pc = get_pc(E)
return pc / math.sqrt(M * M + pc * pc)
def pdf(E, C1=125.89223, alpha=44.00673, C2=-0.000460574, C3=112356, sum=361.781665965026):
pc = get_pc(E)
beta = get_beta(E)
return C1 * math.pow(beta, alpha - 1) * math.pow(pc / (pc + C2), C3) / pc / sum
pdf = np.vectorize(pdf)
with open('debug.txt', 'r') as f:
data = f.readlines()
data = np.array([float(s) for s in data[0].split(', ')[:-1]])
x = np.linspace(220, 1e5, 10000)
plt.plot(x, pdf(x), label="f(x)")
plt.hist(data, bins=np.linspace(220, 1e5, 100), density=True, histtype='step', label="Simulation")
plt.legend()
plt.ylabel('Frequency')
plt.xlabel('Energy/(MeV)')
plt.savefig('docs/simu-gcrp.png', bbox_inches='tight')
def gcr_ion():
def pdf(E, C1=0.01815, C2=0.000297726, C3=0.00184, C4=0.49271, sum=48.4143048):
return C1 * math.exp(-C2 * E) * (1 - math.exp(-C3 * E + C4)) / sum
pdf = np.vectorize(pdf)
with open('debug.txt', 'r') as f:
data = f.readlines()
data = np.array([float(s) for s in data[0].split(', ')[:-1]])
x = np.linspace(220, 1e5, 10000)
plt.plot(x, pdf(x), label="f(x)")
plt.hist(data, bins=np.linspace(220, 1e5, 300), density=True, histtype='step', label="Simulation")
plt.legend()
plt.ylabel('Frequency')
plt.xlabel('Energy/(MeV)')
plt.savefig('docs/simu-gcri.png', bbox_inches='tight')
trap_proton()
gcr_proton()
gcr_ion()