File indexing completed on 2025-02-23 09:20:43
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030 #include "RunAction.hh"
0031
0032 #include "DetectorConstruction.hh"
0033 #include "PrimaryGeneratorAction.hh"
0034
0035 #include "G4AccumulableManager.hh"
0036 #include "G4LogicalVolume.hh"
0037 #include "G4LogicalVolumeStore.hh"
0038 #include "G4Run.hh"
0039 #include "G4RunManager.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4UnitsTable.hh"
0042
0043 namespace B1Con
0044 {
0045
0046
0047
0048 RunAction::RunAction()
0049 {
0050
0051
0052 const G4double milligray = 1.e-3 * gray;
0053 const G4double microgray = 1.e-6 * gray;
0054 const G4double nanogray = 1.e-9 * gray;
0055 const G4double picogray = 1.e-12 * gray;
0056
0057 new G4UnitDefinition("milligray", "milliGy", "Dose", milligray);
0058 new G4UnitDefinition("microgray", "microGy", "Dose", microgray);
0059 new G4UnitDefinition("nanogray", "nanoGy", "Dose", nanogray);
0060 new G4UnitDefinition("picogray", "picoGy", "Dose", picogray);
0061
0062
0063 G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
0064 accumulableManager->RegisterAccumulable(fEdep);
0065 accumulableManager->RegisterAccumulable(fEdep2);
0066 accumulableManager->RegisterAccumulable(&fEdepPerEvent);
0067 }
0068
0069
0070
0071 void RunAction::BeginOfRunAction(const G4Run*)
0072 {
0073
0074 G4RunManager::GetRunManager()->SetRandomNumberStore(false);
0075
0076
0077 G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
0078 accumulableManager->Reset();
0079
0080 if (IsMaster()) {
0081 fDoseTally = new G4ConvergenceTester("DOSE_TALLY");
0082
0083 }
0084 }
0085
0086
0087
0088 void RunAction::EndOfRunAction(const G4Run* run)
0089 {
0090 G4int nofEvents = run->GetNumberOfEvent();
0091 if (nofEvents == 0) return;
0092
0093
0094 G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
0095 accumulableManager->Merge();
0096
0097
0098
0099 G4double edep = fEdep.GetValue();
0100 G4double edep2 = fEdep2.GetValue();
0101
0102 G4double rms = edep2 - edep * edep / nofEvents;
0103 if (rms > 0.)
0104 rms = std::sqrt(rms);
0105 else
0106 rms = 0.;
0107
0108 const auto detConstruction = static_cast<const B1::DetectorConstruction*>(
0109 G4RunManager::GetRunManager()->GetUserDetectorConstruction());
0110 G4double mass = detConstruction->GetScoringVolume()->GetMass();
0111 G4double dose = edep / mass;
0112 G4double rmsDose = rms / mass;
0113
0114
0115
0116
0117 const auto generatorAction = static_cast<const B1::PrimaryGeneratorAction*>(
0118 G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction());
0119 G4String runCondition;
0120 if (generatorAction) {
0121 const G4ParticleGun* particleGun = generatorAction->GetParticleGun();
0122 runCondition += particleGun->GetParticleDefinition()->GetParticleName();
0123 runCondition += " of ";
0124 G4double particleEnergy = particleGun->GetParticleEnergy();
0125 runCondition += G4BestUnit(particleEnergy, "Energy");
0126 }
0127
0128
0129
0130 if (IsMaster()) {
0131 for (auto edepPerEvent : fEdepPerEvent.GetVector()) {
0132 G4double dosePerEvent = edepPerEvent / mass / gray;
0133 fDoseTally->AddScore(dosePerEvent);
0134 }
0135
0136 fDoseTally->ShowResult();
0137 fDoseTally->ShowHistory();
0138 delete fDoseTally;
0139 fDoseTally = nullptr;
0140
0141 G4cout << G4endl << "--------------------End of Global Run-----------------------";
0142 }
0143 else {
0144 G4cout << G4endl << "--------------------End of Local Run------------------------";
0145 }
0146
0147 G4cout << G4endl << " The run consists of " << nofEvents << " " << runCondition << G4endl
0148 << " Cumulated dose per run, in scoring volume : " << G4BestUnit(dose, "Dose")
0149 << " rms = " << G4BestUnit(rmsDose, "Dose") << G4endl
0150 << "------------------------------------------------------------" << G4endl << G4endl;
0151 }
0152
0153
0154
0155 void RunAction::AddEdep(G4double edep)
0156 {
0157 fEdep += edep;
0158 fEdep2 += edep * edep;
0159 fEdepPerEvent.AddValue(edep);
0160 }
0161
0162
0163
0164 }