diff --git a/CMakeLists.txt b/CMakeLists.txt index dc14bbb..f6cce1f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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}) diff --git a/G4.code-workspace b/G4.code-workspace index 220a01e..38cfd23 100644 --- a/G4.code-workspace +++ b/G4.code-workspace @@ -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" } } } \ No newline at end of file diff --git a/README.md b/README.md index e9ee44b..91d5006 100644 --- a/README.md +++ b/README.md @@ -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 @@ $$

-2. 质子 +2. 质子 与电子保持一致，幂律模型修改为$j=j_0(E-\beta)^{-\alpha}$，进行`BIC test`，最终选用幂律模型$j_0=559.76377,\ \alpha=2.40873,\ \beta=-1.53424$。
@@ -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 |
@@ -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`本质上是一种模拟算法，效率很低，但是适应性更广，特别是对于一些复杂分布函数。 -
- -| $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}$ | - -
- -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$已经归一化） +
+ + | $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}$ | + +
+
+* 俘获质子的函数较为复杂，使用`ARM`生成`100000`个随机数需要`12`分钟，效率满足要求； +
+* 其余俘获辐射粒子同样使用`ARM`，生成`100000`个随机数需要`4.5`分钟。 +
+ ## 空间站建模 1. 尺寸与分区[^6]
diff --git a/assets/model.txt b/assets/model.txt index 245715a..bc48252 100644 --- a/assets/model.txt +++ b/assets/model.txt @@ -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 diff --git a/debug.txt b/debug.txt index 8824d44..e69de29 100644 --- a/debug.txt +++ b/debug.txt @@ -1 +0,0 @@ -22210 diff --git a/docs/simu-gcri.png b/docs/simu-gcri.png new file mode 100644 index 0000000..60b7802 Binary files /dev/null and b/docs/simu-gcri.png differ diff --git a/docs/simu-gcrp.png b/docs/simu-gcrp.png new file mode 100644 index 0000000..1293183 Binary files /dev/null and b/docs/simu-gcrp.png differ diff --git a/docs/simu-tp.png b/docs/simu-tp.png new file mode 100644 index 0000000..65068d9 Binary files /dev/null and b/docs/simu-tp.png differ diff --git a/include/ActionInitialization.h b/include/ActionInitialization.h index 7f3f7c7..791b0cd 100644 --- a/include/ActionInitialization.h +++ b/include/ActionInitialization.h @@ -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 diff --git a/include/Particle.h b/include/Particle.h new file mode 100644 index 0000000..77ffb41 --- /dev/null +++ b/include/Particle.h @@ -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 \ No newline at end of file diff --git a/include/PrimaryGeneratorAction.h b/include/PrimaryGeneratorAction.h index a6f37e1..6680458 100644 --- a/include/PrimaryGeneratorAction.h +++ b/include/PrimaryGeneratorAction.h @@ -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 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 diff --git a/main.cpp b/main.cpp index 262f28f..e51f9cb 100644 --- a/main.cpp +++ b/main.cpp @@ -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; diff --git a/src/ActionInitialization.cpp b/src/ActionInitialization.cpp index f563642..ae36d1f 100644 --- a/src/ActionInitialization.cpp +++ b/src/ActionInitialization.cpp @@ -1,8 +1,27 @@ #include "ActionInitialization.h" #include "PrimaryGeneratorAction.h" +#include -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); + else + SetUserAction(new PrimaryGeneratorAction); +} diff --git a/src/DetectorConstruction.cpp b/src/DetectorConstruction.cpp index deedaa8..84cd5af 100644 --- a/src/DetectorConstruction.cpp +++ b/src/DetectorConstruction.cpp @@ -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; } diff --git a/src/PrimaryGeneratorAction.cpp b/src/PrimaryGeneratorAction.cpp index 1bf72f9..c083a63 100644 --- a/src/PrimaryGeneratorAction.cpp +++ b/src/PrimaryGeneratorAction.cpp @@ -1,9 +1,86 @@ #include "PrimaryGeneratorAction.h" -PrimaryGeneratorAction::PrimaryGeneratorAction() { +#include "G4PhysicalConstants.hh" +#include "G4SystemOfUnits.hh" +#include "Randomize.hh" + +PrimaryGeneratorAction::PrimaryGeneratorAction() { fParticleGun = new G4GeneralParticleSource(); } -PrimaryGeneratorAction::~PrimaryGeneratorAction() { delete fParticleGun; } +PrimaryGeneratorAction::PrimaryGeneratorAction() { + fParticleGun = new G4ParticleGun(); +} -void PrimaryGeneratorAction::GeneratePrimaries(G4Event* e) { fParticleGun->GeneratePrimaryVertex(e); } +template +PrimaryGeneratorAction::~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::GeneratePrimaries(G4Event* e) { + fParticleGun->GeneratePrimaryVertex(e); +} + +void PrimaryGeneratorAction::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); +} diff --git a/utils/check.py b/utils/check.py new file mode 100644 index 0000000..386783f --- /dev/null +++ b/utils/check.py @@ -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()