fix: dir vector

2022-05-18 20:01:10 +08:00 · 2022-05-18 20:01:10 +08:00 · 87c330f37a
parent 8f4475d3e1
commit 87c330f37a
67 changed files with 505 additions and 209 deletions
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@ -20,8 +20,6 @@ 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
 )
--- a/G4.code-workspace
+++ b/G4.code-workspace
@ -10,7 +10,116 @@
 			"line": "cpp",
 			"plane": "cpp",
 			"qrot": "cpp",
-			"cmath": "cpp"
+			"cmath": "cpp",
+			"string": "cpp",
+			"vector": "cpp",
+			"algorithm": "cpp",
+			"array": "cpp",
+			"atomic": "cpp",
+			"bit": "cpp",
+			"cctype": "cpp",
+			"cfenv": "cpp",
+			"charconv": "cpp",
+			"chrono": "cpp",
+			"clocale": "cpp",
+			"compare": "cpp",
+			"complex": "cpp",
+			"concepts": "cpp",
+			"condition_variable": "cpp",
+			"csetjmp": "cpp",
+			"csignal": "cpp",
+			"cstdarg": "cpp",
+			"cstddef": "cpp",
+			"cstdint": "cpp",
+			"cstdio": "cpp",
+			"cstdlib": "cpp",
+			"cstring": "cpp",
+			"ctime": "cpp",
+			"cwchar": "cpp",
+			"deque": "cpp",
+			"exception": "cpp",
+			"format": "cpp",
+			"forward_list": "cpp",
+			"fstream": "cpp",
+			"functional": "cpp",
+			"future": "cpp",
+			"initializer_list": "cpp",
+			"iomanip": "cpp",
+			"ios": "cpp",
+			"iosfwd": "cpp",
+			"iostream": "cpp",
+			"istream": "cpp",
+			"iterator": "cpp",
+			"limits": "cpp",
+			"list": "cpp",
+			"locale": "cpp",
+			"map": "cpp",
+			"memory": "cpp",
+			"mutex": "cpp",
+			"new": "cpp",
+			"numeric": "cpp",
+			"optional": "cpp",
+			"ostream": "cpp",
+			"queue": "cpp",
+			"ratio": "cpp",
+			"regex": "cpp",
+			"set": "cpp",
+			"sstream": "cpp",
+			"stack": "cpp",
+			"stdexcept": "cpp",
+			"stop_token": "cpp",
+			"streambuf": "cpp",
+			"system_error": "cpp",
+			"thread": "cpp",
+			"tuple": "cpp",
+			"type_traits": "cpp",
+			"typeinfo": "cpp",
+			"unordered_map": "cpp",
+			"utility": "cpp",
+			"xfacet": "cpp",
+			"xhash": "cpp",
+			"xiosbase": "cpp",
+			"xlocale": "cpp",
+			"xlocbuf": "cpp",
+			"xlocinfo": "cpp",
+			"xlocmes": "cpp",
+			"xlocmon": "cpp",
+			"xlocnum": "cpp",
+			"xloctime": "cpp",
+			"xmemory": "cpp",
+			"xstddef": "cpp",
+			"xstring": "cpp",
+			"xtr1common": "cpp",
+			"xtree": "cpp",
+			"xutility": "cpp",
+			"matrix": "cpp",
+			"vec2f": "cpp",
+			"vec3f": "cpp",
+			"vec4f": "cpp",
+			"sf_vec": "cpp",
+			"sf_vec3f": "cpp",
+			"aida_ntuple": "cpp",
+			"handle": "cpp",
+			"img": "cpp",
+			"mat": "cpp",
+			"mapmanip": "cpp",
+			"basket": "cpp",
+			"branch": "cpp",
+			"dummy": "cpp",
+			"graph": "cpp",
+			"info": "cpp",
+			"iros": "cpp",
+			"leaf": "cpp",
+			"named": "cpp",
+			"obj_array": "cpp",
+			"obj_list": "cpp",
+			"stl_vector": "cpp",
+			"tree_index": "cpp",
+			"vector3": "cpp",
+			"style_parser": "cpp",
+			"vmanip": "cpp",
+			"element": "cpp",
+			"tree": "cpp"
 		}
 	}
 }
--- a/README.md
+++ b/README.md
@ -207,14 +207,96 @@ $$
    double z = rho * cos(theta);
    ```

-2. 发射方向  
-在`Geant4`中，假定粒子的位置矢量为$\vec{r}$，则以$\vec{r}$反方向作为$z$轴正方向建立右手坐标系，粒子的指向为:
+2. 发射方向[^6]  
+在`Geant4`中，假定粒子的位置矢量为$\vec{r}$，则以$\vec{r}$反方向作为$z$轴正方向建立右手坐标系（简称为`P`系），随机生成指向$\vec{d}$，为了保证发射方向是朝向球面内，因此需要限制$\theta$的最大角度为$90^{\circ}$。
+    ```c
+    double u = rand();
+    double v = rand();
+    double phi = 2 * pi * u;
+    double theta = arccos(2 * v - 1);
+    double x = sin(theta) * cos(phi);
+    double y = sin(theta) * sin(phi);
+    double z = cos(theta);
+    ```
+
+随后我们需要将$\vec{d}$转化到`World`坐标系（简称为`W`系）中，已知两个坐标系有两个公共点，在`W`系中表示为原点$O$和发射点$\vec{r}=(x_0,y_0,z_0)$，在`P`系中则为$(0,0,\rho)$和原点$O^\prime$,则有：
 $$
-P_x = -\sin\theta\cos\phi \\
-P_y = -\sin\theta\sin\phi \\
-P_z = -\cos\theta
+\begin{bmatrix}
+x\\ y\\ z
+\end{bmatrix} = \lambda R \begin{bmatrix}
+u\\ v\\ w
+\end{bmatrix} + \begin{bmatrix}
+x_0\\ y_0\\ z_0
+\end{bmatrix}
+\cdots(1)
 $$
-为了保证发射方向是朝向球面内，因此需要限制$\theta$的最大角度为$90^{\circ}$。
+
+其中$\lambda$为尺度比例因子，设置为1；$R$被称为罗德里格矩阵，定义反对称矩阵：
+$$
+S = \begin{bmatrix}
+0,& -c,& -b \\
+c,& 0,& -a \\
+b,& a,& 0
+\end{bmatrix}
+$$
+
+则$R$可表示为$R = (I+S)(I-S)^{-1}$，展开后则有：
+$$
+R = \frac{1}{1+a^2+b^2+c^2}\begin{bmatrix}
+    1 + a^2 - b^2 - c^2, -2c - 2ab, -2b + 2ac \\
+    2c - 2ab, 1 - a^2 + b^2 - c^2, -2a - 2bc \\
+    2b + 2ac, 2a - 2bc, 1 - a^2 - b^2 + c^2
+\end{bmatrix}\cdots(2)
+$$
+
+将两个公共点代入$(1)$，相减则可消去平移项，可得到：
+$$
+\begin{bmatrix}
+x_2-x_1\\ y_2-y_1\\ z_2-z_1
+\end{bmatrix} = R \begin{bmatrix}
+u_2-u_1\\ v_2-v_1\\ w_2-w_1
+\end{bmatrix}
+$$
+
+代入$(2)$则可得到（$u_{21}=u_2-u_1$）：
+$$
+\begin{bmatrix}
+0,& w_{21}+z_{21},& v_{21}+y_{21} \\
+w_{21}+z_{21},& 0,& -u_{21}-x_{21} \\
+-v_{21}-y_{21},& -u_{21}-x_{21},& 0
+\end{bmatrix}\begin{bmatrix}
+a\\ b\\ c
+\end{bmatrix} = \begin{bmatrix}
+u_{21} - x_{21} \\
+v_{21} - y_{21} \\
+w_{21} - z_{21}
+\end{bmatrix}
+$$
+
+代入本题，则有：
+$$
+\begin{bmatrix}
+0,& -\rho+z_0,& y_0 \\
+-\rho+z_0,& 0,& -x_0 \\
+-y_0,& -x_0,& 0
+\end{bmatrix}\begin{bmatrix}
+a\\ b\\ c
+\end{bmatrix} = \begin{bmatrix}
+-x_0 \\
+-y_0 \\
+-\rho-z_0
+\end{bmatrix}
+$$
+
+令$a=1$，则可求得：
+$$
+\begin{cases}
+b = \frac{z_0+\rho-y_0}{x_0} \\
+c = \frac{z_0-\rho+y_0}{x_0}
+\end{cases}
+$$
+
+计算得到$R$，继而根据$(1)$计算得到$\vec{d}$在`W`系中的表达式。

 3. 能量  
 * 俘获电子的函数较为简单，采用`ITM`；（下表中$j_0$已经归一化）
@ -240,8 +322,12 @@ $$
 * 其余`GCR`辐射粒子同样使用`ARM`，生成`100000`个随机数需要`4.5`分钟。
    <div align=center><img src="docs/simu-gcri.png" style="max-width: 50%;"></div>

+## 体模
+使用`Geant4`示例中的`MIRD`人体模型，分为男、女两个性别，
+<div align=center><img src="docs/MIRD.png" style="max-width: 60%;"></div>
+
 ## 空间站建模
-1. 尺寸与分区[^6]  
+1. 尺寸与分区[^7]  
 <div align=center><img src="docs/size.webp" style="max-width: 50%;"></div>

 * 全长：16.6 m
@ -261,8 +347,8 @@ $$
 2. 材料
 * 外壳： 3层 2mm 铝 + 10mm 芳纶 + 5mm 铝
        <div align=center><img src="docs/Shell.png" style="max-width: 40%;"></div>
-    * 5系铝合金[^7]
-    * 泰普龙[^8]
+    * 5系铝合金[^8]
+    * 泰普龙[^9]
        <div align=center><img src="docs/Taparan.png" style="max-width: 30%;"></div>
 * 填充：金属为主，使得总重大致相当
 * 内部：空气
@ -289,6 +375,7 @@ $$
 [^3]: Bourdarie S, Xapsos M. The near-earth space radiation environment[J]. IEEE transactions on nuclear science, 2008, 55(4): 1810-1832.
 [^4]: Matthiä D, Berger T, Mrigakshi A I, et al. A ready-to-use galactic cosmic ray model[J]. Advances in Space Research, 2013, 51(3): 329-338.
 [^5]: Zhao H, Johnston W R, Baker D N, et al. Characterization and evolution of radiation belt electron energy spectra based on the Van Allen Probes measurements[J]. Journal of Geophysical Research: Space Physics, 2019, 124(6): 4217-4232.
-[^6]: [【知识点·航天】“天和”核心舱、“天宫”空间站和新一代载人飞船的最新知识点（干货版）](https://zhuanlan.zhihu.com/p/103709953)
-[^7]: [中铝造为“天和”号核心舱披“铠甲”壮“筋骨”](https://m.thepaper.cn/baijiahao_12484370)
-[^8]: [泰普龙产品相关知识](https://wenku.baidu.com/view/a2ecf93501d8ce2f0066f5335a8102d276a261cd.html)
+[^6]: [三维坐标系之间转换方法](https://zhuanlan.zhihu.com/p/458827599)
+[^7]: [【知识点·航天】“天和”核心舱、“天宫”空间站和新一代载人飞船的最新知识点（干货版）](https://zhuanlan.zhihu.com/p/103709953)
+[^8]: [中铝造为“天和”号核心舱披“铠甲”壮“筋骨”](https://m.thepaper.cn/baijiahao_12484370)
+[^9]: [泰普龙产品相关知识](https://wenku.baidu.com/view/a2ecf93501d8ce2f0066f5335a8102d276a261cd.html)
--- a/debug.mac
+++ b/debug.mac
@ -1,28 +0,0 @@
-# 多线程设置
-/run/numberOfThreads 1
-
-# verbose
-/control/verbose 1
-/run/verbose 0
-/event/verbose 0
-/tracking/verbose 0
-/gps/verbose 0
-
-# 初始化
-/run/initialize
-
-# 电子射程
-/gps/particle e-
-/gps/energy 100 keV
-/gps/position 0 0 5 m
-/gps/ang/type iso
-/gps/ang/maxphi 60 deg
-/gps/ang/maxtheta 20 deg
-
-# 质子射程
-/gps/particle proton
-/gps/energy 50 MeV
-/gps/position 0 0 5 m
-/gps/ang/type iso
-/gps/ang/maxphi 60 deg
-/gps/ang/maxtheta 20 deg
--- a/docs/MIRD.png
+++ b/docs/MIRD.png
--- a/electron.mac
+++ b/electron.mac
@ -1,26 +0,0 @@
-# 多线程设置
-/run/numberOfThreads 12
-
-# verbose
-/control/verbose 1
-/run/verbose 0
-/event/verbose 0
-/tracking/verbose 0
-
-# 初始化
-/run/initialize
-
-# Trapped electron
-/gps/particle e-
-/gps/pos/type Surface
-/gps/pos/shape Sphere
-/gps/pos/radius 15 m
-/gps/ene/type Exp
-/gps/ene/min 100 keV
-/gps/ene/max 10 MeV
-/gps/ene/ezero 0.0824
-/gps/ang/type iso
-/gps/ang/maxtheta 90 deg
-/gps/verbose 0
-
-/run/beamOn 923242
--- a/include/G4BasePhantomBuilder.h
+++ b/include/G4BasePhantomBuilder.h
@ -1,5 +1,5 @@
-#ifndef G4BasePhantomBuilder_h
-#define G4BasePhantomBuilder_h 1
+#ifndef DESCSS_BasePhantomBuilder_h
+#define DESCSS_BasePhantomBuilder_h 1

 #include "G4VPhysicalVolume.hh"

--- a/include/G4FemaleBuilder.h
+++ b/include/G4FemaleBuilder.h
@ -1,5 +1,5 @@
-#ifndef G4FemaleBuilder_h
-#define G4FemaleBuilder_h 1
+#ifndef DESCSS_FemaleBuilder_h
+#define DESCSS_FemaleBuilder_h 1

 #include "G4PhantomBuilder.h"

--- a/include/G4HumanPhantomColour.h
+++ b/include/G4HumanPhantomColour.h
@ -1,5 +1,5 @@
-#ifndef G4HumanPhantomColour_H
-#define G4HumanPhantomColour_H 1
+#ifndef DESCSS_HumanPhantomColour_H
+#define DESCSS_HumanPhantomColour_H 1

 #include "G4Colour.hh"
 #include "globals.hh"
--- a/include/G4HumanPhantomHit.h
+++ b/include/G4HumanPhantomHit.h
@ -1,5 +1,5 @@
-#ifndef G4HumanPhantomHit_h
-#define G4HumanPhantomHit_h 1
+#ifndef DESCSS_HumanPhantomHit_h
+#define DESCSS_HumanPhantomHit_h 1

 #include "G4Allocator.hh"
 #include "G4THitsCollection.hh"
--- a/include/G4HumanPhantomMaterial.h
+++ b/include/G4HumanPhantomMaterial.h
@ -1,5 +1,5 @@
-#ifndef G4HumanPhantomMaterial_H
-#define G4HumanPhantomMaterial_H 1
+#ifndef DESCSS_HumanPhantomMaterial_H
+#define DESCSS_HumanPhantomMaterial_H 1

 #include "globals.hh"
 class G4Material;
--- a/include/G4HumanPhantomMessenger.h
+++ b/include/G4HumanPhantomMessenger.h
@ -1,5 +1,5 @@
-#ifndef G4HumanPhantomMessenger_h
-#define G4HumanPhantomMessenger_h 1
+#ifndef DESCSS_HumanPhantomMessenger_h
+#define DESCSS_HumanPhantomMessenger_h 1

 class DetectorConstruction;

--- a/include/G4HumanPhantomSD.h
+++ b/include/G4HumanPhantomSD.h
@ -1,5 +1,5 @@
-#ifndef G4HumanPhantomSD_h
-#define G4HumanPhantomSD_h 1
+#ifndef DESCSS_HumanPhantomSD_h
+#define DESCSS_HumanPhantomSD_h 1

 #include "G4HumanPhantomHit.h"

--- a/include/G4MIRDBodyFactory.h
+++ b/include/G4MIRDBodyFactory.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDBodyFactory_h
-#define G4MIRDBodyFactory_h 1
+#ifndef DESCSS_MIRDBodyFactory_h
+#define DESCSS_MIRDBodyFactory_h 1

 #include "G4VBodyFactory.h"
 #include "G4VOrgan.h"
--- a/include/G4MIRDBrain.h
+++ b/include/G4MIRDBrain.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDBrain_h
-#define G4MIRDBrain_h 1
+#ifndef DESCSS_MIRDBrain_h
+#define DESCSS_MIRDBrain_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDHead.h
+++ b/include/G4MIRDHead.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDHead_h
-#define G4MIRDHead_h 1
+#ifndef DESCSS_MIRDHead_h
+#define DESCSS_MIRDHead_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDHeart.h
+++ b/include/G4MIRDHeart.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDHeart_h
-#define G4MIRDHeart_h 1
+#ifndef DESCSS_MIRDHeart_h
+#define DESCSS_MIRDHeart_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLeftAdrenal.h
+++ b/include/G4MIRDLeftAdrenal.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLeftAdrenal_h
-#define G4MIRDLeftAdrenal_h 1
+#ifndef DESCSS_MIRDLeftAdrenal_h
+#define DESCSS_MIRDLeftAdrenal_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLeftArmBone.h
+++ b/include/G4MIRDLeftArmBone.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLeftArmBone_h
-#define G4MIRDLeftArmBone_h 1
+#ifndef DESCSS_MIRDLeftArmBone_h
+#define DESCSS_MIRDLeftArmBone_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLeftBreast.h
+++ b/include/G4MIRDLeftBreast.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLeftBreast_h
-#define G4MIRDLeftBreast_h 1
+#ifndef DESCSS_MIRDLeftBreast_h
+#define DESCSS_MIRDLeftBreast_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLeftClavicle.h
+++ b/include/G4MIRDLeftClavicle.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLeftClavicle_h
-#define G4MIRDLeftClavicle_h 1
+#ifndef DESCSS_MIRDLeftClavicle_h
+#define DESCSS_MIRDLeftClavicle_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLeftKidney.h
+++ b/include/G4MIRDLeftKidney.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLeftKidney_h
-#define G4MIRDLeftKidney_h 1
+#ifndef DESCSS_MIRDLeftKidney_h
+#define DESCSS_MIRDLeftKidney_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLeftLeg.h
+++ b/include/G4MIRDLeftLeg.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLeftLeg_h
-#define G4MIRDLeftLeg_h 1
+#ifndef DESCSS_MIRDLeftLeg_h
+#define DESCSS_MIRDLeftLeg_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLeftLegBone.h
+++ b/include/G4MIRDLeftLegBone.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLeftLegBone_h
-#define G4MIRDLeftLegBone_h 1
+#ifndef DESCSS_MIRDLeftLegBone_h
+#define DESCSS_MIRDLeftLegBone_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLeftLung.h
+++ b/include/G4MIRDLeftLung.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLeftLung_h
-#define G4MIRDLeftLung_h 1
+#ifndef DESCSS_MIRDLeftLung_h
+#define DESCSS_MIRDLeftLung_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLeftOvary.h
+++ b/include/G4MIRDLeftOvary.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLeftOvary_h
-#define G4MIRDLeftOvary_h 1
+#ifndef DESCSS_MIRDLeftOvary_h
+#define DESCSS_MIRDLeftOvary_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLeftScapula.h
+++ b/include/G4MIRDLeftScapula.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLeftScapula_h
-#define G4MIRDLeftScapula_h 1
+#ifndef DESCSS_MIRDLeftScapula_h
+#define DESCSS_MIRDLeftScapula_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLeftTeste.h
+++ b/include/G4MIRDLeftTeste.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLeftTeste_h
-#define G4MIRDLeftTeste_h 1
+#ifndef DESCSS_MIRDLeftTeste_h
+#define DESCSS_MIRDLeftTeste_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLiver.h
+++ b/include/G4MIRDLiver.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLiver_h
-#define G4MIRDLiver_h 1
+#ifndef DESCSS_MIRDLiver_h
+#define DESCSS_MIRDLiver_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDLowerLargeIntestine.h
+++ b/include/G4MIRDLowerLargeIntestine.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDLowerLargeIntestine_h
-#define G4MIRDLowerLargeIntestine_h 1
+#ifndef DESCSS_MIRDLowerLargeIntestine_h
+#define DESCSS_MIRDLowerLargeIntestine_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDMaleGenitalia.h
+++ b/include/G4MIRDMaleGenitalia.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDMaleGenitalia_h
-#define G4MIRDMaleGenitalia_h 1
+#ifndef DESCSS_MIRDMaleGenitalia_h
+#define DESCSS_MIRDMaleGenitalia_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDMiddleLowerSpine.h
+++ b/include/G4MIRDMiddleLowerSpine.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDMiddleLowerSpine_h
-#define G4MIRDMiddleLowerSpine_h 1
+#ifndef DESCSS_MIRDMiddleLowerSpine_h
+#define DESCSS_MIRDMiddleLowerSpine_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDPancreas.h
+++ b/include/G4MIRDPancreas.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDPancreas_h
-#define G4MIRDPancreas_h 1
+#ifndef DESCSS_MIRDPancreas_h
+#define DESCSS_MIRDPancreas_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDPelvis.h
+++ b/include/G4MIRDPelvis.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDPelvis_h
-#define G4MIRDPelvis_h 1
+#ifndef DESCSS_MIRDPelvis_h
+#define DESCSS_MIRDPelvis_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDRibCage.h
+++ b/include/G4MIRDRibCage.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDRibCage_h
-#define G4MIRDRibCage_h 1
+#ifndef DESCSS_MIRDRibCage_h
+#define DESCSS_MIRDRibCage_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDRightAdrenal.h
+++ b/include/G4MIRDRightAdrenal.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDRightAdrenal_h
-#define G4MIRDRightAdrenal_h 1
+#ifndef DESCSS_MIRDRightAdrenal_h
+#define DESCSS_MIRDRightAdrenal_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDRightArmBone.h
+++ b/include/G4MIRDRightArmBone.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDRightArmBone_h
-#define G4MIRDRightArmBone_h 1
+#ifndef DESCSS_MIRDRightArmBone_h
+#define DESCSS_MIRDRightArmBone_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDRightBreast.h
+++ b/include/G4MIRDRightBreast.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDRightBreast_h
-#define G4MIRDRightBreast_h 1
+#ifndef DESCSS_MIRDRightBreast_h
+#define DESCSS_MIRDRightBreast_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDRightClavicle.h
+++ b/include/G4MIRDRightClavicle.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDRightClavicle_h
-#define G4MIRDRightClavicle_h 1
+#ifndef DESCSS_MIRDRightClavicle_h
+#define DESCSS_MIRDRightClavicle_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDRightKidney.h
+++ b/include/G4MIRDRightKidney.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDRightKidney_h
-#define G4MIRDRightKidney_h 1
+#ifndef DESCSS_MIRDRightKidney_h
+#define DESCSS_MIRDRightKidney_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDRightLeg.h
+++ b/include/G4MIRDRightLeg.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDRightLeg_h
-#define G4MIRDRightLeg_h 1
+#ifndef DESCSS_MIRDRightLeg_h
+#define DESCSS_MIRDRightLeg_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDRightLegBone.h
+++ b/include/G4MIRDRightLegBone.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDRightLegBone_h
-#define G4MIRDRightLegBone_h 1
+#ifndef DESCSS_MIRDRightLegBone_h
+#define DESCSS_MIRDRightLegBone_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDRightLung.h
+++ b/include/G4MIRDRightLung.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDRightLung_h
-#define G4MIRDRightLung_h 1
+#ifndef DESCSS_MIRDRightLung_h
+#define DESCSS_MIRDRightLung_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDRightOvary.h
+++ b/include/G4MIRDRightOvary.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDRightOvary_h
-#define G4MIRDRightOvary_h 1
+#ifndef DESCSS_MIRDRightOvary_h
+#define DESCSS_MIRDRightOvary_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDRightScapula.h
+++ b/include/G4MIRDRightScapula.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDRightScapula_h
-#define G4MIRDRightScapula_h 1
+#ifndef DESCSS_MIRDRightScapula_h
+#define DESCSS_MIRDRightScapula_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDRightTeste.h
+++ b/include/G4MIRDRightTeste.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDRightTeste_h
-#define G4MIRDRightTeste_h 1
+#ifndef DESCSS_MIRDRightTeste_h
+#define DESCSS_MIRDRightTeste_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDSkull.h
+++ b/include/G4MIRDSkull.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDSkull_h
-#define G4MIRDSkull_h 1
+#ifndef DESCSS_MIRDSkull_h
+#define DESCSS_MIRDSkull_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDSmallIntestine.h
+++ b/include/G4MIRDSmallIntestine.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDSmallIntestine_h
-#define G4MIRDSmallIntestine_h 1
+#ifndef DESCSS_MIRDSmallIntestine_h
+#define DESCSS_MIRDSmallIntestine_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDSpleen.h
+++ b/include/G4MIRDSpleen.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDSpleen_h
-#define G4MIRDSpleen_h 1
+#ifndef DESCSS_MIRDSpleen_h
+#define DESCSS_MIRDSpleen_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDStomach.h
+++ b/include/G4MIRDStomach.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDStomach_h
-#define G4MIRDStomach_h 1
+#ifndef DESCSS_MIRDStomach_h
+#define DESCSS_MIRDStomach_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDThymus.h
+++ b/include/G4MIRDThymus.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDThymus_h
-#define G4MIRDThymus_h 1
+#ifndef DESCSS_MIRDThymus_h
+#define DESCSS_MIRDThymus_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDThyroid.h
+++ b/include/G4MIRDThyroid.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDThyroid_h
-#define G4MIRDThyroid_h 1
+#ifndef DESCSS_MIRDThyroid_h
+#define DESCSS_MIRDThyroid_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDTrunk.h
+++ b/include/G4MIRDTrunk.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDTrunk_h
-#define G4MIRDTrunk_h 1
+#ifndef DESCSS_MIRDTrunk_h
+#define DESCSS_MIRDTrunk_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDUpperLargeIntestine.h
+++ b/include/G4MIRDUpperLargeIntestine.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDUpperLargeIntestine_h
-#define G4MIRDUpperLargeIntestine_h 1
+#ifndef DESCSS_MIRDUpperLargeIntestine_h
+#define DESCSS_MIRDUpperLargeIntestine_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDUpperSpine.h
+++ b/include/G4MIRDUpperSpine.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDUpperSpine_h
-#define G4MIRDUpperSpine_h 1
+#ifndef DESCSS_MIRDUpperSpine_h
+#define DESCSS_MIRDUpperSpine_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDUrinaryBladder.h
+++ b/include/G4MIRDUrinaryBladder.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDUrinaryBladder_h
-#define G4MIRDUrinaryBladder_h 1
+#ifndef DESCSS_MIRDUrinaryBladder_h
+#define DESCSS_MIRDUrinaryBladder_h 1

 #include "G4VOrgan.h"

--- a/include/G4MIRDUterus.h
+++ b/include/G4MIRDUterus.h
@ -1,5 +1,5 @@
-#ifndef G4MIRDUterus_h
-#define G4MIRDUterus_h 1
+#ifndef DESCSS_MIRDUterus_h
+#define DESCSS_MIRDUterus_h 1

 #include "G4VOrgan.h"

--- a/include/G4MaleBuilder.h
+++ b/include/G4MaleBuilder.h
@ -1,5 +1,5 @@
-#ifndef G4MaleBuilder_h
-#define G4MaleBuilder_h 1
+#ifndef DESCSS_MaleBuilder_h
+#define DESCSS_MaleBuilder_h 1

 #include "G4PhantomBuilder.h"

--- a/include/G4PhantomBuilder.h
+++ b/include/G4PhantomBuilder.h
@ -1,5 +1,5 @@
-#ifndef G4PhantomBuilder_h
-#define G4PhantomBuilder_h 1
+#ifndef DESCSS_PhantomBuilder_h
+#define DESCSS_PhantomBuilder_h 1

 #include "G4BasePhantomBuilder.h"

--- a/include/G4VBodyFactory.h
+++ b/include/G4VBodyFactory.h
@ -1,5 +1,5 @@
-#ifndef G4VBodyFactory_h
-#define G4VBodyFactory_h 1
+#ifndef DESCSS_VBodyFactory_h
+#define DESCSS_VBodyFactory_h 1

 #include "G4VPhysicalVolume.hh"
 class G4VPhysicalVolume;
--- a/include/G4VOrgan.h
+++ b/include/G4VOrgan.h
@ -1,5 +1,5 @@
-#ifndef G4VOrgan_h
-#define G4VOrgan_h 1
+#ifndef DESCSS_VOrgan_h
+#define DESCSS_VOrgan_h 1

 #include "G4VPhysicalVolume.hh"

--- a/include/Matrix.h
+++ b/include/Matrix.h
@ -0,0 +1,131 @@
+#ifndef DESCSS_Matrix_h
+#define DESCSS_Matrix_h
+
+template <class T>
+class matrix {
+public:
+    matrix(int n, int m);
+    matrix(const matrix<T>&);
+    ~matrix();
+
+    void Input(T* p);
+    void Output();
+
+    int rows() const { return n; };
+    int cols() const { return m; };
+    T& operator()(int i, int j) const;
+
+    matrix<T>& operator=(const matrix<T>&);
+    matrix<T> operator+() const;
+    matrix<T> operator+(const matrix<T>&) const;
+    matrix<T> operator*(const matrix<T>&) const;
+    matrix<T> operator+(const T&) const;
+    matrix<T> operator*(const T&) const;
+    matrix<T> operator/(const T&) const;
+
+protected:
+    int n, m;
+    T** element;
+};
+
+template <class T>
+matrix<T>::matrix(int n, int m) {
+    if (m < 0 || n < 0) std::cout << "Rows and Cols must be >= 0 " << std::endl;
+    if ((m == 0 || n == 0) && (m != 0 || n != 0))
+        std::cout << "Either both or neither rows and columns should be zero " << std::endl;
+
+    this->m = m, this->n = n;
+    element = new T*[n];
+    for (int i = 0; i < n; i++) element[i] = new T[m];
+}
+
+template <class T>
+matrix<T>::matrix(const matrix<T>& A) {
+    size_t number = A.n * A.m;
+    element = new T*[n];
+    for (int i = 0; i < n; i++) element[i] = new T[m];
+    std::copy(A.element, A.element + number, element);
+}
+
+template <class T>
+matrix<T>::~matrix() {}
+
+template <class T>
+matrix<T>& matrix<T>::operator=(const matrix<T>& A) {
+    if (this != &A) {
+        delete[] element;
+        size_t number = A.n * A.m;
+        element = new T*[n];
+        for (int i = 0; i < n; i++) element[i] = new T[m];
+        std::copy(A.element, A.element + number, element);
+    }
+    return *this;
+}
+
+template <class T>
+void matrix<T>::Input(T* p) {
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) this->element[i][j] = *(p++);
+}
+
+template <class T>
+void matrix<T>::Output() {
+    for (int i = 0; i < n; i++) {
+        for (int j = 0; j < m; j++) std::cout << element[i][j] << " ";
+        std::cout << std::endl;
+    }
+}
+
+template <class T>
+T& matrix<T>::operator()(int i, int j) const {
+    if (i < 1 || i > n || j < 1 || j > m) std::cout << "Matrix Index Out Of Bounds " << std::endl;
+    return element[i - 1][j - 1];
+}
+
+template <class T>
+matrix<T> matrix<T>::operator+(const matrix<T>& A) const {
+    if (n != A.n || m != A.m) std::cout << "Matrix Size is Out of batch " << std::endl;
+    matrix<T> w(n, m);
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) w.element[i][j] = element[i][j] + A.element[i][j];
+    return w;
+}
+
+template <class T>
+matrix<T> matrix<T>::operator+(const T& C) const {
+    matrix<T> w(n, m);
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) w.element[i][j] = element[i][j] + C;
+    return w;
+}
+
+template <class T>
+matrix<T> matrix<T>::operator*(const matrix<T>& A) const {
+    if (m != A.n) std::cout << "Matrix Style is Out of batch " << std::endl;
+
+    matrix<T> w(n, m);
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < A.m; j++) {
+            w.element[i][j] = 0;
+            for (int k = 0; k < m; k++) w.element[i][j] += element[i][k] * A.element[k][j];
+        }
+    return w;
+}
+
+template <class T>
+matrix<T> matrix<T>::operator*(const T& C) const {
+    matrix<T> w(n, m);
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) w.element[i][j] = element[i][j] * C;
+    return w;
+}
+
+template <class T>
+matrix<T> matrix<T>::operator/(const T& C) const {
+    matrix<T> w(n, m);
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) w.element[i][j] = element[i][j] / C;
+    return w;
+}
+
+#endif
--- a/src/G4HumanPhantomHit.cpp
+++ b/src/G4HumanPhantomHit.cpp
@ -29,5 +29,5 @@ G4bool G4HumanPhantomHit::operator==(const G4HumanPhantomHit& right) const { ret
 void G4HumanPhantomHit::Draw() {}

 void G4HumanPhantomHit::Print() {
-    G4cout << "Energy deposit: " << G4BestUnit(edep, "Energy") << "BodyPartID: " << bodyPartID << G4endl;
+    G4cout << "Energy deposit: " << G4BestUnit(edep, "Energy") << ", BodyPartID: " << bodyPartID << G4endl;
 }
--- a/src/G4HumanPhantomSD.cpp
+++ b/src/G4HumanPhantomSD.cpp
@ -28,8 +28,6 @@ G4bool G4HumanPhantomSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {

    G4String bodypartName = aStep->GetPreStepPoint()->GetTouchable()->GetVolume()->GetLogicalVolume()->GetName();

-    // G4cout <<bodypartName <<":" << edep/MeV<< G4endl;
-
    G4HumanPhantomHit* newHit = new G4HumanPhantomHit();
    newHit->SetEdep(edep);
    newHit->SetBodyPartID(bodypartName);
@ -40,7 +38,7 @@ G4bool G4HumanPhantomSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {

 void G4HumanPhantomSD::EndOfEvent(G4HCofThisEvent*) {
    // G4int NbHits = collection->entries();
-    //      G4cout << "\n-------->Hits Collection: in this event they are " << NbHits
-    //             << " hits in the tracker chambers: " << G4endl;
-    //      for (G4int i=0;i<NbHits;i++) (*collection)[i]->Print();
+    // G4cout << "\n-------->Hits Collection: in this event, there are " << NbHits
+    //        << " hits in the tracker chambers: " << G4endl;
+    // for (G4int i = 0; i < NbHits; i++) (*collection)[i]->Print();
 }
--- a/src/G4PhantomBuilder.cpp
+++ b/src/G4PhantomBuilder.cpp
@ -334,6 +334,5 @@ void G4PhantomBuilder::SetMotherVolume(G4VPhysicalVolume* mother) { motherVolume

 void G4PhantomBuilder::SetModel(G4String modelFlag) {
    model = modelFlag;
-
-    if (model == "MIRD" || model == "MIX") body = new G4MIRDBodyFactory();
+    if (model == "MIRD") body = new G4MIRDBodyFactory();
 }
--- a/src/Material.cpp
+++ b/src/Material.cpp
@ -57,6 +57,4 @@ void DefineMaterials() {
    midMaterial->AddElementByMassFraction(nist->FindOrBuildElement("Zn"), .03);
    midMaterial->AddElementByMassFraction(nist->FindOrBuildElement("Cr"), .015);
    midMaterial->AddElementByMassFraction(nist->FindOrBuildElement("Ti"), .015);
-
-    std::cout << *(G4Material::GetMaterialTable()) << std::endl;
 }
--- a/src/PrimaryGeneratorAction.cpp
+++ b/src/PrimaryGeneratorAction.cpp
@ -1,3 +1,4 @@
+#include "Matrix.h"
 #include "PrimaryGeneratorAction.h"

 #include "G4IonTable.hh"
@ -12,7 +13,7 @@ G4double Z, A;
 G4double M = 931.5;
 G4double C1, C2, C3, C4, sum;
 std::string tmp, line, name;
-std::ifstream modelFile("F:/Project/Geant4/DESCSS/assets/model.txt");
+std::ifstream modelFile("assets/model.txt");

 PrimaryGeneratorAction::PrimaryGeneratorAction() {
    while (std::getline(modelFile, line)) {
@ -24,9 +25,9 @@ PrimaryGeneratorAction::PrimaryGeneratorAction() {
            std::getline(ss, tmp, ','), C1 = std::stof(tmp);
            std::getline(ss, tmp, ','), C2 = std::stof(tmp);
            std::getline(ss, tmp, ','), C3 = std::stof(tmp);
-            if (tmp == "TE")
+            if (name == "TE")
                sum = C3;
-            else if (tmp == "TP")
+            else if (name == "TP")
                std::getline(ss, tmp, ','), sum = std::stof(tmp);
            else {
                std::getline(ss, tmp, ','), C4 = std::stof(tmp);
@ -38,7 +39,9 @@ PrimaryGeneratorAction::PrimaryGeneratorAction() {

    G4ParticleDefinition* ion;
    fParticleGun = new G4ParticleGun();
-    if (Z == 1)
+    if (Z == -1)
+        ion = G4ParticleTable::GetParticleTable()->FindParticle("e-");
+    else if (Z == 1)
        ion = G4ParticleTable::GetParticleTable()->FindParticle("proton");
    else if (Z == 2)
        ion = G4ParticleTable::GetParticleTable()->FindParticle("alpha");
@ -55,12 +58,13 @@ G4double randomPhi() {
    return 2 * pi * v;
 }

-G4double randomTheta() {
+G4double randomTheta(G4double minU = -1, G4double maxU = 1) {
    G4double u = G4UniformRand();
-    return std::acos(2 * u - 1);
+    u = u * (maxU - minU) + minU;
+    return std::acos(u);
 }

-G4ThreeVector randomPos(G4double rho) {
+G4ThreeVector randomPos(G4double rho = 15. * m) {
    G4double phi = randomPhi();
    G4double theta = randomTheta();

@ -71,20 +75,47 @@ G4ThreeVector randomPos(G4double rho) {
    return G4ThreeVector(x, y, z);
 }

-G4ThreeVector randomDir() {
+G4ThreeVector randomDir(G4ThreeVector pos, G4double maxTheta = 90. * deg) {
    G4double phi = randomPhi();
-    G4double theta = randomTheta();
+    G4double theta = randomTheta(cos(maxTheta), 1);

-    G4double x = -std::sin(theta) * std::cos(phi);
-    G4double y = -std::sin(theta) * std::sin(phi);
-    G4double z = -std::cos(theta);
+    G4double u = -std::sin(theta) * std::cos(phi);
+    G4double v = -std::sin(theta) * std::sin(phi);
+    G4double w = -std::cos(theta);
+    G4double mD[3][1] = {{u}, {v}, {w}};
+    matrix<G4double> D(3, 1);
+    D.Input(&mD[0][0]);

-    return G4ThreeVector(x, y, z);
+    G4double rho = pos.getRho();
+    G4double x0 = pos.getX();
+    G4double y0 = pos.getY();
+    G4double z0 = pos.getZ();
+    G4double mO[3][1] = {{x0}, {y0}, {z0}};
+    matrix<G4double> O(3, 1);
+    O.Input(&mO[0][0]);
+
+    G4double a = 1;
+    G4double b = (z0 + rho - y0) / x0;
+    G4double c = (z0 - rho + y0) / x0;
+
+    G4double mR[3][3] = {{1 + a * a - b * b - c * c, -2 * c - 2 * a * b, -2 * b + 2 * a * c},
+                         {2 * c - 2 * a * b, 1 - a * a + b * b - c * c, -2 * a - 2 * b * c},
+                         {2 * b + 2 * a * c, 2 * a - 2 * b * c, 1 - a * a - b * b + c * c}};
+    matrix<G4double> R(3, 3);
+    R.Input(&mR[0][0]);
+    R = R / (1 + a * a + b * b + c * c);
+
+    matrix<G4double> M(3, 1);
+    M = R * D, M = M + O;
+
+    M.Output();
+
+    return G4ThreeVector(M(1, 1), M(2, 1), M(3, 1));
 }

 G4double ITM(G4double Emin, G4double Emax, G4double (*f)(G4double)) {
-    G4double x = G4UniformRand();
-    return (*f)(x);
+    G4double y = G4UniformRand();
+    return (*f)(y);
 }

 G4double ARM(G4double Emin, G4double Emax, G4double (*f)(G4double)) {
@ -98,8 +129,8 @@ G4double ARM(G4double Emin, G4double Emax, G4double (*f)(G4double)) {
 }

 G4double pdfTElectronInv(G4double y) {
-    G4double C = - C2 * C1 * exp(- 0.1 / C2) / sum;
-    return - C2 * log((C - y) * sum / C1 / C2);
+    G4double C = C2 * C1 * exp(-0.1 / C2) / sum;
+    return -C2 * log((C - y) * sum / C1 / C2);
 }

 G4double pdfTProtonInv(G4double y) {
@ -117,9 +148,8 @@ G4double pdfGCRIon(G4double E) { return C1 * exp(-C2 * E) * (1 - exp(-C3 * E + C

 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* e) {
    G4double E;
-    G4double R = 15. * m;
-    G4ThreeVector pos = randomPos(R);
-    G4ThreeVector dir = randomDir();
+    G4ThreeVector pos = randomPos();
+    G4ThreeVector dir = randomDir(pos);

    if (particleType == "TE")
        E = ITM(0.1, 10, pdfTElectronInv);