fix: scroing;add: headers and sources, root path

.vscode/c_cpp_properties.json

@@ -5,7 +5,8 @@
             "includePath": [
                 "${default}",
                 "${workspaceFolder}/include",
-                "D:/Geant4/dist/include/Geant4"
+                "D:/Geant4/dist/include/Geant4",
+                "D:/ROOT/include"
             ],
             "defines": [
                 "_DEBUG",

auto.mac

@@ -12,13 +12,13 @@
 
 # 定义 scoring 网格
 /score/create/boxMesh water_box
-/score/mesh/boxSize 638. 630. 555. mm
+/score/mesh/boxSize 319. 315. 277.5 mm
 /score/mesh/nBin 50 50 40
 
 # 定义需要获取的物理量与过滤器
-/score/quantity/energyDesposit eDep
-/score/quantity/doseDesposit Dose
-# 过滤器，本例需要考虑次级粒子，不需要
+# 过滤器，本例需要考虑次级粒子
+/score/quantity/energyDeposit eDep
+/score/quantity/doseDeposit Dose
 # /score/filter/particle gammaFilter gamma
 /score/close
 

include/ActionInitialization.h

@@ -0,0 +1,15 @@
+#ifndef B0ActionInitialization_h
+#define B0ActionInitialization_h
+
+#include "G4VUserActionInitialization.hh"
+
+class ActionInitialization : public G4VUserActionInitialization {
+public:
+    ActionInitialization();
+    ~ActionInitialization() override;
+
+    void BuildForMaster() const override;
+    void Build() const override;
+};
+
+#endif

include/DetectorConstruction.h

@@ -0,0 +1,22 @@
+#ifndef B0DetectorConstruction_h
+#define B0DetectorConstruction_h
+
+#include "G4VUserDetectorConstruction.hh"
+#include "globals.hh"
+
+class G4VPhysicalVolume;
+class G4LogicalVolume;
+
+class DetectorConstruction : public G4VUserDetectorConstruction {
+public:
+    DetectorConstruction();
+    ~DetectorConstruction() override;
+
+    G4VPhysicalVolume* Construct() override;
+    G4LogicalVolume* GetScoringVolume() const { return fScoringVolume; }
+
+protected:
+    G4LogicalVolume* fScoringVolume = nullptr;
+};
+
+#endif

include/EventAction.h

@@ -0,0 +1,24 @@
+#ifndef B0EventAction_h
+#define B0EventAction_h
+
+#include "G4UserEventAction.hh"
+#include "globals.hh"
+
+class RunAction;
+
+class EventAction : public G4UserEventAction {
+public:
+    EventAction(RunAction* r);
+    ~EventAction() override;
+
+    void BeginOfEventAction(const G4Event*) override;
+    void EndOfEventAction(const G4Event*) override;
+
+    void AddEdep(G4double edep) { fEdep += edep; }
+
+private:
+    RunAction* fRunAction = nullptr;
+    G4double fEdep = 0.;
+};
+
+#endif

include/Geometry.h

@@ -1,11 +0,0 @@
-#ifndef GEOMETRY
-#define GEOMETRY
-
-#include "G4VUserDetectorConstruction.hh"
-
-class DetectorConstruction : public G4VUserDetectorConstruction {
-public:
-    virtual G4VPhysicalVolume* Construct();
-};
-
-#endif

include/ParticleSource.h

@@ -1,19 +0,0 @@
-#ifndef SOURCE
-#define SOURCE
-
-#include "G4ParticleGun.hh"
-#include "G4VUserPrimaryGeneratorAction.hh"
-
-class G4Event;
-
-class GeneratorAction : public G4VUserPrimaryGeneratorAction {
-public:
-    GeneratorAction();
-    ~GeneratorAction();
-    virtual void GeneratePrimaries(G4Event*);
-
-private:
-    G4ParticleGun* particleGun;
-};
-
-#endif

include/PrimaryGeneratorAction.h

@@ -0,0 +1,20 @@
+#ifndef B0PrimaryGeneratorAction_h
+#define B0PrimaryGeneratorAction_h
+
+#include "G4ParticleGun.hh"
+#include "G4VUserPrimaryGeneratorAction.hh"
+
+class G4Event;
+
+class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction {
+public:
+    PrimaryGeneratorAction();
+    ~PrimaryGeneratorAction();
+    virtual void GeneratePrimaries(G4Event*);
+    const G4ParticleGun* GetParticleGun() const { return fParticleGun; }
+
+private:
+    G4ParticleGun* fParticleGun;
+};
+
+#endif

include/RunAction.h

@@ -0,0 +1,25 @@
+#ifndef B0RunAction_h
+#define B0RunAction_h
+
+#include "G4Accumulable.hh"
+#include "G4UserRunAction.hh"
+#include "globals.hh"
+
+class G4Run;
+
+class RunAction : public G4UserRunAction {
+public:
+    RunAction();
+    ~RunAction() override;
+
+    void BeginOfRunAction(const G4Run*) override;
+    void EndOfRunAction(const G4Run*) override;
+
+    void AddEdep(G4double edep);
+
+private:
+    G4double fEdep = 0.;
+    G4double fEdep2 = 0.;
+};
+
+#endif

include/SteppingAction.h

@@ -0,0 +1,22 @@
+#ifndef B0SteppingAction_h
+#define B0SteppingAction_h
+
+#include "G4UserSteppingAction.hh"
+#include "globals.hh"
+
+class G4LogicalVolume;
+class EventAction;
+
+class SteppingAction : public G4UserSteppingAction {
+public:
+    SteppingAction(EventAction* e);
+    ~SteppingAction() override;
+
+    void UserSteppingAction(const G4Step*) override;
+
+private:
+    EventAction* fEventAction = nullptr;
+    G4LogicalVolume* fScoringVolume = nullptr;
+};
+
+#endif

main.cpp

@@ -1,30 +1,25 @@
+#include "ActionInitialization.h"
+#include "DetectorConstruction.h"
 #include "G4MTRunManager.hh"
 #include "G4ScoringManager.hh"
 #include "G4UIExecutive.hh"
 #include "G4UImanager.hh"
-#include "G4VUserActionInitialization.hh"
 #include "G4VisExecutive.hh"
-#include "Geometry.h"
-#include "ParticleSource.h"
+#include "PrimaryGeneratorAction.h"
 #include "QGSP_BERT.hh"
 
-
-class ActionInitialzation : public G4VUserActionInitialization {
-public:
-    virtual void Build() const { SetUserAction(new GeneratorAction); }
-};
-
 int main(int argc, char** argv) {
     G4UIExecutive* ui = nullptr;
     if (argc == 1) ui = new G4UIExecutive(argc, argv);
 
-    auto runManager = new G4MTRunManager;
-    auto visManager = new G4VisExecutive;
-    auto UIManager = G4UImanager::GetUIpointer();
+    G4MTRunManager* runManager = new G4MTRunManager;
+    G4VisExecutive* visManager = new G4VisExecutive;
+    G4UImanager* UIManager = G4UImanager::GetUIpointer();
+    G4ScoringManager::GetScoringManager();
 
     runManager->SetUserInitialization(new DetectorConstruction());
-    runManager->SetUserInitialization(new QGSP_BERT());
-    runManager->SetUserInitialization(new ActionInitialzation());
+    runManager->SetUserInitialization(new QGSP_BERT);
+    runManager->SetUserInitialization(new ActionInitialization());
 
     visManager->Initialize();
     if (!ui) {

src/ActionInitialization.cpp

@@ -0,0 +1,27 @@
+#include "ActionInitialization.h"
+
+#include "EventAction.h"
+#include "PrimaryGeneratorAction.h"
+#include "RunAction.h"
+#include "SteppingAction.h"
+
+ActionInitialization::ActionInitialization() {}
+
+ActionInitialization::~ActionInitialization() {}
+
+void ActionInitialization::BuildForMaster() const {
+    RunAction* runAction = new RunAction;
+    SetUserAction(runAction);
+}
+
+void ActionInitialization::Build() const {
+    SetUserAction(new PrimaryGeneratorAction);
+
+    RunAction* runAction = new RunAction;
+    SetUserAction(runAction);
+
+    EventAction* eventAction = new EventAction(runAction);
+    SetUserAction(eventAction);
+
+    SetUserAction(new SteppingAction(eventAction));
+}

src/DetectorConstruction.cpp

@@ -0,0 +1,22 @@
+#include "DetectorConstruction.h"
+
+#include "G4Box.hh"
+#include "G4LogicalVolume.hh"
+#include "G4NistManager.hh"
+#include "G4PVPlacement.hh"
+#include "G4SystemOfUnits.hh"
+
+DetectorConstruction::DetectorConstruction() {}
+
+DetectorConstruction::~DetectorConstruction() {}
+
+G4VPhysicalVolume* DetectorConstruction::Construct() {
+    G4NistManager* nist = G4NistManager::Instance();
+
+    G4Box* solid_world = new G4Box("world", 638 / 2 * mm, 630 / 2 * mm, 555. / 2 * mm);
+    G4LogicalVolume* logic_world = new G4LogicalVolume(solid_world, nist->FindOrBuildMaterial("G4_WATER"), "world");
+    G4VPhysicalVolume* physics_world =
+        new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logic_world, "world", 0, false, 0, true);
+
+    return physics_world;
+}

src/EventAction.cpp

@@ -0,0 +1,13 @@
+#include "EventAction.h"
+
+#include "G4Event.hh"
+#include "G4RunManager.hh"
+#include "RunAction.h"
+
+EventAction::EventAction(RunAction* r) : fRunAction(r) {}
+
+EventAction::~EventAction() {}
+
+void EventAction::BeginOfEventAction(const G4Event* event) { fEdep = 0.; }
+
+void EventAction::EndOfEventAction(const G4Event* event) { fRunAction->AddEdep(fEdep); }

src/Geometry.cpp

@@ -1,18 +0,0 @@
-#include "Geometry.h"
-
-#include "G4Box.hh"
-#include "G4LogicalVolume.hh"
-#include "G4NistManager.hh"
-#include "G4PVPlacement.hh"
-#include "G4SystemOfUnits.hh"
-
-G4VPhysicalVolume* DetectorConstruction::Construct() {
-    // Solid
-    auto solid_world = new G4Box("world", 638 / 2 * mm, 630 / 2 * mm, 555. / 2 * mm);
-    // Logic
-    auto nist = G4NistManager::Instance();
-    auto logic_world = new G4LogicalVolume(solid_world, nist->FindOrBuildMaterial("G4_WATER"), "world");
-    // Physics
-    auto physics_world = new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logic_world, "world", 0, false, 0, true);
-    return physics_world;
-}

src/ParticleSource.cpp

@@ -1,18 +0,0 @@
-#include "ParticleSource.h"
-
-#include "G4ParticleTable.hh"
-#include "G4SystemOfUnits.hh"
-
-GeneratorAction::GeneratorAction() {
-    auto table = G4ParticleTable::GetParticleTable();
-
-    particleGun = new G4ParticleGun(1);
-    particleGun->SetParticleDefinition(table->FindParticle("gamma"));
-    particleGun->SetParticlePosition(G4ThreeVector(0, 0, 555. / 2 * mm));
-    particleGun->SetParticleMomentumDirection(G4ThreeVector(0, 0, -1));
-    particleGun->SetParticleEnergy(1 * MeV);
-}
-
-GeneratorAction::~GeneratorAction() { delete particleGun; }
-
-void GeneratorAction::GeneratePrimaries(G4Event* e) { particleGun->GeneratePrimaryVertex(e); }

src/PrimaryGeneratorAction.cpp

@@ -0,0 +1,18 @@
+#include "PrimaryGeneratorAction.h"
+
+#include "G4ParticleTable.hh"
+#include "G4SystemOfUnits.hh"
+
+PrimaryGeneratorAction::PrimaryGeneratorAction() {
+    auto table = G4ParticleTable::GetParticleTable();
+
+    fParticleGun = new G4ParticleGun(1);
+    fParticleGun->SetParticleDefinition(table->FindParticle("gamma"));
+    fParticleGun->SetParticlePosition(G4ThreeVector(0, 0, 555. / 2 * mm));
+    fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0, 0, -1));
+    fParticleGun->SetParticleEnergy(1 * MeV);
+}
+
+PrimaryGeneratorAction::~PrimaryGeneratorAction() { delete fParticleGun; }
+
+void PrimaryGeneratorAction::GeneratePrimaries(G4Event* e) { fParticleGun->GeneratePrimaryVertex(e); }

src/RunAction.cpp

@@ -0,0 +1,81 @@
+#include "RunAction.h"
+
+#include "DetectorConstruction.h"
+#include "G4AccumulableManager.hh"
+#include "G4LogicalVolume.hh"
+#include "G4LogicalVolumeStore.hh"
+#include "G4Run.hh"
+#include "G4RunManager.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4UnitsTable.hh"
+#include "PrimaryGeneratorAction.h"
+
+RunAction::RunAction() {
+    const G4double milligray = 1.e-3 * gray;
+    const G4double microgray = 1.e-6 * gray;
+    const G4double nanogray = 1.e-9 * gray;
+    const G4double picogray = 1.e-12 * gray;
+
+    new G4UnitDefinition("milligray", "mGy", "Dose", milligray);
+    new G4UnitDefinition("microgray", "uGy", "Dose", microgray);
+    new G4UnitDefinition("nanogray", "nGy", "Dose", nanogray);
+    new G4UnitDefinition("picogray", "pGy", "Dose", picogray);
+}
+
+RunAction::~RunAction() {}
+
+void RunAction::BeginOfRunAction(const G4Run*) {
+    G4RunManager::GetRunManager()->SetRandomNumberStore(false);
+
+    G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
+    accumulableManager->Reset();
+}
+
+void RunAction::EndOfRunAction(const G4Run* run) {
+    G4int nofEvents = run->GetNumberOfEvent();
+    if (nofEvents == 0) return;
+
+    G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
+    accumulableManager->Merge();
+
+    G4double rms = fEdep2 - fEdep * fEdep / nofEvents;
+    if (rms > 0.)
+        rms = std::sqrt(rms);
+    else
+        rms = 0.;
+
+    const DetectorConstruction* detConstruction =
+        static_cast<const DetectorConstruction*>(G4RunManager::GetRunManager()->GetUserDetectorConstruction());
+    G4double mass = detConstruction->GetScoringVolume()->GetMass();
+    G4double dose = fEdep / mass;
+    G4double rmsDose = rms / mass;
+
+    G4String runCondition;
+    const PrimaryGeneratorAction* generatorAction =
+        static_cast<const PrimaryGeneratorAction*>(G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction());
+    if (generatorAction) {
+        const G4ParticleGun* particleGun = generatorAction->GetParticleGun();
+        runCondition += particleGun->GetParticleDefinition()->GetParticleName();
+        runCondition += " of ";
+        G4double particleEnergy = particleGun->GetParticleEnergy();
+        runCondition += G4BestUnit(particleEnergy, "Energy");
+    }
+
+    // Print
+    if (IsMaster()) {
+        G4cout << G4endl << "--------------------End of Global Run-----------------------";
+    } else {
+        G4cout << G4endl << "--------------------End of Local Run------------------------";
+    }
+
+    G4cout << G4endl << " The run consists of " << nofEvents << " " << runCondition << G4endl
+           << " Cumulated dose per run, in scoring volume : " << G4BestUnit(dose, "Dose")
+           << " rms = " << G4BestUnit(rmsDose, "Dose") << G4endl
+           << "------------------------------------------------------------" << G4endl << G4endl;
+}
+
+void RunAction::AddEdep(G4double edep) {
+    // 累加沉积能量
+    fEdep += edep;
+    fEdep2 += edep * edep;
+}

src/SteppingAction.cpp

@@ -0,0 +1,30 @@
+#include "SteppingAction.h"
+#include "DetectorConstruction.h"
+#include "EventAction.h"
+
+#include "G4Event.hh"
+#include "G4LogicalVolume.hh"
+#include "G4RunManager.hh"
+#include "G4Step.hh"
+
+SteppingAction::SteppingAction(EventAction* e) : fEventAction(e) {}
+
+SteppingAction::~SteppingAction() {}
+
+void SteppingAction::UserSteppingAction(const G4Step* step) {
+    if (!fScoringVolume) {
+        const DetectorConstruction* detConstruction =
+            static_cast<const DetectorConstruction*>(G4RunManager::GetRunManager()->GetUserDetectorConstruction());
+        fScoringVolume = detConstruction->GetScoringVolume();
+    }
+
+    // 获取当前的 Volume
+    G4LogicalVolume* volume = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetLogicalVolume();
+
+    // 检查是否在计算的 Volume （Shape 2）
+    if (volume != fScoringVolume) return;
+
+    // 获取当前 step 的沉积能量
+    G4double edepStep = step->GetTotalEnergyDeposit();
+    fEventAction->AddEdep(edepStep);
+}

vis.mac

@@ -16,7 +16,7 @@
 /vis/viewer/set/autoRefresh false
 /vis/verbose errors
 
-# 绘制集合体
+# 绘制几何体
 /vis/drawVolume
 
 # 绘制轨迹 类型为平滑