Merge branch 'dev'

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
 )

G4.code-workspace

@@ -119,7 +119,8 @@
 			"style_parser": "cpp",
 			"vmanip": "cpp",
 			"element": "cpp",
-			"tree": "cpp"
+			"tree": "cpp",
+			"vertices": "cpp"
 		}
 	}
 }

README.md

@@ -296,7 +296,7 @@ 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$已经归一化）

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

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

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<std::string, G4bool> sensitivities;
     G4HumanPhantomMaterial* material;
     G4HumanPhantomMessenger* messenger;

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

include/Matrix.h

@@ -41,10 +41,11 @@ matrix<T>::matrix(int n, int m) {
 
 template <class T>
 matrix<T>::matrix(const matrix<T>& A) {
-    size_t number = A.n * A.m;
+    n = A.n, m = 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);
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) element[i][j] = A.element[i][j];
 }
 
 template <class T>
@@ -54,10 +55,11 @@ template <class T>
 matrix<T>& matrix<T>::operator=(const matrix<T>& A) {
     if (this != &A) {
         delete[] element;
-        size_t number = A.n * A.m;
+        n = A.n, m = 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);
+        for (int i = 0; i < n; i++)
+            for (int j = 0; j < m; j++) element[i][j] = A.element[i][j];
     }
     return *this;
 }

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;
 };
 

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;

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);
+    }
 }

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<const DetectorConstruction*>(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 x = rho * sin(theta) * cos(phi);
-    G4double y = rho * std::sin(theta) * std::sin(phi);
+    G4double y = rho * sin(theta) * sin(phi);
-    G4double z = rho * std::cos(theta);
+    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 u = sin(theta) * cos(phi);
-    G4double v = -std::sin(theta) * std::sin(phi);
+    G4double v = sin(theta) * sin(phi);
-    G4double w = -std::cos(theta);
+    G4double w = cos(theta);
     G4double mD[3][1] = {{u}, {v}, {w}};
     matrix<G4double> 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<G4double> 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<G4double> M(3, 1);
-    M = R * D, M = M + O;
+    M = R * D;
-
-    M.Output();
 
     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);