diff --git a/README.md b/README.md index 91d5006..e289b0a 100644 --- a/README.md +++ b/README.md @@ -236,7 +236,7 @@ $$ 为了保证发射方向是朝向球面内，因此需要限制$\theta$的最大角度为$90^{\circ}$。 3. 能量 -* `GCR`质子的函数较为简单，采用`ITM`；（下表中$j_0$已经归一化） +* 俘获质子的函数较为简单，采用`ITM`；（下表中$j_0$已经归一化）
| $F(E)$ | $\int F(E)\mathrm{d}E$ | $C(x)=\int_a^xF(u)\mathrm{d}u$ | $C^{-1}(y)$ | @@ -245,9 +245,9 @@ $$
-* 俘获质子的函数较为复杂，使用`ARM`生成`100000`个随机数需要`12`分钟，效率满足要求； +* `GCR`质子的函数较为复杂，使用`ARM`生成`100000`个随机数需要`12`分钟，效率满足要求；
-* 其余俘获辐射粒子同样使用`ARM`，生成`100000`个随机数需要`4.5`分钟。 +* 其余`GCR`辐射粒子同样使用`ARM`，生成`100000`个随机数需要`4.5`分钟。
## 空间站建模 diff --git a/include/Particle.h b/include/Particle.h deleted file mode 100644 index 1752809..0000000 --- a/include/Particle.h +++ /dev/null @@ -1,73 +0,0 @@ -#ifndef DESCSS_Particle_h -#define DESCSS_Particle_h - -#include "G4IonTable.hh" -#include "G4ParticleDefinition.hh" -#include "G4ParticleTable.hh" -#include "globals.hh" - -class Particle { -public: - Particle(){}; - Particle(G4int Z, G4int A) { - if (Z == 1) - ion = G4ParticleTable::GetParticleTable()->FindParticle("proton"); - else if (Z == 2) - ion = G4ParticleTable::GetParticleTable()->FindParticle("alpha"); - else - ion = G4IonTable::GetIonTable()->GetIon(Z, A, 0); - this->Z = Z, this->A = A; - }; - virtual ~Particle(){}; - -public: - void setC1(G4double v) { C1 = v; }; - void setC2(G4double v) { C2 = v; }; - void setC3(G4double v) { C3 = v; }; - void setC4(G4double v) { C4 = v; }; - void setSum(G4double v) { sum = v; }; - virtual G4double pdf(G4double E) { return 0; }; - -public: - G4ParticleDefinition* ion; - -protected: - G4int Z, A; - G4double sum; - G4double C1, C2, C3, C4; -}; - -class TProton : public Particle { -public: - using Particle::Particle; - G4double pdf(G4double y) { - G4double C = C1 * pow(50 - C2, 1 - C3) / sum / (1 - C3); - return pow(sum * (C + y) * (1 - C3) / C1, 1 / (1 - C3)) + C2; - } -}; - -class GCRProton : public Particle { -public: - using Particle::Particle; - G4double pdf(G4double E) { - G4double pc = getPc(E); - G4double beta = getBeta(E); - return C1 * pow(beta, C2 - 1) * pow(pc / (pc + C3), C4) / pc / sum; - } - -private: - G4double M = 931.5; - G4double getPc(G4double E) { return sqrt(2 * M * E + E * E); } - G4double getBeta(G4double E) { - G4double pc = getPc(E); - return pc / sqrt(M * M + pc * pc); - } -}; - -class GCRIon : public Particle { -public: - using Particle::Particle; - G4double pdf(G4double E) { return C1 * exp(-C2 * E) * (1 - exp(-C3 * E + C4)) / sum; } -}; - -#endif \ No newline at end of file diff --git a/include/PrimaryGeneratorAction.h b/include/PrimaryGeneratorAction.h index 1ff5fa8..ebdd835 100644 --- a/include/PrimaryGeneratorAction.h +++ b/include/PrimaryGeneratorAction.h @@ -1,8 +1,6 @@ #ifndef DESCSS_PrimaryGeneratorAction_h #define DESCSS_PrimaryGeneratorAction_h -#include "Particle.h" - #include "G4GeneralParticleSource.hh" #include "G4ParticleGun.hh" #include "G4VUserPrimaryGeneratorAction.hh" diff --git a/src/PrimaryGeneratorAction.cpp b/src/PrimaryGeneratorAction.cpp index d76c726..6da5a8e 100644 --- a/src/PrimaryGeneratorAction.cpp +++ b/src/PrimaryGeneratorAction.cpp @@ -1,13 +1,15 @@ -#include "Particle.h" #include "PrimaryGeneratorAction.h" +#include "G4IonTable.hh" +#include "G4ParticleDefinition.hh" +#include "G4ParticleTable.hh" #include "G4PhysicalConstants.hh" #include "G4SystemOfUnits.hh" #include "Randomize.hh" -Particle fParticle; G4String particleType; 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"); @@ -26,25 +28,26 @@ 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 == "TP") { + if (tmp == "TP") std::getline(ss, tmp, ','), sum = std::stof(tmp); - fParticle = TProton(Z, A); - } else { + else { std::getline(ss, tmp, ','), C4 = std::stof(tmp); std::getline(ss, tmp, ','), sum = std::stof(tmp); - if (tmp == "GCR_H") - fParticle = GCRProton(Z, A); - else - fParticle = GCRIon(Z, A); - fParticle.setC4(C4); } - fParticle.setC1(C1), fParticle.setC2(C2), fParticle.setC3(C3), fParticle.setSum(sum); break; } } + G4ParticleDefinition* ion; fParticleGun = new G4ParticleGun(); - fParticleGun->SetParticleDefinition(fParticle.ion); + if (Z == 1) + ion = G4ParticleTable::GetParticleTable()->FindParticle("proton"); + else if (Z == 2) + ion = G4ParticleTable::GetParticleTable()->FindParticle("alpha"); + else + ion = G4IonTable::GetIonTable()->GetIon(Z, A, 0.); + fParticleGun->SetParticleCharge(Z * eplus); + fParticleGun->SetParticleDefinition(ion); } template @@ -99,7 +102,18 @@ G4double ARM(G4double Emin, G4double Emax, G4double (*f)(G4double)) { return x; } -G4double pdf(G4double E) { return fParticle.pdf(E); } +G4double pdfTProtonInv(G4double y) { + G4double C = C1 * pow(50 - C2, 1 - C3) / sum / (1 - C3); + return pow(sum * (C + y) * (1 - C3) / C1, 1 / (1 - C3)) + C2; +} + +G4double pdfGCRProton(G4double E) { + G4double pc = sqrt(2 * M * E + E * E); + G4double beta = pc / sqrt(M * M + pc * pc); + return C1 * pow(beta, C2 - 1) * pow(pc / (pc + C3), C4) / pc / sum; +} + +G4double pdfGCRIon(G4double E) { return C1 * exp(-C2 * E) * (1 - exp(-C3 * E + C4)) / sum; } void PrimaryGeneratorAction::GeneratePrimaries(G4Event* e) { fParticleGun->GeneratePrimaryVertex(e); @@ -112,9 +126,11 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* e) { G4ThreeVector dir = randomDir(); if (particleType == "TP") - E = ARM(50, 1000, pdf); + E = ITM(50, 1000, pdfTProtonInv); + else if (particleType == "GCR_H") + E = ARM(220, 1e5, pdfGCRProton); else - E = ARM(220, 1e5, pdf); + E = ARM(220, 1e5, pdfGCRIon); fParticleGun->SetParticleEnergy(E); fParticleGun->SetParticlePosition(pos);