From 49db68f847c766198b0649e41e464136c8e4c5be Mon Sep 17 00:00:00 2001 From: YiHui Liu Date: Tue, 17 May 2022 00:34:06 +0800 Subject: [PATCH] add: Particle;change: model format --- G4.code-workspace | 3 +- assets/model.txt | 58 ++++++++++++++++---------------- include/ActionInitialization.h | 2 -- include/Particle.h | 22 +++++++----- include/PrimaryGeneratorAction.h | 5 ++- src/ActionInitialization.cpp | 3 +- src/PrimaryGeneratorAction.cpp | 53 ++++++++++++++++++++++++----- 7 files changed, 94 insertions(+), 52 deletions(-) diff --git a/G4.code-workspace b/G4.code-workspace index 38cfd23..ef0c72e 100644 --- a/G4.code-workspace +++ b/G4.code-workspace @@ -110,7 +110,8 @@ "cid": "cpp", "mem": "cpp", "sout": "cpp", - "typedefs": "cpp" + "typedefs": "cpp", + "senum": "cpp" } } } \ No newline at end of file diff --git a/assets/model.txt b/assets/model.txt index bc48252..a507c08 100644 --- a/assets/model.txt +++ b/assets/model.txt @@ -1,29 +1,29 @@ -559.76377, -1.53424, 2.40873, 1.5156509 -125.89223, 44.00673, -0.000460574, 112356, 361.781666 -0.01815, 0.000297726, 0.00184, 0.49271, 48.4143048 -0.000184794, 0.000389118, 0.0021, 0.49406, 0.3655716 -0.0000953693, 0.000386546, 0.002, 0.48222, 0.1883217 -0.000245835, 0.000383147, 0.00182, 0.45495, 0.4814399 -0.000551454, 0.000297607, 0.00184, 0.49310, 1.4715932 -0.000177397, 0.000347376, 0.00161, 0.43464, 0.382109 -0.000515269, 0.000297659, 0.00184, 0.49294, 1.3747775 -0.0000121025, 0.000307953, 0.00208, 0.53260, 0.0316204 -0.0000855093, 0.000299579, 0.00188, 0.50022, 0.2271729 -0.0000200236, 0.000306074, 0.00204, 0.52582, 0.0525426 -0.000107752, 0.000283677, 0.00197, 0.52149, 0.3078042 -0.0000207575, 0.000304806, 0.00201, 0.52089, 0.0546116 -0.0000809795, 0.000281074, 0.00191, 0.51232, 0.2327364 -0.00000442291, 0.000303868, 0.00198, 0.51717, 0.0116488 -0.0000162654, 0.000281067, 0.00191, 0.51235, 0.0467482 -0.00000463438, 0.00033962, 0.00178, 0.46921, 0.0104706 -0.00000956898, 0.000351681, 0.00199, 0.49986, 0.0211592 -0.00000578307, 0.000336363, 0.00172, 0.45907, 0.0131355 -0.0000109853, 0.000266092, 0.00154, 0.42213, 0.0327014 -0.00000291072, 0.000344777, 0.00187, 0.48297, 0.006517 -0.0000104814, 0.000348052, 0.00193, 0.49105, 0.0233399 -0.0000051388, 0.000351427, 0.00199, 0.49888, 0.0113753 -0.0000098941, 0.000346998, 0.00191, 0.48874, 0.022068 -0.00000723814, 0.000349723, 0.00196, 0.49543, 0.0160699 -0.0000592702, 0.000279954, 0.00175, 0.45382, 0.1696584 -0.00000020753, 0.000283367, 0.0018, 0.46005, 0.0005883 -0.00000284117, 0.000275123, 0.00168, 0.44440, 0.0082463 +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 +GCR_Li, 3, 7, 0.000184794, 0.000389118, 0.0021, 0.49406, 0.3655716 +GCR_Be, 4, 9, 0.0000953693, 0.000386546, 0.002, 0.48222, 0.1883217 +GCR_B, 5, 10, 0.000245835, 0.000383147, 0.00182, 0.45495, 0.4814399 +GCR_C, 6, 12, 0.000551454, 0.000297607, 0.00184, 0.49310, 1.4715932 +GCR_N, 7, 14, 0.000177397, 0.000347376, 0.00161, 0.43464, 0.382109 +GCR_O, 8, 16, 0.000515269, 0.000297659, 0.00184, 0.49294, 1.3747775 +GCR_F, 9, 19, 0.0000121025, 0.000307953, 0.00208, 0.53260, 0.0316204 +GCR_Ne, 10, 20, 0.0000855093, 0.000299579, 0.00188, 0.50022, 0.2271729 +GCR_Na, 11, 23, 0.0000200236, 0.000306074, 0.00204, 0.52582, 0.0525426 +GCR_Mg, 12, 24, 0.000107752, 0.000283677, 0.00197, 0.52149, 0.3078042 +GCR_Al, 13, 27, 0.0000207575, 0.000304806, 0.00201, 0.52089, 0.0546116 +GCR_Si, 14, 28, 0.0000809795, 0.000281074, 0.00191, 0.51232, 0.2327364 +GCR_P, 15, 31, 0.00000442291, 0.000303868, 0.00198, 0.51717, 0.0116488 +GCR_S, 16, 32, 0.0000162654, 0.000281067, 0.00191, 0.51235, 0.0467482 +GCR_Cl, 17, 36, 0.00000463438, 0.00033962, 0.00178, 0.46921, 0.0104706 +GCR_Ar, 18, 40, 0.00000956898, 0.000351681, 0.00199, 0.49986, 0.0211592 +GCR_K, 19, 39, 0.00000578307, 0.000336363, 0.00172, 0.45907, 0.0131355 +GCR_Ca, 20, 40, 0.0000109853, 0.000266092, 0.00154, 0.42213, 0.0327014 +GCR_Sc, 21, 45, 0.00000291072, 0.000344777, 0.00187, 0.48297, 0.006517 +GCR_Ti, 22, 48, 0.0000104814, 0.000348052, 0.00193, 0.49105, 0.0233399 +GCR_V, 23, 51, 0.0000051388, 0.000351427, 0.00199, 0.49888, 0.0113753 +GCR_Cr, 24, 52, 0.0000098941, 0.000346998, 0.00191, 0.48874, 0.022068 +GCR_Mn, 25, 55, 0.00000723814, 0.000349723, 0.00196, 0.49543, 0.0160699 +GCR_Fe, 26, 56, 0.0000592702, 0.000279954, 0.00175, 0.45382, 0.1696584 +GCR_Co, 27, 59, 0.00000020753, 0.000283367, 0.0018, 0.46005, 0.0005883 +GCR_Ni, 28, 58, 0.00000284117, 0.000275123, 0.00168, 0.44440, 0.0082463 diff --git a/include/ActionInitialization.h b/include/ActionInitialization.h index 791b0cd..8225b3d 100644 --- a/include/ActionInitialization.h +++ b/include/ActionInitialization.h @@ -3,8 +3,6 @@ #include "G4VUserActionInitialization.hh" -extern G4String particleType; - class ActionInitialization : public G4VUserActionInitialization { public: ActionInitialization(); diff --git a/include/Particle.h b/include/Particle.h index 77ffb41..1752809 100644 --- a/include/Particle.h +++ b/include/Particle.h @@ -8,35 +8,39 @@ class Particle { public: - Particle(G4String name, G4int Z, G4int A) { + 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->name = name; 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); + virtual G4double pdf(G4double E) { return 0; }; + +public: + G4ParticleDefinition* ion; protected: - G4ParticleDefinition* ion; G4int Z, A; - G4String name; G4double sum; G4double C1, C2, C3, C4; }; class TProton : public Particle { public: - TProton(); - G4double pdf(G4double E) { + 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; } @@ -44,7 +48,7 @@ public: class GCRProton : public Particle { public: - GCRProton(); + using Particle::Particle; G4double pdf(G4double E) { G4double pc = getPc(E); G4double beta = getBeta(E); @@ -62,7 +66,7 @@ private: class GCRIon : public Particle { public: - GCRIon(); + using Particle::Particle; G4double pdf(G4double E) { return C1 * exp(-C2 * E) * (1 - exp(-C3 * E + C4)) / sum; } }; diff --git a/include/PrimaryGeneratorAction.h b/include/PrimaryGeneratorAction.h index 6680458..1ff5fa8 100644 --- a/include/PrimaryGeneratorAction.h +++ b/include/PrimaryGeneratorAction.h @@ -1,6 +1,8 @@ #ifndef DESCSS_PrimaryGeneratorAction_h #define DESCSS_PrimaryGeneratorAction_h +#include "Particle.h" + #include "G4GeneralParticleSource.hh" #include "G4ParticleGun.hh" #include "G4VUserPrimaryGeneratorAction.hh" @@ -10,6 +12,8 @@ class G4GeneralParticleSource; class G4ParticleGun; class G4Event; +extern G4String particleType; + template class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction { public: @@ -20,7 +24,6 @@ public: const T* GetParticleGun() const { return fParticleGun; } private: - G4String particleType; T* fParticleGun = nullptr; }; diff --git a/src/ActionInitialization.cpp b/src/ActionInitialization.cpp index ae36d1f..e141cdd 100644 --- a/src/ActionInitialization.cpp +++ b/src/ActionInitialization.cpp @@ -4,8 +4,7 @@ #include "Randomize.hh" -G4String particleType; - +extern G4String particleType; class G4GeneralParticleSource; class G4ParticleGun; diff --git a/src/PrimaryGeneratorAction.cpp b/src/PrimaryGeneratorAction.cpp index c083a63..d76c726 100644 --- a/src/PrimaryGeneratorAction.cpp +++ b/src/PrimaryGeneratorAction.cpp @@ -1,15 +1,50 @@ +#include "Particle.h" #include "PrimaryGeneratorAction.h" #include "G4PhysicalConstants.hh" #include "G4SystemOfUnits.hh" #include "Randomize.hh" +Particle fParticle; +G4String particleType; +G4double Z, A; +G4double C1, C2, C3, C4, sum; +std::string tmp, line, name; +std::ifstream modelFile("F:/Project/Geant4/DESCSS/assets/model.txt"); + PrimaryGeneratorAction::PrimaryGeneratorAction() { fParticleGun = new G4GeneralParticleSource(); } PrimaryGeneratorAction::PrimaryGeneratorAction() { + while (std::getline(modelFile, line)) { + std::stringstream ss(line); + std::getline(ss, name, ','); + if (name == particleType) { + std::getline(ss, tmp, ','), Z = std::stof(tmp); + std::getline(ss, tmp, ','), A = std::stof(tmp); + 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") { + std::getline(ss, tmp, ','), sum = std::stof(tmp); + fParticle = TProton(Z, A); + } 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; + } + } + fParticleGun = new G4ParticleGun(); + fParticleGun->SetParticleDefinition(fParticle.ion); } template @@ -17,11 +52,6 @@ PrimaryGeneratorAction::~PrimaryGeneratorAction() { delete fParticleGun; } -std::string line; -G4double M = 931.5; -G4double C1, C2, C3, C4, sum; -std::ifstream modelFile("F:/Project/Geant4/DESCSS/assets/model.txt"); - G4double randomPhi() { G4double v = G4UniformRand(); return 2 * pi * v; @@ -69,15 +99,22 @@ G4double ARM(G4double Emin, G4double Emax, G4double (*f)(G4double)) { return x; } +G4double pdf(G4double E) { return fParticle.pdf(E); } + void PrimaryGeneratorAction::GeneratePrimaries(G4Event* e) { fParticleGun->GeneratePrimaryVertex(e); } void PrimaryGeneratorAction::GeneratePrimaries(G4Event* e) { - G4double E = 0; + G4double E; G4double R = 15. * m; - G4ThreeVector pos = G4ThreeVector(); // randomPos(R); - G4ThreeVector dir = G4ThreeVector(); // randomDir(); + G4ThreeVector pos = randomPos(R); + G4ThreeVector dir = randomDir(); + + if (particleType == "TP") + E = ARM(50, 1000, pdf); + else + E = ARM(220, 1e5, pdf); fParticleGun->SetParticleEnergy(E); fParticleGun->SetParticlePosition(pos);