add: SD, custom command

This commit is contained in:
liuyihui 2022-05-17 21:14:44 +08:00
parent e1e512fa18
commit 8da1bce0dc
12 changed files with 403 additions and 57 deletions

View File

@ -6,7 +6,9 @@
],
"settings": {
"files.associations": {
"*.icc": "cpp"
"*.icc": "cpp",
"line": "cpp",
"plane": "cpp"
}
}
}

View File

@ -6,7 +6,6 @@
/run/verbose 0
/event/verbose 0
/tracking/verbose 0
/gps/verbose 0
# 初始化
/run/initialize
@ -22,4 +21,6 @@
/gps/ene/ezero 0.0824
/gps/ang/type iso
/gps/ang/maxtheta 90 deg
/gps/verbose 0
/run/beamOn 923242

View File

@ -7,6 +7,7 @@
class G4VPhysicalVolume;
class G4LogicalVolume;
class G4HumanPhantomMaterial;
class G4HumanPhantomMessenger;
class DetectorConstruction : public G4VUserDetectorConstruction {
public:
@ -15,6 +16,12 @@ public:
G4VPhysicalVolume* Construct() override;
public:
void SetBodyPartSensitivity(G4String, G4bool);
void SetPhantomSex(G4String);
void SetPhantomModel(G4String);
void ConstructSDandField();
private:
void ConstructSectionSphere(G4LogicalVolume* fMotherLogical, G4double zBias);
void ConstructSectionCons(G4String name, G4LogicalVolume* fMotherLogical, G4double zBias, G4double pRmax1,
@ -24,10 +31,11 @@ private:
void ConstructHumanPhantom(G4VPhysicalVolume* fMotherPhysics);
private:
G4HumanPhantomMaterial* material;
G4String model;
G4String sex;
G4String sex = "Male";
G4String model = "MIRD";
std::map<std::string, G4bool> sensitivities;
G4HumanPhantomMaterial* material;
G4HumanPhantomMessenger* messenger;
};
#endif

View File

@ -0,0 +1,49 @@
#ifndef G4HumanPhantomHit_h
#define G4HumanPhantomHit_h 1
#include "G4Allocator.hh"
#include "G4THitsCollection.hh"
#include "G4ThreeVector.hh"
#include "G4VHit.hh"
#include "tls.hh" // FOR MT
class G4HumanPhantomHit : public G4VHit {
public:
G4HumanPhantomHit();
~G4HumanPhantomHit();
G4HumanPhantomHit(const G4HumanPhantomHit&);
const G4HumanPhantomHit& operator=(const G4HumanPhantomHit&);
G4bool operator==(const G4HumanPhantomHit&) const;
inline void* operator new(size_t);
inline void operator delete(void*);
void Draw();
void Print();
public:
void SetBodyPartID(G4String bodyPartName) { bodyPartID = bodyPartName; };
void SetEdep(G4double de) { edep = de; };
G4String GetBodyPartID() { return bodyPartID; };
G4double GetEdep() { return edep; };
private:
G4String bodyPartID;
G4double edep;
};
typedef G4THitsCollection<G4HumanPhantomHit> G4HumanPhantomHitsCollection;
extern G4ThreadLocal G4Allocator<G4HumanPhantomHit>* G4HumanPhantomHitAllocator;
inline void* G4HumanPhantomHit::operator new(size_t) {
if (!G4HumanPhantomHitAllocator) G4HumanPhantomHitAllocator = new G4Allocator<G4HumanPhantomHit>;
return (void*)G4HumanPhantomHitAllocator->MallocSingle();
}
inline void G4HumanPhantomHit::operator delete(void* aHit) {
G4HumanPhantomHitAllocator->FreeSingle((G4HumanPhantomHit*)aHit);
}
#endif

View File

@ -0,0 +1,40 @@
#ifndef G4HumanPhantomMessenger_h
#define G4HumanPhantomMessenger_h 1
class DetectorConstruction;
class G4UIcommand;
class G4UIdirectory;
class G4UIcmdWithAString;
class G4UIcmdWithoutParameter;
#include "G4UImessenger.hh"
#include "globals.hh"
#include <iostream>
class G4HumanPhantomMessenger : public G4UImessenger {
public:
G4HumanPhantomMessenger(DetectorConstruction* myUsrPhtm);
~G4HumanPhantomMessenger();
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

View File

@ -0,0 +1,22 @@
#ifndef G4HumanPhantomSD_h
#define G4HumanPhantomSD_h 1
#include "G4HumanPhantomHit.h"
#include "G4VSensitiveDetector.hh"
class G4Step;
class G4HumanPhantomSD : public G4VSensitiveDetector {
public:
G4HumanPhantomSD(const G4String& name, const G4String& hitsCollectionName);
~G4HumanPhantomSD();
void Initialize(G4HCofThisEvent*);
G4bool ProcessHits(G4Step*, G4TouchableHistory*);
void EndOfEvent(G4HCofThisEvent*);
private:
G4HumanPhantomHitsCollection* collection;
};
#endif

View File

@ -3,6 +3,8 @@
#include "G4CustomFemaleBuilder.h"
#include "G4FemaleBuilder.h"
#include "G4HumanPhantomMaterial.h"
#include "G4HumanPhantomMessenger.h"
#include "G4HumanPhantomSD.h"
#include "G4MaleBuilder.h"
#include "Material.h"
@ -11,6 +13,7 @@
#include "G4IntersectionSolid.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4SDManager.hh"
#include "G4Sphere.hh"
#include "G4SubtractionSolid.hh"
#include "G4SystemOfUnits.hh"
@ -18,9 +21,15 @@
#include "G4Tubs.hh"
#include "G4UnionSolid.hh"
DetectorConstruction::DetectorConstruction() { material = new G4HumanPhantomMaterial(); }
DetectorConstruction::DetectorConstruction() {
material = new G4HumanPhantomMaterial();
messenger = new G4HumanPhantomMessenger(this);
}
DetectorConstruction::~DetectorConstruction() { delete material; }
DetectorConstruction::~DetectorConstruction() {
delete material;
delete messenger;
}
void DetectorConstruction::ConstructSectionSphere(G4LogicalVolume* fMotherLogical, G4double zBias) {
G4double pRadius = 2.8 * m / 2;
@ -226,37 +235,22 @@ void DetectorConstruction::ConstructSectionBig(G4LogicalVolume* fMotherLogical,
void DetectorConstruction::ConstructHumanPhantom(G4VPhysicalVolume* fMotherPhysics) {
material->DefineMaterials();
G4BasePhantomBuilder* builder = 0;
if (sex == "Female") {
if (model == "MIX")
builder = new G4CustomFemaleBuilder;
else {
builder = new G4FemaleBuilder;
}
builder->SetModel(model);
G4cout << model << " " << sex << G4endl;
} else if (sex == "Male") {
if (sex == "Female")
builder = new G4FemaleBuilder;
else
builder = new G4MaleBuilder;
builder->SetModel(model);
if (model == "MIX") {
G4cout << "Custom Male is not available!!! MIRD model is selected !" << G4endl;
model = "MIRD";
builder->SetModel(model);
}
}
builder->SetModel(model);
builder->SetMotherVolume(fMotherPhysics);
// the argument indicates the sensitivity of the volume
builder->BuildHead("black", false, sensitivities["Head"]);
builder->BuildHead("blue", false, sensitivities["Head"]);
builder->BuildSkull("orange", false, sensitivities["Skull"]);
builder->BuildBrain("yellow", true, sensitivities["Brain"]);
if (model != "MIRDHead") {
// builder->SetModel(model);
builder->BuildTrunk("yellow", false, sensitivities["Trunk"]);
builder->BuildLeftLeg("yellow", false, sensitivities["LeftLeg"]);
@ -270,17 +264,15 @@ void DetectorConstruction::ConstructHumanPhantom(G4VPhysicalVolume* fMotherPhysi
builder->BuildUpperSpine("yellow", true, sensitivities["UpperSpine"]);
if (model == "MIRD" || model == "MIX") {
builder->BuildLeftScapula("grey", true, sensitivities["LeftScapula"]);
builder->BuildRightScapula("grey", true, sensitivities["RightScapula"]);
builder->BuildLeftAdrenal("yellow", true, sensitivities["LeftAdrenal"]);
builder->BuildRightAdrenal("yellow", true, sensitivities["RightAdrenal"]);
builder->BuildThymus("orange", true, sensitivities["Thymus"]);
builder->BuildLeftClavicle("grey", true, sensitivities["LeftClavicle"]);
builder->BuildRightClavicle("grey", true, sensitivities["RightClavicle"]);
builder->BuildSmallIntestine("orange", true, sensitivities["SmallIntestine"]);
builder->BuildRibCage("grey", true, sensitivities["RibCage"]);
}
builder->BuildLeftScapula("grey", true, sensitivities["LeftScapula"]);
builder->BuildRightScapula("grey", true, sensitivities["RightScapula"]);
builder->BuildLeftAdrenal("yellow", true, sensitivities["LeftAdrenal"]);
builder->BuildRightAdrenal("yellow", true, sensitivities["RightAdrenal"]);
builder->BuildThymus("orange", true, sensitivities["Thymus"]);
builder->BuildLeftClavicle("grey", true, sensitivities["LeftClavicle"]);
builder->BuildRightClavicle("grey", true, sensitivities["RightClavicle"]);
builder->BuildSmallIntestine("orange", true, sensitivities["SmallIntestine"]);
builder->BuildRibCage("grey", true, sensitivities["RibCage"]);
builder->BuildMiddleLowerSpine("yellow", true, sensitivities["MiddleLowerSpine"]);
@ -311,7 +303,7 @@ void DetectorConstruction::ConstructHumanPhantom(G4VPhysicalVolume* fMotherPhysi
if (model == "MIRD") {
builder->BuildLeftBreast("purple", true, sensitivities["LeftBreast"]);
builder->BuildRightBreast("purple", true, sensitivities["RightBreast"]);
} else if (model == "MIX") {
} else {
builder->BuildVoxelLeftBreast("purple", false, sensitivities["LeftBreast"]);
builder->BuildVoxelRightBreast("purple", false, sensitivities["RightBreast"]);
}
@ -323,7 +315,7 @@ void DetectorConstruction::ConstructHumanPhantom(G4VPhysicalVolume* fMotherPhysi
builder->BuildRightTeste("purple", true, sensitivities["RightTeste"]);
}
}
G4VPhysicalVolume* result = builder->GetPhantom();
delete builder;
}
@ -336,19 +328,6 @@ G4VPhysicalVolume* DetectorConstruction::Construct() {
G4LogicalVolume* logicWorld = new G4LogicalVolume(solidWorld, G4Material::GetMaterial("Vacuum"), "World");
G4VPhysicalVolume* physicsWorld = new G4PVPlacement(0, G4ThreeVector(), logicWorld, "World", 0, false, 0, true);
// Debug for particle range
// G4Box* solidPlane1 = new G4Box("Plane1", 10 * m, 10 * m, 2 * mm);
// G4LogicalVolume* logicPlane1 = new G4LogicalVolume(solidPlane1, G4Material::GetMaterial("AluminumAlloySeries5"),
// "Plane1"); new G4PVPlacement(0, G4ThreeVector(), logicPlane1, "Plane1", logicWorld, false, 0, true);
// G4Box* solidPlane2 = new G4Box("Plane2", 10 * m, 10 * m, 10 * mm);
// G4LogicalVolume* logicPlane2 = new G4LogicalVolume(solidPlane2, G4Material::GetMaterial("Taparan"), "Plane2");
// new G4PVPlacement(0, G4ThreeVector(0, 0, -6), logicPlane2, "Plane2", logicWorld, false, 0, true);
// G4Box* solidPlane3 = new G4Box("Plane3", 10 * m, 10 * m, 5 * mm);
// G4LogicalVolume* logicPlane3 = new G4LogicalVolume(solidPlane3, G4Material::GetMaterial("AluminumAlloySeries5"),
// "Plane3"); new G4PVPlacement(0, G4ThreeVector(0, 0, -14.5), logicPlane3, "Plane3", logicWorld, false, 0, true);
// 节点舱
ConstructSectionSphere(logicWorld, (5.18 + 0.33) * m);
// 过渡段 1
@ -361,6 +340,87 @@ G4VPhysicalVolume* DetectorConstruction::Construct() {
0.815 * m / 2);
// 大柱段
ConstructSectionBig(logicWorld, 0.815 * m);
// 体模
G4RotationMatrix* rotationMatrix = new G4RotationMatrix();
rotationMatrix->rotateX(90. * deg);
G4Box* solidPhantom = new G4Box("Phantom", 2. * m, 1. * m, 1. * m);
G4LogicalVolume* logicPhantom = new G4LogicalVolume(solidPhantom, G4Material::GetMaterial("Vacuum"), "Phantom");
G4VPhysicalVolume* physicsPhantom = new G4PVPlacement(rotationMatrix, G4ThreeVector(0, 0, -1 * m), logicPhantom,
"Phantom", logicWorld, false, 0, true);
ConstructHumanPhantom(physicsPhantom);
return physicsWorld;
}
void DetectorConstruction::SetBodyPartSensitivity(G4String, G4bool) {
G4cout << "This method is not currently working !!!!" << G4endl;
}
void DetectorConstruction::SetPhantomSex(G4String newSex) {
sex = newSex;
if (sex == "Male") G4cout << ">> Male Phantom will be built." << G4endl;
if (sex == "Female") G4cout << ">> Female Phantom will be built." << G4endl;
if ((sex != "Female") && (sex != "Male")) G4cout << sex << " can not be defined!" << G4endl;
}
void DetectorConstruction::SetPhantomModel(G4String newModel) {
model = newModel;
if (model == "MIRD") G4cout << " >> Phantom " << model << " will be built." << G4endl;
if (model == "MIRDHead") G4cout << " >> Phantom " << model << " will be built." << G4endl;
if ((model != "MIRD") && (model != "MIRDHead")) G4cout << model << " can not be defined!" << G4endl;
}
void DetectorConstruction::ConstructSDandField() {
G4HumanPhantomSD* SD = new G4HumanPhantomSD("SD", "HumanPhantomCollection");
G4SDManager::GetSDMpointer()->AddNewDetector(SD);
SetSensitiveDetector("logicalHead", SD);
SetSensitiveDetector("logicalSkull", SD);
SetSensitiveDetector("logicalBrain", SD);
if (model != "MIRDHead") {
SetSensitiveDetector("logicalTrunk", SD);
SetSensitiveDetector("logicalLeftLeg", SD);
SetSensitiveDetector("logicalRightLeg", SD);
SetSensitiveDetector("logicalLeftArmBone", SD);
SetSensitiveDetector("logicalRightArmBone", SD);
SetSensitiveDetector("logicalLeftLegBone", SD);
SetSensitiveDetector("logicalRightLegBone", SD);
SetSensitiveDetector("logicalUpperSpine", SD);
SetSensitiveDetector("logicalLeftScapula", SD);
SetSensitiveDetector("logicalRightScapula", SD);
SetSensitiveDetector("logicalLeftAdrenal", SD);
SetSensitiveDetector("logicalRightAdrenal", SD);
SetSensitiveDetector("logicalThymus", SD);
SetSensitiveDetector("logicalLeftClavicle", SD);
SetSensitiveDetector("logicalRightClavicle", SD);
SetSensitiveDetector("logicalSmallIntestine", SD);
SetSensitiveDetector("logicalRibCage", SD);
SetSensitiveDetector("logicalMiddleLowerSpine", SD);
SetSensitiveDetector("logicalStomach", SD);
SetSensitiveDetector("logicalUpperLargeIntestine", SD);
SetSensitiveDetector("logicalLowerLargeIntestine", SD);
SetSensitiveDetector("logicalSpleen", SD);
SetSensitiveDetector("logicalPancreas", SD);
SetSensitiveDetector("logicalLeftKidney", SD);
SetSensitiveDetector("logicalRightKidney", SD);
SetSensitiveDetector("logicalUrinaryBladder", SD);
if (sex == "Female") {
SetSensitiveDetector("logicalLeftOvary", SD);
SetSensitiveDetector("logicalRightOvary", SD);
SetSensitiveDetector("logicalUterus", SD);
SetSensitiveDetector("logicalLeftBreast", SD);
SetSensitiveDetector("logicalRightBreast", SD);
}
if (sex == "Male") {
SetSensitiveDetector("logicalMaleGenitalia", SD);
SetSensitiveDetector("logicalLeftTeste", SD);
SetSensitiveDetector("logicalRightTeste", SD);
}
}
}

33
src/G4HumanPhantomHit.cpp Normal file
View File

@ -0,0 +1,33 @@
#include "G4HumanPhantomHit.h"
#include "G4Circle.hh"
#include "G4Colour.hh"
#include "G4UnitsTable.hh"
#include "G4VVisManager.hh"
#include "G4VisAttributes.hh"
// THIS IS NECESSARY FOR MT MODE
G4ThreadLocal G4Allocator<G4HumanPhantomHit>* G4HumanPhantomHitAllocator = 0;
G4HumanPhantomHit::G4HumanPhantomHit() {}
G4HumanPhantomHit::~G4HumanPhantomHit() {}
G4HumanPhantomHit::G4HumanPhantomHit(const G4HumanPhantomHit& right) : G4VHit() {
bodyPartID = right.bodyPartID;
edep = right.edep;
}
const G4HumanPhantomHit& G4HumanPhantomHit::operator=(const G4HumanPhantomHit& right) {
bodyPartID = right.bodyPartID;
edep = right.edep;
return *this;
}
G4bool G4HumanPhantomHit::operator==(const G4HumanPhantomHit& right) const { return (this == &right) ? true : false; }
void G4HumanPhantomHit::Draw() {}
void G4HumanPhantomHit::Print() {
G4cout << "Energy deposit: " << G4BestUnit(edep, "Energy") << "BodyPartID: " << bodyPartID << G4endl;
}

View File

@ -0,0 +1,85 @@
#include "DetectorConstruction.h"
#include "G4HumanPhantomMessenger.h"
#include "G4RunManager.hh"
#include "G4UIcmdWithAString.hh"
#include "G4UIcmdWithoutParameter.hh"
#include "G4UIdirectory.hh"
#include "globals.hh"
G4HumanPhantomMessenger::G4HumanPhantomMessenger(DetectorConstruction* myUsrPhtm)
: myUserPhantom(myUsrPhtm), bps(false) {
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);
modelCmd->SetDefaultValue("MIRD");
modelCmd->SetCandidates("MIRD MIRDHead");
modelCmd->AvailableForStates(G4State_PreInit, G4State_Idle);
sexCmd = new G4UIcmdWithAString("/phantom/setPhantomSex", this);
sexCmd->SetGuidance("Set sex of Phantom: Male or Female.");
sexCmd->SetParameterName("phantomSex", true);
sexCmd->SetDefaultValue("Female");
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);
}
G4HumanPhantomMessenger::~G4HumanPhantomMessenger() {
delete modelCmd;
delete sexCmd;
delete bodypartCmd;
delete endCmd;
delete phantomDir;
delete bpDir;
}
void G4HumanPhantomMessenger::SetNewValue(G4UIcommand* command, G4String newValue) {
if (command == modelCmd) {
myUserPhantom->SetPhantomModel(newValue);
}
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);
}

46
src/G4HumanPhantomSD.cpp Normal file
View File

@ -0,0 +1,46 @@
#include "G4HumanPhantomSD.h"
#include "G4HCofThisEvent.hh"
#include "G4SDManager.hh"
#include "G4Step.hh"
#include "G4ios.hh"
G4HumanPhantomSD::G4HumanPhantomSD(const G4String& name, const G4String& hitsCollectionName)
: G4VSensitiveDetector(name) {
collectionName.insert(hitsCollectionName);
}
G4HumanPhantomSD::~G4HumanPhantomSD() {}
void G4HumanPhantomSD::Initialize(G4HCofThisEvent* HCE) {
collection = new G4HumanPhantomHitsCollection(SensitiveDetectorName, collectionName[0]);
static G4int HCID = -1;
if (HCID < 0) {
HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
}
HCE->AddHitsCollection(HCID, collection);
}
G4bool G4HumanPhantomSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
G4double edep = aStep->GetTotalEnergyDeposit();
if (edep == 0.) return false;
G4String bodypartName = aStep->GetPreStepPoint()->GetTouchable()->GetVolume()->GetLogicalVolume()->GetName();
// G4cout <<bodypartName <<":" << edep/MeV<< G4endl;
G4HumanPhantomHit* newHit = new G4HumanPhantomHit();
newHit->SetEdep(edep);
newHit->SetBodyPartID(bodypartName);
collection->insert(newHit);
return true;
}
void G4HumanPhantomSD::EndOfEvent(G4HCofThisEvent*) {
// G4int NbHits = collection->entries();
// G4cout << "\n-------->Hits Collection: in this event they are " << NbHits
// << " hits in the tracker chambers: " << G4endl;
// for (G4int i=0;i<NbHits;i++) (*collection)[i]->Print();
}

View File

@ -1,7 +1,6 @@
#include "G4VoxelBreastFactory.hh"
#include "G4VoxelLeftBreast.hh"
#include "G4VoxelRightBreast.hh"
#include "G4VoxelBreastFactory.h"
#include "G4VoxelLeftBreast.h"
#include "G4VoxelRightBreast.h"
G4VoxelBreastFactory::G4VoxelBreastFactory() {
// Map with name of the organ and pointer to the MIRDOrgan class

View File

@ -33,6 +33,7 @@
# 设置颜色
/vis/geometry/set/visibility World 0 false
/vis/geometry/set/visibility Phantom 0 false
/vis/geometry/set/visibility BigWorld 0 false
/vis/geometry/set/colour BigMid 0 0.631373 0.686275 0.733333 0.584314
/vis/geometry/set/colour BigAir 0 0.631373 0.686275 0.733333 0.584314