diff --git a/CMakeLists.txt b/CMakeLists.txt index f6cce1f..4ce8e36 100644 --- 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 ) diff --git a/G4.code-workspace b/G4.code-workspace index 3c470bf..79e8674 100644 --- 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" } } } \ No newline at end of file diff --git a/README.md b/README.md index 9041dd6..36c889d 100644 --- 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`分钟。
+## 体模 +使用`Geant4`示例中的`MIRD`人体模型,分为男、女两个性别, +
+ ## 空间站建模 -1. 尺寸与分区[^6] +1. 尺寸与分区[^7]
* 全长:16.6 m @@ -261,8 +347,8 @@ $$ 2. 材料 * 外壳: 3层 2mm 铝 + 10mm 芳纶 + 5mm 铝
- * 5系铝合金[^7] - * 泰普龙[^8] + * 5系铝合金[^8] + * 泰普龙[^9]
* 填充:金属为主,使得总重大致相当 * 内部:空气 @@ -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) diff --git a/debug.mac b/debug.mac deleted file mode 100644 index 6ade0ff..0000000 --- a/debug.mac +++ /dev/null @@ -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 diff --git a/docs/MIRD.png b/docs/MIRD.png new file mode 100644 index 0000000..fce6bf7 Binary files /dev/null and b/docs/MIRD.png differ diff --git a/electron.mac b/electron.mac deleted file mode 100644 index 2bd3993..0000000 --- a/electron.mac +++ /dev/null @@ -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 diff --git a/include/G4BasePhantomBuilder.h b/include/G4BasePhantomBuilder.h index d9a9e41..dd6e8f5 100644 --- 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" diff --git a/include/G4FemaleBuilder.h b/include/G4FemaleBuilder.h index d03d5bc..a7fa691 100644 --- 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" diff --git a/include/G4HumanPhantomColour.h b/include/G4HumanPhantomColour.h index f1b4ff0..fef4efa 100644 --- 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" diff --git a/include/G4HumanPhantomHit.h b/include/G4HumanPhantomHit.h index 1efa78e..6720115 100644 --- 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" diff --git a/include/G4HumanPhantomMaterial.h b/include/G4HumanPhantomMaterial.h index 2676361..e837b01 100644 --- 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; diff --git a/include/G4HumanPhantomMessenger.h b/include/G4HumanPhantomMessenger.h index fc24452..1547472 100644 --- 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; diff --git a/include/G4HumanPhantomSD.h b/include/G4HumanPhantomSD.h index 1d895eb..a78ad52 100644 --- 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" diff --git a/include/G4MIRDBodyFactory.h b/include/G4MIRDBodyFactory.h index 0b98542..24e38f8 100644 --- 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" diff --git a/include/G4MIRDBrain.h b/include/G4MIRDBrain.h index c9b570d..bdc65a8 100644 --- 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" diff --git a/include/G4MIRDHead.h b/include/G4MIRDHead.h index 68dfd28..5e58e9d 100644 --- 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" diff --git a/include/G4MIRDHeart.h b/include/G4MIRDHeart.h index d35e93f..f28ea49 100644 --- 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" diff --git a/include/G4MIRDLeftAdrenal.h b/include/G4MIRDLeftAdrenal.h index 551969c..014a80c 100644 --- 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" diff --git a/include/G4MIRDLeftArmBone.h b/include/G4MIRDLeftArmBone.h index 6813bd2..517f7bc 100644 --- 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" diff --git a/include/G4MIRDLeftBreast.h b/include/G4MIRDLeftBreast.h index 66ce674..5f19dbe 100644 --- 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" diff --git a/include/G4MIRDLeftClavicle.h b/include/G4MIRDLeftClavicle.h index c00211f..c05bd8c 100644 --- 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" diff --git a/include/G4MIRDLeftKidney.h b/include/G4MIRDLeftKidney.h index 46fec50..fd7d202 100644 --- 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" diff --git a/include/G4MIRDLeftLeg.h b/include/G4MIRDLeftLeg.h index ac4172f..c164181 100644 --- 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" diff --git a/include/G4MIRDLeftLegBone.h b/include/G4MIRDLeftLegBone.h index 5fb4516..6372733 100644 --- 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" diff --git a/include/G4MIRDLeftLung.h b/include/G4MIRDLeftLung.h index f44fb49..e869aff 100644 --- 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" diff --git a/include/G4MIRDLeftOvary.h b/include/G4MIRDLeftOvary.h index bab1c5e..e889619 100644 --- 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" diff --git a/include/G4MIRDLeftScapula.h b/include/G4MIRDLeftScapula.h index a81afa6..c0dd1cd 100644 --- 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" diff --git a/include/G4MIRDLeftTeste.h b/include/G4MIRDLeftTeste.h index a7ce69c..50404e0 100644 --- 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" diff --git a/include/G4MIRDLiver.h b/include/G4MIRDLiver.h index 1b56d39..6633e05 100644 --- 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" diff --git a/include/G4MIRDLowerLargeIntestine.h b/include/G4MIRDLowerLargeIntestine.h index 90424b5..aef2508 100644 --- 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" diff --git a/include/G4MIRDMaleGenitalia.h b/include/G4MIRDMaleGenitalia.h index e25c358..2c93267 100644 --- 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" diff --git a/include/G4MIRDMiddleLowerSpine.h b/include/G4MIRDMiddleLowerSpine.h index 0e0b5b7..5814e8e 100644 --- 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" diff --git a/include/G4MIRDPancreas.h b/include/G4MIRDPancreas.h index 758a486..696bb84 100644 --- 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" diff --git a/include/G4MIRDPelvis.h b/include/G4MIRDPelvis.h index 0571fcd..ce59818 100644 --- 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" diff --git a/include/G4MIRDRibCage.h b/include/G4MIRDRibCage.h index 670846a..113cd1c 100644 --- 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" diff --git a/include/G4MIRDRightAdrenal.h b/include/G4MIRDRightAdrenal.h index 2480de9..bb19031 100644 --- 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" diff --git a/include/G4MIRDRightArmBone.h b/include/G4MIRDRightArmBone.h index 466849f..2da88fd 100644 --- 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" diff --git a/include/G4MIRDRightBreast.h b/include/G4MIRDRightBreast.h index 6953853..1fb9818 100644 --- 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" diff --git a/include/G4MIRDRightClavicle.h b/include/G4MIRDRightClavicle.h index 228eef4..a3e6499 100644 --- 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" diff --git a/include/G4MIRDRightKidney.h b/include/G4MIRDRightKidney.h index fbed525..2825573 100644 --- 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" diff --git a/include/G4MIRDRightLeg.h b/include/G4MIRDRightLeg.h index bb6280f..f25accd 100644 --- 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" diff --git a/include/G4MIRDRightLegBone.h b/include/G4MIRDRightLegBone.h index 315de24..f9752bb 100644 --- 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" diff --git a/include/G4MIRDRightLung.h b/include/G4MIRDRightLung.h index 4918734..aa3e932 100644 --- 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" diff --git a/include/G4MIRDRightOvary.h b/include/G4MIRDRightOvary.h index 71b39d4..90371a9 100644 --- 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" diff --git a/include/G4MIRDRightScapula.h b/include/G4MIRDRightScapula.h index 6dbc6bd..3d081db 100644 --- 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" diff --git a/include/G4MIRDRightTeste.h b/include/G4MIRDRightTeste.h index 0345efd..f2e9f7d 100644 --- 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" diff --git a/include/G4MIRDSkull.h b/include/G4MIRDSkull.h index 91ec2c5..bd31673 100644 --- 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" diff --git a/include/G4MIRDSmallIntestine.h b/include/G4MIRDSmallIntestine.h index e15650e..24c4ea5 100644 --- 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" diff --git a/include/G4MIRDSpleen.h b/include/G4MIRDSpleen.h index 3698adf..25b3e38 100644 --- 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" diff --git a/include/G4MIRDStomach.h b/include/G4MIRDStomach.h index 43fd9de..1c70d75 100644 --- 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" diff --git a/include/G4MIRDThymus.h b/include/G4MIRDThymus.h index 4692347..ce54541 100644 --- 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" diff --git a/include/G4MIRDThyroid.h b/include/G4MIRDThyroid.h index 0651a89..2b24fe6 100644 --- 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" diff --git a/include/G4MIRDTrunk.h b/include/G4MIRDTrunk.h index c07abf0..e2e825b 100644 --- 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" diff --git a/include/G4MIRDUpperLargeIntestine.h b/include/G4MIRDUpperLargeIntestine.h index b21019c..3e05b93 100644 --- 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" diff --git a/include/G4MIRDUpperSpine.h b/include/G4MIRDUpperSpine.h index 26d5ca0..6752dd2 100644 --- 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" diff --git a/include/G4MIRDUrinaryBladder.h b/include/G4MIRDUrinaryBladder.h index dde98a9..1541722 100644 --- 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" diff --git a/include/G4MIRDUterus.h b/include/G4MIRDUterus.h index d31c845..4bf3b6c 100644 --- 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" diff --git a/include/G4MaleBuilder.h b/include/G4MaleBuilder.h index 96f3401..1f369f2 100644 --- 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" diff --git a/include/G4PhantomBuilder.h b/include/G4PhantomBuilder.h index 2fd7f20..ea4eeba 100644 --- 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" diff --git a/include/G4VBodyFactory.h b/include/G4VBodyFactory.h index 666658c..4cdb4c1 100644 --- 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; diff --git a/include/G4VOrgan.h b/include/G4VOrgan.h index 5c31dc3..398eab8 100644 --- 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" diff --git a/include/Matrix.h b/include/Matrix.h new file mode 100644 index 0000000..3888598 --- /dev/null +++ b/include/Matrix.h @@ -0,0 +1,131 @@ +#ifndef DESCSS_Matrix_h +#define DESCSS_Matrix_h + +template +class matrix { +public: + matrix(int n, int m); + matrix(const matrix&); + ~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& operator=(const matrix&); + matrix operator+() const; + matrix operator+(const matrix&) const; + matrix operator*(const matrix&) const; + matrix operator+(const T&) const; + matrix operator*(const T&) const; + matrix operator/(const T&) const; + +protected: + int n, m; + T** element; +}; + +template +matrix::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 +matrix::matrix(const matrix& 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 +matrix::~matrix() {} + +template +matrix& matrix::operator=(const matrix& 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 +void matrix::Input(T* p) { + for (int i = 0; i < n; i++) + for (int j = 0; j < m; j++) this->element[i][j] = *(p++); +} + +template +void matrix::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 +T& matrix::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 +matrix matrix::operator+(const matrix& A) const { + if (n != A.n || m != A.m) std::cout << "Matrix Size is Out of batch " << std::endl; + matrix 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 +matrix matrix::operator+(const T& C) const { + matrix 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 +matrix matrix::operator*(const matrix& A) const { + if (m != A.n) std::cout << "Matrix Style is Out of batch " << std::endl; + + matrix 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 +matrix matrix::operator*(const T& C) const { + matrix 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 +matrix matrix::operator/(const T& C) const { + matrix 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 \ No newline at end of file diff --git a/src/G4HumanPhantomHit.cpp b/src/G4HumanPhantomHit.cpp index a9ff19a..418fa22 100644 --- 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; } diff --git a/src/G4HumanPhantomSD.cpp b/src/G4HumanPhantomSD.cpp index 22bcdba..a331e6e 100644 --- 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 <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;iPrint(); + // 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(); } diff --git a/src/G4PhantomBuilder.cpp b/src/G4PhantomBuilder.cpp index 1e3814b..7ed160b 100644 --- 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(); } diff --git a/src/Material.cpp b/src/Material.cpp index e1c840d..5f585e7 100644 --- 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; } diff --git a/src/PrimaryGeneratorAction.cpp b/src/PrimaryGeneratorAction.cpp index 7bf8dab..bdbb803 100644 --- 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 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 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 R(3, 3); + R.Input(&mR[0][0]); + R = R / (1 + a * a + b * b + c * c); + + matrix 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);