diff --git a/G4.code-workspace b/G4.code-workspace index e447f0e..3c470bf 100644 --- a/G4.code-workspace +++ b/G4.code-workspace @@ -8,7 +8,9 @@ "files.associations": { "*.icc": "cpp", "line": "cpp", - "plane": "cpp" + "plane": "cpp", + "qrot": "cpp", + "cmath": "cpp" } } } \ No newline at end of file diff --git a/README.md b/README.md index e289b0a..9041dd6 100644 --- a/README.md +++ b/README.md @@ -185,26 +185,7 @@ $$ -### 电子 -由于该模型较为简单，使用`Geant4`的`General Particle Source`（`GPS`）可以较为方便地生成源， -```macro -/gps/source/clear -/gps/source/add 1 -/gps/particle e- -/gps/pos/type Surface -/gps/pos/shape Sphere -/gps/pos/radius 15 m -/gps/ene/type Exp -/gps/ene/min 100 keV -/gps/ene/max 10 MeV -/gps/ene/ezero 0.0824 -/gps/ang/type iso -/gps/ang/maxtheta 90 deg -/run/beamOn 923242 -``` - -### 质子与俘获辐射 -由于模型较为复杂，且需要发射的粒子种类繁多，因此该部分的粒子使用`G4ParticleGun`生成，关键是生成各向同性（在球面上均匀分布）的位置、方向和符合能谱的能量分布。 +由于模型较为复杂，且需要发射的粒子种类繁多，因此使用`G4ParticleGun`生成各类粒子，关键是生成各向同性（在球面上均匀分布）的位置、方向和符合能谱的能量分布。 1. `Inverse Transform Method`（`ITM`） 对于一个在$[a, b]$的分布，设其概率密度函数(PDF)为$F(x)$，对其积分可获得累积分布函数(CDF)，记为$y=C(x)=\int_a^x F(u)\mathrm{d}u$，其反函数为$C^{-1}(y)$。设$y$为$[0, 1)$上均匀分布的随机数，则$x=C^{-1}(y)$满足$F(x)$的分布。 @@ -236,6 +217,15 @@ $$ 为了保证发射方向是朝向球面内，因此需要限制$\theta$的最大角度为$90^{\circ}$。 3. 能量 +* 俘获电子的函数较为简单，采用`ITM`；（下表中$j_0$已经归一化） +
+ + | $F(E)$ | $\int F(E)\mathrm{d}E$ | $C(x)=\int_a^xF(u)\mathrm{d}u$ | $C^{-1}(y)$ | + |:--------------------------:|:------------------------------------------:|:----------------------------------------------------------------------------:|:---------------------------------------------------------------------------------------------------:| + | $j=j_0e^{-E/E_0}$ | $-j_0E_0e^{-E/E_0}$ | $j_0E_0e^{-a/E_0}-j_0E_0e^{-x/E_0}$ | $-E_0\ln(\frac{C-y}{j_0E_0}),\ C=j_0E_0e^{-a/E_0}$ | + +
+
* 俘获质子的函数较为简单，采用`ITM`；（下表中$j_0$已经归一化）
diff --git a/assets/model.txt b/assets/model.txt index a507c08..d27adb9 100644 --- a/assets/model.txt +++ b/assets/model.txt @@ -1,3 +1,4 @@ +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 diff --git a/docs/simu-te.png b/docs/simu-te.png new file mode 100644 index 0000000..48eb078 Binary files /dev/null and b/docs/simu-te.png differ diff --git a/include/G4CustomFemaleBuilder.h b/include/G4CustomFemaleBuilder.h deleted file mode 100644 index 562c8bf..0000000 --- a/include/G4CustomFemaleBuilder.h +++ /dev/null @@ -1,17 +0,0 @@ -#ifndef G4CustomFemaleBuilder_h -#define G4CustomFemaleBuilder_h 1 - -#include "G4PhantomBuilder.h" - -class G4CustomFemaleBuilder : public G4PhantomBuilder { -public: - G4CustomFemaleBuilder(); - ~G4CustomFemaleBuilder(); - - void BuildVoxelLeftBreast(const G4String&, G4bool, G4bool); - void BuildVoxelRightBreast(const G4String&, G4bool, G4bool); - void BuildLeftOvary(const G4String&, G4bool, G4bool); - void BuildRightOvary(const G4String&, G4bool, G4bool); - void BuildUterus(const G4String&, G4bool, G4bool); -}; -#endif diff --git a/include/G4HumanPhantomMaterial.h b/include/G4HumanPhantomMaterial.h index 221fe9c..2676361 100644 --- a/include/G4HumanPhantomMaterial.h +++ b/include/G4HumanPhantomMaterial.h @@ -14,12 +14,6 @@ public: G4Material* GetMaterial(G4String); // returns the material private: - // G4Material* matW; - // G4Material* matplexiglass; - // G4Material* matPb; - // G4Material* matir192; - // G4Material* Titanium; - // G4Material* matAir; G4Material* matH2O; G4Material* soft; G4Material* skeleton; @@ -27,8 +21,6 @@ private: G4Material* adipose; G4Material* glandular; G4Material* adipose_glandular; - // G4Material*Vacuum; - - // G4Material* muscle; + G4Material* muscle; }; #endif diff --git a/include/G4HumanPhantomMessenger.h b/include/G4HumanPhantomMessenger.h index de6cd3b..fc24452 100644 --- a/include/G4HumanPhantomMessenger.h +++ b/include/G4HumanPhantomMessenger.h @@ -20,21 +20,14 @@ public: void SetNewValue(G4UIcommand* command, G4String newValue); - void AddBodyPart(G4String); // Set Body Parts Sensitivity - private: DetectorConstruction* myUserPhantom; G4UIdirectory* phantomDir; - G4UIdirectory* bpDir; G4UIcmdWithAString* modelCmd; G4UIcmdWithAString* sexCmd; - G4UIcmdWithAString* bodypartCmd; G4UIcmdWithoutParameter* endCmd; - - G4String bodypart; - G4bool bps; }; #endif diff --git a/include/G4VoxelBreastFactory.h b/include/G4VoxelBreastFactory.h deleted file mode 100644 index fbabc33..0000000 --- a/include/G4VoxelBreastFactory.h +++ /dev/null @@ -1,21 +0,0 @@ -#ifndef G4VoxelBreastFactory_h -#define G4VoxelBreastFactory_h 1 - -#include "G4VBodyFactory.h" -#include "G4VOrgan.h" - -#include -class G4VBodyFactory; -class G4VPhysicalVolume; -class G4VOrgan; - -class G4VoxelBreastFactory : public G4VBodyFactory { -public: - G4VoxelBreastFactory(); - ~G4VoxelBreastFactory(); - G4VPhysicalVolume* CreateOrgan(const G4String&, G4VPhysicalVolume*, const G4String&, G4bool, G4bool); - -private: - std::map organ; -}; -#endif diff --git a/include/G4VoxelLeftBreast.h b/include/G4VoxelLeftBreast.h deleted file mode 100644 index e10cd4b..0000000 --- a/include/G4VoxelLeftBreast.h +++ /dev/null @@ -1,16 +0,0 @@ -#ifndef G4VoxelLeftBreast_h -#define G4VoxelLeftBreast_h 1 - -#include "G4VOrgan.h" - -#include "G4VPhysicalVolume.hh" - -class G4VPhysicalVolume; - -class G4VoxelLeftBreast : public G4VOrgan { -public: - G4VoxelLeftBreast(); - ~G4VoxelLeftBreast(); - G4VPhysicalVolume* Construct(const G4String&, G4VPhysicalVolume*, const G4String&, G4bool, G4bool); -}; -#endif diff --git a/include/G4VoxelRightBreast.h b/include/G4VoxelRightBreast.h deleted file mode 100644 index 477f188..0000000 --- a/include/G4VoxelRightBreast.h +++ /dev/null @@ -1,16 +0,0 @@ -#ifndef G4VoxelRightBreast_h -#define G4VoxelRightBreast_h 1 - -#include "G4VOrgan.h" - -#include "G4VPhysicalVolume.hh" - -class G4VPhysicalVolume; - -class G4VoxelRightBreast : public G4VOrgan { -public: - G4VoxelRightBreast(); - ~G4VoxelRightBreast(); - G4VPhysicalVolume* Construct(const G4String&, G4VPhysicalVolume*, const G4String&, G4bool, G4bool); -}; -#endif diff --git a/include/PrimaryGeneratorAction.h b/include/PrimaryGeneratorAction.h index ebdd835..3abed07 100644 --- a/include/PrimaryGeneratorAction.h +++ b/include/PrimaryGeneratorAction.h @@ -1,28 +1,25 @@ #ifndef DESCSS_PrimaryGeneratorAction_h #define DESCSS_PrimaryGeneratorAction_h -#include "G4GeneralParticleSource.hh" #include "G4ParticleGun.hh" #include "G4VUserPrimaryGeneratorAction.hh" #include "globals.hh" -class G4GeneralParticleSource; class G4ParticleGun; class G4Event; extern G4String particleType; -template class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction { public: PrimaryGeneratorAction(); ~PrimaryGeneratorAction(); virtual void GeneratePrimaries(G4Event*); - const T* GetParticleGun() const { return fParticleGun; } + const G4ParticleGun* GetParticleGun() const { return fParticleGun; } private: - T* fParticleGun = nullptr; + G4ParticleGun* fParticleGun = nullptr; }; #endif diff --git a/src/ActionInitialization.cpp b/src/ActionInitialization.cpp index e141cdd..a9822a1 100644 --- a/src/ActionInitialization.cpp +++ b/src/ActionInitialization.cpp @@ -4,7 +4,6 @@ #include "Randomize.hh" -extern G4String particleType; class G4GeneralParticleSource; class G4ParticleGun; @@ -19,8 +18,5 @@ ActionInitialization::~ActionInitialization() {} void ActionInitialization::BuildForMaster() const {} void ActionInitialization::Build() const { - if (particleType == "TE") - SetUserAction(new PrimaryGeneratorAction); - else - SetUserAction(new PrimaryGeneratorAction); + SetUserAction(new PrimaryGeneratorAction); } diff --git a/src/DetectorConstruction.cpp b/src/DetectorConstruction.cpp index 25d5e51..0f48dbf 100644 --- a/src/DetectorConstruction.cpp +++ b/src/DetectorConstruction.cpp @@ -1,6 +1,5 @@ #include "DetectorConstruction.h" #include "G4BasePhantomBuilder.h" -#include "G4CustomFemaleBuilder.h" #include "G4FemaleBuilder.h" #include "G4HumanPhantomMaterial.h" #include "G4HumanPhantomMessenger.h" diff --git a/src/G4CustomFemaleBuilder.cpp b/src/G4CustomFemaleBuilder.cpp deleted file mode 100644 index 5850890..0000000 --- a/src/G4CustomFemaleBuilder.cpp +++ /dev/null @@ -1,62 +0,0 @@ -#include "G4CustomFemaleBuilder.h" -#include "G4VBodyFactory.h" -#include "G4VoxelBreastFactory.h" - -G4CustomFemaleBuilder::G4CustomFemaleBuilder() {} - -G4CustomFemaleBuilder::~G4CustomFemaleBuilder() - -{ - delete body; -} - -void G4CustomFemaleBuilder::BuildUterus(const G4String& colourName, G4bool solidVis, G4bool sensitivity) - -{ - if (trunkVolume == 0) - G4Exception("G4CustomFemaleBuilder::BuildUterus()", "human_phantom0001", FatalException, - "The trunk volume is missing !!!!!"); - - body->CreateOrgan("Uterus", trunkVolume, colourName, solidVis, sensitivity); -} - -void G4CustomFemaleBuilder::BuildLeftOvary(const G4String& colourName, G4bool solidVis, G4bool sensitivity) - -{ - if (trunkVolume == 0) - G4Exception("G4CustomFemaleBuilder::BuildLeftOvary()", "human_phantom0002", FatalException, - "The trunk volume is missing !!!!!"); - - body->CreateOrgan("LeftOvary", trunkVolume, colourName, solidVis, sensitivity); -} - -void G4CustomFemaleBuilder::BuildRightOvary(const G4String& colourName, G4bool solidVis, G4bool sensitivity) - -{ - if (trunkVolume == 0) - G4Exception("G4CustomFemaleBuilder::BuildRightOvary()", "human_phantom0003", FatalException, - "The trunk volume is missing !!!!!"); - - body->CreateOrgan("RightOvary", trunkVolume, colourName, solidVis, sensitivity); -} - -void G4CustomFemaleBuilder::BuildVoxelLeftBreast(const G4String& colourName, G4bool solidVis, G4bool sensitivity) { - G4cout << "BuildVoxelLeftBreast" << G4endl; - if (motherVolume == 0) - G4Exception("G4CustomFemaleBuilder::BuildVoxelLeftBreast()", "human_phantom0004", FatalException, - "The world volume is missing !!!!!"); - - G4VBodyFactory* customBody = new G4VoxelBreastFactory(); - customBody->CreateOrgan("LeftBreast", motherVolume, colourName, solidVis, sensitivity); - delete customBody; -} - -void G4CustomFemaleBuilder::BuildVoxelRightBreast(const G4String& colourName, G4bool solidVis, G4bool sensitivity) { - if (motherVolume == 0) - G4Exception("G4CustomFemaleBuilder::BuildVoxelRightBreast()", "human_phantom0005", FatalException, - "The world volume is missing !!!!!"); - - G4VBodyFactory* customBody = new G4VoxelBreastFactory(); - customBody->CreateOrgan("RightBreast", motherVolume, colourName, solidVis, sensitivity); - delete customBody; -} diff --git a/src/G4HumanPhantomMaterial.cpp b/src/G4HumanPhantomMaterial.cpp index 120f892..c530481 100644 --- a/src/G4HumanPhantomMaterial.cpp +++ b/src/G4HumanPhantomMaterial.cpp @@ -170,12 +170,6 @@ void G4HumanPhantomMaterial::DefineMaterials() { adipose_glandular = new G4Material("adipose_glandular", d, 2); adipose_glandular->AddMaterial(adipose, 0.5); adipose_glandular->AddMaterial(glandular, 0.5); - - // Air - d = 1.290 * mg / cm3; - G4Material* matAir = new G4Material("Air", d, 2); - matAir->AddElement(elN, 0.7); - matAir->AddElement(elO, 0.3); } G4Material* G4HumanPhantomMaterial::GetMaterial(G4String material) { diff --git a/src/G4HumanPhantomMessenger.cpp b/src/G4HumanPhantomMessenger.cpp index 870e211..bd9658b 100644 --- a/src/G4HumanPhantomMessenger.cpp +++ b/src/G4HumanPhantomMessenger.cpp @@ -7,14 +7,10 @@ #include "G4UIdirectory.hh" #include "globals.hh" -G4HumanPhantomMessenger::G4HumanPhantomMessenger(DetectorConstruction* myUsrPhtm) - : myUserPhantom(myUsrPhtm), bps(false) { +G4HumanPhantomMessenger::G4HumanPhantomMessenger(DetectorConstruction* myUsrPhtm) : myUserPhantom(myUsrPhtm) { phantomDir = new G4UIdirectory("/phantom/"); phantomDir->SetGuidance("Set Your Phantom."); - bpDir = new G4UIdirectory("/bodypart/"); - bpDir->SetGuidance("Add Body Part to Phantom"); - modelCmd = new G4UIcmdWithAString("/phantom/setPhantomModel", this); modelCmd->SetGuidance("Set sex of Phantom: MIRD, MIRDHead."); modelCmd->SetParameterName("phantomModel", true); @@ -29,11 +25,6 @@ G4HumanPhantomMessenger::G4HumanPhantomMessenger(DetectorConstruction* myUsrPhtm sexCmd->SetCandidates("Male Female"); sexCmd->AvailableForStates(G4State_PreInit, G4State_Idle); - bodypartCmd = new G4UIcmdWithAString("/bodypart/addBodyPart", this); - bodypartCmd->SetGuidance("Add a Body Part to Phantom"); - bodypartCmd->SetParameterName("bpName", true); - bodypartCmd->AvailableForStates(G4State_PreInit, G4State_Idle); - endCmd = new G4UIcmdWithoutParameter("/phantom/buildNewPhantom", this); endCmd->SetGuidance("Build your Phantom."); endCmd->AvailableForStates(G4State_PreInit, G4State_Idle); @@ -42,10 +33,8 @@ G4HumanPhantomMessenger::G4HumanPhantomMessenger(DetectorConstruction* myUsrPhtm G4HumanPhantomMessenger::~G4HumanPhantomMessenger() { delete modelCmd; delete sexCmd; - delete bodypartCmd; delete endCmd; delete phantomDir; - delete bpDir; } void G4HumanPhantomMessenger::SetNewValue(G4UIcommand* command, G4String newValue) { @@ -55,31 +44,7 @@ void G4HumanPhantomMessenger::SetNewValue(G4UIcommand* command, G4String newValu if (command == sexCmd) { myUserPhantom->SetPhantomSex(newValue); } - if (command == bodypartCmd) { - AddBodyPart(newValue); - } if (command == endCmd) { G4cout << " ****************>>>> NEW PHANTOM CONSTRUCTION <<<<***************** " << G4endl; } } - -void G4HumanPhantomMessenger::AddBodyPart(G4String newBodyPartSensitivity) { - char* str = new char[newBodyPartSensitivity.length() + 1]; - - strcpy(str, newBodyPartSensitivity.c_str()); - - std::string bpart = strtok(str, " "); - - std::string sensitivity = strtok(NULL, " "); - - if (sensitivity == "yes") { - bps = true; - } else { - bps = false; - } - - G4cout << " >>> Body Part = " << bpart << "\n" - << " >>> Sensitivity = " << sensitivity << G4endl; - - myUserPhantom->SetBodyPartSensitivity(bpart, bps); -} diff --git a/src/G4MIRDBodyFactory.cpp b/src/G4MIRDBodyFactory.cpp index b43c2f6..778e418 100644 --- a/src/G4MIRDBodyFactory.cpp +++ b/src/G4MIRDBodyFactory.cpp @@ -1,5 +1,4 @@ #include "G4MIRDBodyFactory.h" - #include "G4MIRDBrain.h" #include "G4MIRDHead.h" #include "G4MIRDHeart.h" diff --git a/src/G4VoxelBreastFactory.cpp b/src/G4VoxelBreastFactory.cpp deleted file mode 100644 index e64740c..0000000 --- a/src/G4VoxelBreastFactory.cpp +++ /dev/null @@ -1,24 +0,0 @@ -#include "G4VoxelBreastFactory.h" -#include "G4VoxelLeftBreast.h" -#include "G4VoxelRightBreast.h" - -G4VoxelBreastFactory::G4VoxelBreastFactory() { - // Map with name of the organ and pointer to the MIRDOrgan class - // organ["ParameterisedRightBreast"] = new G4ParameterisedRightBreast(); - organ["LeftBreast"] = new G4VoxelLeftBreast(); - organ["RightBreast"] = new G4VoxelRightBreast(); -} - -G4VoxelBreastFactory::~G4VoxelBreastFactory() { - delete organ["RightBreast"]; - organ["RightBreast"] = 0; - - delete organ["LeftBreast"]; - organ["LeftBreast"] = 0; -} - -G4VPhysicalVolume* G4VoxelBreastFactory::CreateOrgan(const G4String& organ_name, G4VPhysicalVolume* motherVolume, - const G4String& colourName, G4bool visAttribute, - G4bool sensitivity) { - return organ[organ_name]->Construct(organ_name, motherVolume, colourName, visAttribute, sensitivity); -} diff --git a/src/G4VoxelLeftBreast.cpp b/src/G4VoxelLeftBreast.cpp deleted file mode 100644 index c7a3703..0000000 --- a/src/G4VoxelLeftBreast.cpp +++ /dev/null @@ -1,72 +0,0 @@ -#include "G4HumanPhantomColour.h" -#include "G4HumanPhantomMaterial.h" -#include "G4VoxelLeftBreast.h" - -#include "G4LogicalVolume.hh" -#include "G4Material.hh" -#include "G4PVPlacement.hh" -#include "G4PVReplica.hh" -#include "G4RotationMatrix.hh" -#include "G4SDManager.hh" -#include "G4SystemOfUnits.hh" -#include "G4ThreeVector.hh" -#include "G4Tubs.hh" -#include "G4VPhysicalVolume.hh" -#include "G4VisAttributes.hh" -#include "globals.hh" - -G4VoxelLeftBreast::G4VoxelLeftBreast() {} - -G4VoxelLeftBreast::~G4VoxelLeftBreast() {} - -G4VPhysicalVolume* G4VoxelLeftBreast::Construct(const G4String& volumeName, G4VPhysicalVolume* mother, - const G4String& colourName, G4bool wireFrame, G4bool) - -{ - G4HumanPhantomMaterial* material = new G4HumanPhantomMaterial(); - G4Material* adipose = material->GetMaterial("adipose"); - G4Material* adipose_glandular = material->GetMaterial("adipose_glandular"); - delete material; - - G4double rmin = 0. * cm; - G4double rmax = 6. * cm; - G4double zz = 6. * cm; - G4double startPhi = 0. * degree; - G4double spanningPhi = 180. * degree; - - G4Tubs* breast = new G4Tubs("out_breast", rmin, rmax, zz / 2., startPhi, spanningPhi); - - G4LogicalVolume* breast_log = new G4LogicalVolume(breast, adipose, "logicalOut" + volumeName, 0, 0, 0); - rmax = 5.5 * cm; - zz = 5. * cm; - G4Tubs* innerBreast = new G4Tubs("inner_breast", rmin, rmax, zz / 2., startPhi, spanningPhi); - - G4LogicalVolume* innerBreast_log = - new G4LogicalVolume(innerBreast, adipose_glandular, "logical" + volumeName, 0, 0, 0); - - G4RotationMatrix* matrix = new G4RotationMatrix(); - // matrix -> rotateX(-90.* degree); - // matrix -> rotateY(180.* degree); - matrix->rotateZ(18. * degree); - - G4VPhysicalVolume* physBreast = new G4PVPlacement(matrix, G4ThreeVector(10. * cm, 8.7 * cm, 52. * cm), - "physicalVoxelLeftBreast", breast_log, mother, false, 0, true); - - G4VPhysicalVolume* physInnerBreast = - new G4PVPlacement(0, G4ThreeVector(), "LeftBreast", innerBreast_log, physBreast, false, 0, true); - - G4HumanPhantomColour* colourPointer = new G4HumanPhantomColour(); - G4Colour colour = colourPointer->GetColour(colourName); - delete colourPointer; - - G4VisAttributes* BreastVisAtt = new G4VisAttributes(colour); - BreastVisAtt->SetForceSolid(wireFrame); - breast_log->SetVisAttributes(BreastVisAtt); - - G4VisAttributes* innerBreastVisAtt = new G4VisAttributes(colour); - innerBreastVisAtt->SetForceSolid(false); - innerBreast_log->SetVisAttributes(innerBreastVisAtt); - - G4cout << "Voxel Left Breast created !!!!!! This model has to be refined" << G4endl; - return physInnerBreast; -} diff --git a/src/G4VoxelRightBreast.cpp b/src/G4VoxelRightBreast.cpp deleted file mode 100644 index ff0a602..0000000 --- a/src/G4VoxelRightBreast.cpp +++ /dev/null @@ -1,74 +0,0 @@ -#include "G4HumanPhantomColour.h" -#include "G4HumanPhantomMaterial.h" -#include "G4VoxelRightBreast.h" - -#include "G4LogicalVolume.hh" -#include "G4Material.hh" -#include "G4PVPlacement.hh" -#include "G4PVReplica.hh" -#include "G4RotationMatrix.hh" -#include "G4SDManager.hh" -#include "G4SystemOfUnits.hh" -#include "G4ThreeVector.hh" -#include "G4Tubs.hh" -#include "G4VPhysicalVolume.hh" -#include "G4VisAttributes.hh" -#include "globals.hh" - -G4VoxelRightBreast::G4VoxelRightBreast() {} - -G4VoxelRightBreast::~G4VoxelRightBreast() {} - -G4VPhysicalVolume* G4VoxelRightBreast::Construct(const G4String& volumeName, G4VPhysicalVolume* mother, - const G4String& colourName, G4bool wireFrame, G4bool) - -{ - G4cout << "Construct " << volumeName << " with mother volume " << mother->GetName() << G4endl; - G4HumanPhantomMaterial* material = new G4HumanPhantomMaterial(); - G4Material* adipose = material->GetMaterial("adipose"); - G4Material* adipose_glandular = material->GetMaterial("adipose_glandular"); - delete material; - - G4double rmin = 0. * cm; - G4double rmax = 6. * cm; - G4double zz = 6. * cm; - G4double startPhi = 0. * degree; - G4double spanningPhi = 180. * degree; - - G4Tubs* breast = new G4Tubs("out_breast", rmin, rmax, zz / 2., startPhi, spanningPhi); - - G4LogicalVolume* breast_log = new G4LogicalVolume(breast, adipose, "logicalOut" + volumeName, 0, 0, 0); - rmax = 5.5 * cm; - zz = 5. * cm; - G4Tubs* innerBreast = new G4Tubs("inner_breast", rmin, rmax, zz / 2., startPhi, spanningPhi); - - G4LogicalVolume* innerBreast_log = - new G4LogicalVolume(innerBreast, adipose_glandular, "logical" + volumeName, 0, 0, 0); - - G4RotationMatrix* matrix = new G4RotationMatrix(); - matrix->rotateX(0. * degree); - matrix->rotateY(0. * degree); - matrix->rotateZ(-18. * degree); - - G4VPhysicalVolume* physBreast = new G4PVPlacement(matrix, G4ThreeVector(-10. * cm, 8.7 * cm, 52. * cm), - "physicalVoxelRightBreast", breast_log, mother, false, 0, true); - - G4VPhysicalVolume* physInnerBreast = - new G4PVPlacement(0, G4ThreeVector(), "RightBreast", innerBreast_log, physBreast, false, 0, true); - - G4HumanPhantomColour* colourPointer = new G4HumanPhantomColour(); - G4Colour colour = colourPointer->GetColour(colourName); - delete colourPointer; - - G4VisAttributes* BreastVisAtt = new G4VisAttributes(colour); - BreastVisAtt->SetForceSolid(wireFrame); - breast_log->SetVisAttributes(BreastVisAtt); - - G4VisAttributes* innerBreastVisAtt = new G4VisAttributes(colour); - innerBreastVisAtt->SetForceSolid(false); - innerBreast_log->SetVisAttributes(innerBreastVisAtt); - - G4cout << "Voxel Right Breast created !!!!!! This model must be refined!" << G4endl; - - return physInnerBreast; -} diff --git a/src/PrimaryGeneratorAction.cpp b/src/PrimaryGeneratorAction.cpp index 6da5a8e..7bf8dab 100644 --- a/src/PrimaryGeneratorAction.cpp +++ b/src/PrimaryGeneratorAction.cpp @@ -14,11 +14,7 @@ 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() { +PrimaryGeneratorAction::PrimaryGeneratorAction() { while (std::getline(modelFile, line)) { std::stringstream ss(line); std::getline(ss, name, ','); @@ -28,7 +24,9 @@ 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 == "TE") + sum = C3; + else if (tmp == "TP") std::getline(ss, tmp, ','), sum = std::stof(tmp); else { std::getline(ss, tmp, ','), C4 = std::stof(tmp); @@ -46,14 +44,11 @@ PrimaryGeneratorAction::PrimaryGeneratorAction() { ion = G4ParticleTable::GetParticleTable()->FindParticle("alpha"); else ion = G4IonTable::GetIonTable()->GetIon(Z, A, 0.); - fParticleGun->SetParticleCharge(Z * eplus); + fParticleGun->SetParticleCharge(Z * eplus); fParticleGun->SetParticleDefinition(ion); } -template -PrimaryGeneratorAction::~PrimaryGeneratorAction() { - delete fParticleGun; -} +PrimaryGeneratorAction::~PrimaryGeneratorAction() { delete fParticleGun; } G4double randomPhi() { G4double v = G4UniformRand(); @@ -102,6 +97,11 @@ G4double ARM(G4double Emin, G4double Emax, G4double (*f)(G4double)) { return x; } +G4double pdfTElectronInv(G4double y) { + G4double C = - C2 * C1 * exp(- 0.1 / C2) / sum; + return - C2 * log((C - y) * sum / C1 / C2); +} + 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; @@ -115,17 +115,15 @@ G4double pdfGCRProton(G4double E) { G4double pdfGCRIon(G4double E) { return C1 * exp(-C2 * E) * (1 - exp(-C3 * E + C4)) / sum; } -void PrimaryGeneratorAction::GeneratePrimaries(G4Event* e) { - fParticleGun->GeneratePrimaryVertex(e); -} - -void PrimaryGeneratorAction::GeneratePrimaries(G4Event* e) { +void PrimaryGeneratorAction::GeneratePrimaries(G4Event* e) { G4double E; G4double R = 15. * m; G4ThreeVector pos = randomPos(R); G4ThreeVector dir = randomDir(); - if (particleType == "TP") + if (particleType == "TE") + E = ITM(0.1, 10, pdfTElectronInv); + else if (particleType == "TP") E = ITM(50, 1000, pdfTProtonInv); else if (particleType == "GCR_H") E = ARM(220, 1e5, pdfGCRProton); diff --git a/utils/check.py b/utils/check.py index 386783f..d724ff9 100644 --- a/utils/check.py +++ b/utils/check.py @@ -11,9 +11,6 @@ def trap_proton(): def pdf_int(E, C1=559.76377, C2=-1.53424, C3=2.40873, sum=1.5156509): return C1 * math.pow(E - C2, 1 - C3) / (1 - C3) / sum - def pdf_sum(E): - return pdf_int(E) - pdf_int(50) - def pdf_inv(y, C1=559.76377, C2=-1.53424, C3=2.40873, sum=1.5156509): return pow(sum * (pdf_int(50) + y) * (1 - C3) / C1, 1 / (1 - C3)) + C2