try to fix bug
This commit is contained in:
parent
87c330f37a
commit
a8d6d82050
|
|
@ -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
|
||||
)
|
||||
|
|
|
|||
42
README.md
42
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$已经归一化）
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
};
|
||||
|
||||
|
|
|
|||
6
main.cpp
6
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;
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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 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<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.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);
|
||||
|
|
|
|||
Loading…
Reference in New Issue