fox
/
G4-DESCSS
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

fix: docs;change: particle pdf

This commit is contained in:
liuyihui 2022-05-17 01:08:42 +08:00
parent 49db68f847
commit 93bcbabded
4 changed files with 34 additions and 93 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
README.md
View File

@ -236,7 +236,7 @@ $$
为了保证发射方向是朝向球面内，因此需要限制$\theta$的最大角度为$90^{\circ}$。
3. 能量
* `GCR`质子的函数较为简单，采用`ITM`；（下表中$j_0$已经归一化）
* 俘获质子的函数较为简单，采用`ITM`；（下表中$j_0$已经归一化）
<div align=center>
| $F(E)$ | $\int F(E)\mathrm{d}E$ | $C(x)=\int_a^xF(u)\mathrm{d}u$ | $C^{-1}(y)$ |
@ -245,9 +245,9 @@ $$
</div>
<div align=center><img src="docs/simu-tp.png" style="max-width: 50%;"></div>
* 俘获质子的函数较为复杂，使用`ARM`生成`100000`个随机数需要`12`分钟，效率满足要求；
* `GCR`质子的函数较为复杂，使用`ARM`生成`100000`个随机数需要`12`分钟，效率满足要求；
<div align=center><img src="docs/simu-gcrp.png" style="max-width: 50%;"></div>
* 其余俘获辐射粒子同样使用`ARM`，生成`100000`个随机数需要`4.5`分钟。
* 其余`GCR`辐射粒子同样使用`ARM`，生成`100000`个随机数需要`4.5`分钟。
<div align=center><img src="docs/simu-gcri.png" style="max-width: 50%;"></div>
## 空间站建模

73
include/Particle.h
View File

@ -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

2
include/PrimaryGeneratorAction.h
View File

@ -1,8 +1,6 @@
#ifndef DESCSS_PrimaryGeneratorAction_h
#define DESCSS_PrimaryGeneratorAction_h
#include "Particle.h"
#include "G4GeneralParticleSource.hh"
#include "G4ParticleGun.hh"
#include "G4VUserPrimaryGeneratorAction.hh"

46
src/PrimaryGeneratorAction.cpp
View File

@ -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<G4ParticleGun>::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 <typename T>
@ -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<G4GeneralParticleSource>::GeneratePrimaries(G4Event* e) {
fParticleGun->GeneratePrimaryVertex(e);
@ -112,9 +126,11 @@ void PrimaryGeneratorAction<G4ParticleGun>::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);