diff --git a/CMakeLists.txt b/CMakeLists.txt index 4ce8e36..9538295 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -20,6 +20,7 @@ add_executable(DESCSS main.cpp ${sources} ${headers}) target_link_libraries(DESCSS ${Geant4_LIBRARIES}) set(DESCSS_SCRIPTS + electron.mac vis.mac assets/model.txt ) diff --git a/README.md b/README.md index 36c889d..3e740b4 100644 --- a/README.md +++ b/README.md @@ -222,11 +222,11 @@ $$ 随后我们需要将$\vec{d}$转化到`World`坐标系（简称为`W`系）中，已知两个坐标系有两个公共点，在`W`系中表示为原点$O$和发射点$\vec{r}=(x_0,y_0,z_0)$，在`P`系中则为$(0,0,\rho)$和原点$O^\prime$,则有： $$ \begin{bmatrix} -x\\ y\\ z +x\\\ y\\\ z \end{bmatrix} = \lambda R \begin{bmatrix} -u\\ v\\ w +u\\\ v\\\ w \end{bmatrix} + \begin{bmatrix} -x_0\\ y_0\\ z_0 +x_0\\\ y_0\\\ z_0 \end{bmatrix} \cdots(1) $$ @@ -234,8 +234,8 @@ $$ 其中$\lambda$为尺度比例因子，设置为1；$R$被称为罗德里格矩阵，定义反对称矩阵： $$ S = \begin{bmatrix} -0,& -c,& -b \\ -c,& 0,& -a \\ +0,& -c,& -b \\\ +c,& 0,& -a \\\ b,& a,& 0 \end{bmatrix} $$ @@ -243,8 +243,8 @@ $$ 则$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 \\ + 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) $$ @@ -252,23 +252,23 @@ $$ 将两个公共点代入$(1)$，相减则可消去平移项，可得到： $$ \begin{bmatrix} -x_2-x_1\\ y_2-y_1\\ z_2-z_1 +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 +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} \\ +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 +a\\\ b\\\ c \end{bmatrix} = \begin{bmatrix} -u_{21} - x_{21} \\ -v_{21} - y_{21} \\ +u_{21} - x_{21} \\\ +v_{21} - y_{21} \\\ w_{21} - z_{21} \end{bmatrix} $$ @@ -276,14 +276,14 @@ $$ 代入本题，则有： $$ \begin{bmatrix} -0,& -\rho+z_0,& y_0 \\ --\rho+z_0,& 0,& -x_0 \\ +0,& -\rho+z_0,& y_0 \\\ +-\rho+z_0,& 0,& -x_0 \\\ -y_0,& -x_0,& 0 \end{bmatrix}\begin{bmatrix} -a\\ b\\ c +a\\\ b\\\ c \end{bmatrix} = \begin{bmatrix} --x_0 \\ --y_0 \\ +-x_0 \\\ +-y_0 \\\ -\rho-z_0 \end{bmatrix} $$ @@ -291,12 +291,12 @@ $$ 令$a=1$，则可求得： $$ \begin{cases} -b = \frac{z_0+\rho-y_0}{x_0} \\ +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`系中的表达式。 +计算得到$R$，继而根据$(1)$计算得到$\vec{d}$在`W`系中的表达式（最终无需平移，即不需要$+\begin{bmatrix}x_0\\\ y_0\\\ z_0\end{bmatrix}$）。 3. 能量 * 俘获电子的函数较为简单，采用`ITM`；（下表中$j_0$已经归一化） diff --git a/assets/model.txt b/assets/model.txt index d27adb9..76f5b5e 100644 --- a/assets/model.txt +++ b/assets/model.txt @@ -1,4 +1,4 @@ -TE, 1, -1, 1333680.43009, 0.0824, 32652.96179377461 +TE, -1, 1, 1333680.43009, 0.0824, 32652.96179377461 TP, 1, 1, 559.76377, -1.53424, 2.40873, 1.5156509 GCR_H, 1, 1, 125.89223, 44.00673, -0.000460574, 112356, 361.781666 GCR_He, 2, 4, 0.01815, 0.000297726, 0.00184, 0.49271, 48.4143048 diff --git a/electron.mac b/electron.mac new file mode 100644 index 0000000..1ca1a78 --- /dev/null +++ b/electron.mac @@ -0,0 +1,15 @@ +# 多线程设置 +/run/numberOfThreads 1 + +# verbose +/control/verbose 2 +/run/verbose 0 +/event/verbose 0 +/tracking/verbose 0 +/vis/verbose errors + +# 初始化 +/run/initialize + +/type/set TE +/run/beamOn 50 \ No newline at end of file diff --git a/include/DetectorConstruction.h b/include/DetectorConstruction.h index e71d068..3aede27 100644 --- a/include/DetectorConstruction.h +++ b/include/DetectorConstruction.h @@ -20,7 +20,9 @@ public: void SetBodyPartSensitivity(G4String, G4bool); void SetPhantomSex(G4String); void SetPhantomModel(G4String); + void SetParticleType(G4String s) { particleType = s; }; void ConstructSDandField(); + G4String GetParticleType() const { return particleType; }; private: void ConstructSectionSphere(G4LogicalVolume* fMotherLogical, G4double zBias); @@ -33,6 +35,7 @@ private: private: G4String sex = "Male"; G4String model = "MIRD"; + G4String particleType = "TE"; std::map sensitivities; G4HumanPhantomMaterial* material; G4HumanPhantomMessenger* messenger; diff --git a/include/G4HumanPhantomMessenger.h b/include/G4HumanPhantomMessenger.h index 1547472..5b2ce99 100644 --- a/include/G4HumanPhantomMessenger.h +++ b/include/G4HumanPhantomMessenger.h @@ -2,6 +2,7 @@ #define DESCSS_HumanPhantomMessenger_h 1 class DetectorConstruction; +class PrimaryGeneratorAction; class G4UIcommand; class G4UIdirectory; @@ -24,10 +25,12 @@ private: DetectorConstruction* myUserPhantom; G4UIdirectory* phantomDir; + G4UIdirectory* particleTypeDir; G4UIcmdWithAString* modelCmd; G4UIcmdWithAString* sexCmd; G4UIcmdWithoutParameter* endCmd; + G4UIcmdWithAString* typeCmd; }; #endif diff --git a/include/PrimaryGeneratorAction.h b/include/PrimaryGeneratorAction.h index 3abed07..b0f48a3 100644 --- a/include/PrimaryGeneratorAction.h +++ b/include/PrimaryGeneratorAction.h @@ -1,6 +1,7 @@ #ifndef DESCSS_PrimaryGeneratorAction_h #define DESCSS_PrimaryGeneratorAction_h +#include "G4ParticleDefinition.hh" #include "G4ParticleGun.hh" #include "G4VUserPrimaryGeneratorAction.hh" #include "globals.hh" @@ -8,17 +9,18 @@ class G4ParticleGun; class G4Event; -extern G4String particleType; - class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction { public: PrimaryGeneratorAction(); - ~PrimaryGeneratorAction(); + + G4ParticleDefinition* DefineParticle(); virtual void GeneratePrimaries(G4Event*); const G4ParticleGun* GetParticleGun() const { return fParticleGun; } private: + G4int Z; + G4String particleType = "TE"; G4ParticleGun* fParticleGun = nullptr; }; diff --git a/main.cpp b/main.cpp index e51f9cb..262f28f 100644 --- a/main.cpp +++ b/main.cpp @@ -9,15 +9,9 @@ #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/G4HumanPhantomMessenger.cpp b/src/G4HumanPhantomMessenger.cpp index bd9658b..0079034 100644 --- a/src/G4HumanPhantomMessenger.cpp +++ b/src/G4HumanPhantomMessenger.cpp @@ -7,10 +7,13 @@ #include "G4UIdirectory.hh" #include "globals.hh" -G4HumanPhantomMessenger::G4HumanPhantomMessenger(DetectorConstruction* myUsrPhtm) : myUserPhantom(myUsrPhtm) { +G4HumanPhantomMessenger::G4HumanPhantomMessenger(DetectorConstruction* myUsrPhtm) : myUserPhantom(myUsrPhtm){ phantomDir = new G4UIdirectory("/phantom/"); phantomDir->SetGuidance("Set Your Phantom."); + phantomDir = new G4UIdirectory("/type/"); + phantomDir->SetGuidance("Set Your Particle Type"); + modelCmd = new G4UIcmdWithAString("/phantom/setPhantomModel", this); modelCmd->SetGuidance("Set sex of Phantom: MIRD, MIRDHead."); modelCmd->SetParameterName("phantomModel", true); @@ -28,6 +31,15 @@ G4HumanPhantomMessenger::G4HumanPhantomMessenger(DetectorConstruction* myUsrPhtm endCmd = new G4UIcmdWithoutParameter("/phantom/buildNewPhantom", this); endCmd->SetGuidance("Build your Phantom."); endCmd->AvailableForStates(G4State_PreInit, G4State_Idle); + + typeCmd = new G4UIcmdWithAString("/type/set", this); + typeCmd->SetGuidance("Set type of Particle: TE, TP, GCR_H, GCR_He etc"); + typeCmd->SetParameterName("particleType", true); + typeCmd->SetDefaultValue("TE"); + typeCmd->SetCandidates( + "TE TP GCR_H GCR_He GCR_Li GCR_Be GCR_B GCR_C GCR_N GCR_O GCR_F GCR_Ne GCR_Na GCR_Mg GCR_Al GCR_Si GCR_P GCR_S " + "GCR_Cl GCR_Ar GCR_K GCR_Ca GCR_Sc GCR_Ti GCR_V GCR_Cr GCR_Mn GCR_Fe GCR_Co GCR_Ni"); + typeCmd->AvailableForStates(G4State_PreInit, G4State_Idle); } G4HumanPhantomMessenger::~G4HumanPhantomMessenger() { @@ -47,4 +59,7 @@ void G4HumanPhantomMessenger::SetNewValue(G4UIcommand* command, G4String newValu if (command == endCmd) { G4cout << " ****************>>>> NEW PHANTOM CONSTRUCTION <<<<***************** " << G4endl; } + if (command == typeCmd) { + myUserPhantom->SetParticleType(newValue); + } } diff --git a/src/PrimaryGeneratorAction.cpp b/src/PrimaryGeneratorAction.cpp index bdbb803..f032be7 100644 --- a/src/PrimaryGeneratorAction.cpp +++ b/src/PrimaryGeneratorAction.cpp @@ -1,3 +1,4 @@ +#include "DetectorConstruction.h" #include "Matrix.h" #include "PrimaryGeneratorAction.h" @@ -5,17 +6,30 @@ #include "G4ParticleDefinition.hh" #include "G4ParticleTable.hh" #include "G4PhysicalConstants.hh" +#include "G4RunManager.hh" #include "G4SystemOfUnits.hh" #include "Randomize.hh" -G4String particleType; -G4double Z, A; G4double M = 931.5; G4double C1, C2, C3, C4, sum; -std::string tmp, line, name; std::ifstream modelFile("assets/model.txt"); PrimaryGeneratorAction::PrimaryGeneratorAction() { + G4ParticleDefinition* ion = DefineParticle(); + fParticleGun->SetParticleCharge(this->Z * eplus); + fParticleGun->SetParticleDefinition(ion); + + const DetectorConstruction* detConstruction = + static_cast(G4RunManager::GetRunManager()->GetUserDetectorConstruction()); + particleType = detConstruction->GetParticleType(); +} + +PrimaryGeneratorAction::~PrimaryGeneratorAction() { delete fParticleGun; } + +G4ParticleDefinition* PrimaryGeneratorAction::DefineParticle() { + G4double Z, A; + std::string tmp, line, name; + while (std::getline(modelFile, line)) { std::stringstream ss(line); std::getline(ss, name, ','); @@ -47,11 +61,10 @@ PrimaryGeneratorAction::PrimaryGeneratorAction() { ion = G4ParticleTable::GetParticleTable()->FindParticle("alpha"); else ion = G4IonTable::GetIonTable()->GetIon(Z, A, 0.); - fParticleGun->SetParticleCharge(Z * eplus); - fParticleGun->SetParticleDefinition(ion); -} -PrimaryGeneratorAction::~PrimaryGeneratorAction() { delete fParticleGun; } + this->Z = Z; + return ion; +} G4double randomPhi() { G4double v = G4UniformRand(); @@ -61,35 +74,35 @@ G4double randomPhi() { G4double randomTheta(G4double minU = -1, G4double maxU = 1) { G4double u = G4UniformRand(); u = u * (maxU - minU) + minU; - return std::acos(u); + return acos(u); } -G4ThreeVector randomPos(G4double rho = 15. * m) { +G4ThreeVector randomPos(G4double rho = 12. * m) { 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); + G4double x = rho * sin(theta) * cos(phi); + G4double y = rho * sin(theta) * sin(phi); + G4double z = rho * cos(theta); return G4ThreeVector(x, y, z); } -G4ThreeVector randomDir(G4ThreeVector pos, G4double maxTheta = 90. * deg) { +G4ThreeVector randomDir(G4ThreeVector pos, G4double maxTheta = 30. * deg) { G4double phi = randomPhi(); G4double theta = randomTheta(cos(maxTheta), 1); - G4double u = -std::sin(theta) * std::cos(phi); - G4double v = -std::sin(theta) * std::sin(phi); - G4double w = -std::cos(theta); + G4double u = sin(theta) * cos(phi); + G4double v = sin(theta) * sin(phi); + G4double w = cos(theta); G4double mD[3][1] = {{u}, {v}, {w}}; matrix D(3, 1); D.Input(&mD[0][0]); - G4double rho = pos.getRho(); G4double x0 = pos.getX(); G4double y0 = pos.getY(); G4double z0 = pos.getZ(); + G4double rho = sqrt(x0 * x0 + y0 * y0 + z0 * z0); G4double mO[3][1] = {{x0}, {y0}, {z0}}; matrix O(3, 1); O.Input(&mO[0][0]); @@ -106,9 +119,7 @@ G4ThreeVector randomDir(G4ThreeVector pos, G4double maxTheta = 90. * deg) { R = R / (1 + a * a + b * b + c * c); matrix M(3, 1); - M = R * D, M = M + O; - - M.Output(); + M = R * D; return G4ThreeVector(M(1, 1), M(2, 1), M(3, 1)); } @@ -160,7 +171,7 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* e) { else E = ARM(220, 1e5, pdfGCRIon); - fParticleGun->SetParticleEnergy(E); + fParticleGun->SetParticleEnergy(E * MeV); fParticleGun->SetParticlePosition(pos); fParticleGun->SetParticleMomentumDirection(dir); fParticleGun->GeneratePrimaryVertex(e);