| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file electromagnetic/TestEm7/src/RunActio
27 /// \brief Implementation of the RunAction cla
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 //....oooOO0OOooo........oooOO0OOooo........oo
32
33 #include "RunAction.hh"
34
35 #include "DetectorConstruction.hh"
36 #include "PhysicsList.hh"
37 #include "PrimaryGeneratorAction.hh"
38 #include "StepMax.hh"
39
40 #include "G4Run.hh"
41 #include "G4RunManager.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4UnitsTable.hh"
44 #include "G4ios.hh"
45 #include "Randomize.hh"
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48
49 RunAction::RunAction(DetectorConstruction* det
50 : G4UserRunAction(),
51 fAnalysisManager(0),
52 fDetector(det),
53 fPhysics(phys),
54 fKinematic(kin),
55 fTallyEdep(new G4double[kMaxTally]),
56 fProjRange(0.),
57 fProjRange2(0.),
58 fEdeptot(0.),
59 fEniel(0.),
60 fNbPrimarySteps(0),
61 fRange(0)
62 {
63 // Book predefined histograms
64 BookHisto();
65 }
66
67 //....oooOO0OOooo........oooOO0OOooo........oo
68
69 RunAction::~RunAction()
70 {
71 delete[] fTallyEdep;
72 }
73
74 //....oooOO0OOooo........oooOO0OOooo........oo
75
76 void RunAction::BeginOfRunAction(const G4Run*
77 {
78 G4cout << "### Run " << aRun->GetRunID() <<
79
80 if (!fAnalysisManager) {
81 BookHisto();
82 }
83
84 CLHEP::HepRandom::showEngineStatus();
85
86 // initialize projected range, tallies, Ebea
87 //
88 fNbPrimarySteps = 0;
89 fRange = 0;
90 fProjRange = fProjRange2 = 0.;
91 fEdeptot = fEniel = 0.;
92 for (G4int j = 0; j < kMaxTally; ++j) {
93 fTallyEdep[j] = 0.;
94 }
95 fKinematic->ResetEbeamCumul();
96
97 if (fAnalysisManager->IsActive()) {
98 fAnalysisManager->OpenFile();
99
100 // histogram "1" is defined by the length
101 // zoomed histograms are defined by UI com
102 G4double length = fDetector->GetAbsorSizeX
103 G4double stepMax = fPhysics->GetStepMaxPro
104 G4int nbBins = 100;
105 if (stepMax < DBL_MAX) {
106 G4int nb = (G4int)(0.5 + length / stepMa
107 nbBins = std::min(std::max(nbBins, nb),
108 }
109 fAnalysisManager->SetH1(1, nbBins, 0., len
110 }
111 }
112
113 //....oooOO0OOooo........oooOO0OOooo........oo
114
115 void RunAction::EndOfRunAction(const G4Run* aR
116 {
117 G4int nbofEvents = aRun->GetNumberOfEvent();
118 if (nbofEvents == 0) return;
119
120 // run conditions
121 //
122 const G4Material* material = fDetector->GetA
123 G4double density = material->GetDensity();
124
125 G4String particle = fKinematic->GetParticleG
126 G4double energy = fKinematic->GetParticleGun
127 G4cout << "\n The run consists of " << nbofE
128 << G4BestUnit(energy, "Energy") << "
129 << G4BestUnit(fDetector->GetAbsorSize
130 << " (density: " << G4BestUnit(densit
131
132 // compute projected range and straggling
133 //
134 if (fRange > 0) {
135 fProjRange /= fRange;
136 fProjRange2 /= fRange;
137 }
138 G4double rms = fProjRange2 - fProjRange * fP
139 if (rms > 0.)
140 rms = std::sqrt(rms);
141 else
142 rms = 0.;
143
144 G4double nstep = G4double(fNbPrimarySteps) /
145
146 G4cout.precision(6);
147 G4cout << "\n Projected Range= " << G4BestUn
148 << " rms= " << G4BestUnit(rms, "Len
149 G4cout << " Mean number of primary steps = "
150
151 // compute energy deposition and niel
152 //
153 fEdeptot /= nbofEvents;
154 G4cout << " Total energy deposit= " << G4Bes
155 fEniel /= nbofEvents;
156 G4cout << " niel energy deposit = " << G4Bes
157
158 // print dose in tallies
159 //
160 G4int tallyNumber = fDetector->GetTallyNumbe
161 if (tallyNumber > 0) {
162 G4double Ebeam = fKinematic->GetEbeamCumul
163 G4cout << "\n-----------------------------
164 G4cout << " Cumulated Doses : \tEdep
165 for (G4int j = 0; j < tallyNumber; ++j) {
166 G4double Edep = fTallyEdep[j], ratio = 1
167 G4double tallyMass = fDetector->GetTally
168 G4double Dose = Edep / tallyMass;
169 G4cout << " tally " << j << ": \t \t" <<
170 << " % \t" << G4BestUnit(Dose, "D
171 }
172 G4cout << "\n-----------------------------
173 G4cout << G4endl;
174 }
175
176 if (fAnalysisManager->IsActive()) {
177 // normalize histograms
178 //
179 for (G4int j = 1; j < 3; ++j) {
180 G4double binWidth = fAnalysisManager->Ge
181 G4double fac = (mm / MeV) / (nbofEvents
182 fAnalysisManager->ScaleH1(j, fac);
183 }
184 fAnalysisManager->ScaleH1(3, 1. / nbofEven
185
186 // save histograms
187 fAnalysisManager->Write();
188 fAnalysisManager->CloseFile();
189 }
190
191 // show Rndm status
192 //
193 CLHEP::HepRandom::showEngineStatus();
194 }
195
196 //....oooOO0OOooo........oooOO0OOooo........oo
197
198 void RunAction::BookHisto()
199 {
200 // Create or get analysis manager
201 // The choice of analysis technology is done
202 // in HistoManager.hh
203 fAnalysisManager = G4AnalysisManager::Instan
204 fAnalysisManager->SetDefaultFileType("root")
205 fAnalysisManager->SetFileName("testem7");
206 fAnalysisManager->SetVerboseLevel(1);
207 fAnalysisManager->SetActivation(true); // e
208
209 // Define histograms start values
210 const G4int kMaxHisto = 4;
211 const G4String id[] = {"h0", "h1", "h2", "h3
212 const G4String title[] = {
213 "dummy", // 0
214 "Edep (MeV/mm) along absorber ", // 1
215 "Edep (MeV/mm) along absorber zoomed", //
216 "projectile range" // 3
217 };
218
219 // Default values (to be reset via /analysis
220 G4int nbins = 100;
221 G4double vmin = 0.;
222 G4double vmax = 100.;
223
224 // Create all histograms as inactivated
225 // as we have not yet set nbins, vmin, vmax
226 for (G4int k = 0; k < kMaxHisto; ++k) {
227 G4int ih = fAnalysisManager->CreateH1(id[k
228 G4bool activ = false;
229 if (k == 1) activ = true;
230 fAnalysisManager->SetH1Activation(ih, acti
231 }
232 }
233
234 //....oooOO0OOooo........oooOO0OOooo........oo
235