| 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 medical/GammaTherapy/src/Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 //....oooOO0OOooo........oooOO0OOooo........oo
32
33 #include "Run.hh"
34
35 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh"
37 #include "PrimaryGeneratorAction.hh"
38
39 #include "G4EmCalculator.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4Track.hh"
42 #include "G4UnitsTable.hh"
43 #include "G4VPhysicalVolume.hh"
44
45 #include <iomanip>
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48
49 Run::Run(DetectorConstruction* det, HistoManag
50 : G4Run(), fDetector(det), fHistoMgr(histoMg
51 {
52 fSumR = 0;
53 fNevt = fNelec = fNposit = fNgam = fNstep =
54 0;
55
56 fAnalysisManager = G4AnalysisManager::Instan
57
58 fHistoId = fHistoMgr->GetHistoIdentifiers();
59 fNHisto = fHistoId.size();
60
61 fCheckVolume = fDetector->GetCheckVolume();
62 fGasVolume = fDetector->GetGasVolume();
63 fPhantom = fDetector->GetPhantom();
64 fTarget1 = fDetector->GetTarget1();
65 fTarget2 = fDetector->GetTarget2();
66
67 fNBinsR = fDetector->GetNumberDivR();
68 fNBinsZ = fDetector->GetNumberDivZ();
69
70 fScoreZ = fDetector->GetScoreZ();
71 fAbsorberZ = fDetector->GetAbsorberZ();
72 fAbsorberR = fDetector->GetAbsorberR();
73 fMaxEnergy = fDetector->GetMaxEnergy();
74
75 fNBinsE = fDetector->GetNumberDivE();
76 fMaxEnergy = fDetector->GetMaxEnergy();
77
78 fStepZ = fAbsorberZ / (G4double)fNBinsZ;
79 fStepR = fAbsorberR / (G4double)fNBinsR;
80 fStepE = fMaxEnergy / (G4double)fNBinsE;
81 fScoreBin = (G4int)(fScoreZ / fStepZ + 0.5);
82
83 fVerbose = fDetector->GetVerbose();
84
85 fGamma = G4Gamma::Gamma();
86 fElectron = G4Electron::Electron();
87 fPositron = G4Positron::Positron();
88
89 G4cout << " " << fNBinsR << " bins R ste
90 G4cout << " " << fNBinsZ << " bins Z ste
91 G4cout << " " << fNBinsE << " bins E ste
92 G4cout << " " << fScoreBin << "-th bin in
93
94 fVolumeR.clear();
95 fEdep.clear();
96 fGammaE.clear();
97
98 fVolumeR.resize(fNBinsR, 0.0);
99 fEdep.resize(fNBinsR, 0.0);
100 fGammaE.resize(fNBinsE, 0.0);
101
102 G4double r1 = 0.0;
103 G4double r2 = fStepR;
104 for (G4int i = 0; i < fNBinsR; ++i) {
105 fVolumeR[i] = cm * cm / (CLHEP::pi * (r2 *
106 r1 = r2;
107 r2 += fStepR;
108 }
109
110 if (fAnalysisManager->GetFileName() != "") f
111 }
112
113 //....oooOO0OOooo........oooOO0OOooo........oo
114
115 Run::~Run() {}
116
117 //....oooOO0OOooo........oooOO0OOooo........oo
118
119 void Run::Merge(const G4Run* run)
120 {
121 const Run* localRun = static_cast<const Run*
122
123 fNevt += localRun->fNevt;
124
125 fNelec += localRun->fNelec;
126 fNgam += localRun->fNgam;
127 fNposit += localRun->fNposit;
128
129 fNgamTar += localRun->fNgamTar;
130 fNgamPh += localRun->fNgamPh;
131 fNeTar += localRun->fNeTar;
132 fNePh += localRun->fNePh;
133
134 fNstep += localRun->fNstep;
135 fNstepTarget += localRun->fNstepTarget;
136
137 fSumR += localRun->fSumR;
138
139 for (int i = 0; i < (int)fEdep.size(); i++)
140 fEdep[i] += localRun->fEdep[i];
141 for (int i = 0; i < (int)fGammaE.size(); i++
142 fGammaE[i] += localRun->fGammaE[i];
143
144 G4Run::Merge(run);
145 }
146
147 //....oooOO0OOooo........oooOO0OOooo........oo
148
149 void Run::EndOfRun()
150 {
151 G4cout << "Histo: End of run actions are sta
152
153 // average
154 G4cout << "=================================
155 G4double x = (G4double)fNevt;
156 if (fNevt > 0) {
157 x = 1.0 / x;
158 }
159 G4double xe = x * (G4double)fNelec;
160 G4double xg = x * (G4double)fNgam;
161 G4double xp = x * (G4double)fNposit;
162 G4double xs = x * (G4double)fNstep;
163 G4double xph = x * (G4double)fNgamPh;
164 G4double xes = x * (G4double)fNstepTarget;
165 G4double xgt = x * (G4double)fNgamTar;
166 G4double xet = x * (G4double)fNeTar;
167 G4double xphe = x * (G4double)fNePh;
168
169 G4cout << "Number of events
170 << G4endl;
171 G4cout << std::setprecision(4) << "Average n
172 G4cout << std::setprecision(4) << "Average n
173 G4cout << std::setprecision(4) << "Average n
174 G4cout << std::setprecision(4) << "Average n
175 G4cout << std::setprecision(4) << "Average n
176 << G4endl;
177 G4cout << std::setprecision(4) << "Average n
178 << G4endl;
179 G4cout << std::setprecision(4) << "Average n
180 << G4endl;
181 G4cout << std::setprecision(4) << "Average n
182 << G4endl;
183 G4cout << std::setprecision(4) << "Average n
184 << G4endl;
185 G4cout << std::setprecision(4) << "Total fGa
186 << x * fSumR / MeV << " MeV " << G4en
187 G4cout << "=================================
188 G4cout << G4endl;
189 G4cout << G4endl;
190
191 G4double sum = 0.0;
192 for (G4int i = 0; i < fNBinsR; i++) {
193 fEdep[i] *= x;
194 sum += fEdep[i];
195 }
196
197 if (fAnalysisManager) {
198 if (fAnalysisManager->IsActive()) {
199 // normalise histograms
200 for (G4int i = 0; i < fNHisto; i++) {
201 fAnalysisManager->GetH1(fHistoId[i])->
202 }
203 G4double nr = fEdep[0];
204 if (nr > 0.0) {
205 nr = 1. / nr;
206 }
207 fAnalysisManager->GetH1(fHistoId[0])->sc
208
209 nr = sum * fStepR;
210 if (nr > 0.0) {
211 nr = 1. / nr;
212 }
213 fAnalysisManager->GetH1(fHistoId[1])->sc
214
215 fAnalysisManager->GetH1(fHistoId[3])
216 ->scale(1000.0 * cm3 / (CLHEP::pi * fA
217 fAnalysisManager->GetH1(fHistoId[4])->sc
218
219 // Write histogram file
220 if (!fAnalysisManager->Write()) {
221 G4Exception("Histo::Save()", "hist01",
222 }
223 G4cout << "### Histo::Save: Histograms a
224 if (fAnalysisManager->CloseFile() && fVe
225 G4cout << " File is cl
226 }
227 }
228 }
229 }
230
231 //....oooOO0OOooo........oooOO0OOooo........oo
232
233 void Run::AddTargetPhoton(const G4DynamicParti
234 {
235 ++fNgamTar;
236 if (fAnalysisManager) {
237 fAnalysisManager->FillH1(fHistoId[5], p->G
238 }
239 }
240
241 //....oooOO0OOooo........oooOO0OOooo........oo
242
243 void Run::AddPhantomPhoton(const G4DynamicPart
244 {
245 ++fNgamPh;
246 }
247
248 //....oooOO0OOooo........oooOO0OOooo........oo
249
250 void Run::AddTargetElectron(const G4DynamicPar
251 {
252 ++fNeTar;
253 if (fAnalysisManager) {
254 fAnalysisManager->FillH1(fHistoId[8], p->G
255 }
256 }
257
258 //....oooOO0OOooo........oooOO0OOooo........oo
259
260 void Run::AddPhantomElectron(const G4DynamicPa
261 {
262 ++fNePh;
263 if (fAnalysisManager) {
264 fAnalysisManager->FillH1(fHistoId[7], p->G
265 }
266 }
267
268 //....oooOO0OOooo........oooOO0OOooo........oo
269
270 void Run::ScoreNewTrack(const G4Track* aTrack)
271 {
272 // Save primary parameters
273 const G4ParticleDefinition* particle = aTrac
274 G4int pid = aTrack->GetParentID();
275 G4VPhysicalVolume* pv = aTrack->GetVolume();
276 const G4DynamicParticle* dp = aTrack->GetDyn
277
278 // primary particle
279 if (0 == pid) {
280 ++fNevt;
281 if (fVerbose) {
282 G4ThreeVector pos = aTrack->GetVertexPos
283 G4ThreeVector dir = aTrack->GetMomentumD
284 G4cout << "TrackingAction: Primary " <<
285 << " Ekin(MeV)= " << aTrack->GetK
286 << "; dir= " << dir << G4endl;
287 }
288 // secondary electron
289 }
290 else if (0 < pid && particle == fElectron) {
291 if (fVerbose) {
292 G4cout << "TrackingAction: Secondary ele
293 }
294 AddElectron();
295 if (pv == fPhantom) {
296 AddPhantomElectron(dp);
297 }
298 else if (pv == fTarget1 || pv == fTarget2)
299 AddTargetElectron(dp);
300 }
301
302 // secondary positron
303 }
304 else if (0 < pid && particle == fPositron) {
305 if (fVerbose) {
306 G4cout << "TrackingAction: Secondary pos
307 }
308 AddPositron();
309
310 // secondary gamma
311 }
312 else if (0 < pid && particle == fGamma) {
313 if (fVerbose) {
314 G4cout << "TrackingAction: Secondary gam
315 << " E= " << aTrack->GetKineticEn
316 }
317 AddPhoton();
318 if (pv == fPhantom) {
319 AddPhantomPhoton(dp);
320 }
321 else if (pv == fTarget1 || pv == fTarget2)
322 AddTargetPhoton(dp);
323 }
324 }
325 }
326
327 //....oooOO0OOooo........oooOO0OOooo........oo
328
329 void Run::AddPhantomGamma(G4double e, G4double
330 {
331 e /= MeV;
332 fSumR += e;
333 G4int bin = (G4int)(e / fStepE);
334 if (bin >= fNBinsE) {
335 bin = fNBinsE - 1;
336 }
337 fGammaE[bin] += e;
338 G4int bin1 = (G4int)(r / fStepR);
339 if (bin1 >= fNBinsR) {
340 bin1 = fNBinsR - 1;
341 }
342 if (fAnalysisManager) {
343 G4AnalysisManager::Instance()->FillH1(fHis
344 G4AnalysisManager::Instance()->FillH1(fHis
345 }
346 }
347
348 //....oooOO0OOooo........oooOO0OOooo........oo
349
350 void Run::AddPhantomStep(G4double edep, G4doub
351 G4double r0, G4double
352 {
353 ++fNstep;
354 G4int nzbin = (G4int)(z0 / fStepZ);
355 if (fVerbose) {
356 G4cout << "Histo: edep(MeV)= " << edep / M
357 << " z1= " << z1 << " r2= " << r2 <
358 << G4endl;
359 }
360 if (nzbin == fScoreBin) {
361 G4int bin = (G4int)(r0 / fStepR);
362 if (bin >= fNBinsR) {
363 bin = fNBinsR - 1;
364 }
365 double w = edep * fVolumeR[bin];
366 fEdep[bin] += w;
367
368 if (fAnalysisManager) {
369 G4AnalysisManager::Instance()->FillH1(fH
370 G4AnalysisManager::Instance()->FillH1(fH
371 G4AnalysisManager::Instance()->FillH1(fH
372 }
373 }
374 G4int bin1 = (G4int)(z1 / fStepZ);
375 if (bin1 >= fNBinsZ) {
376 bin1 = fNBinsZ - 1;
377 }
378 G4int bin2 = (G4int)(z2 / fStepZ);
379 if (bin2 >= fNBinsZ) {
380 bin2 = fNBinsZ - 1;
381 }
382 if (bin1 == bin2) {
383 if (fAnalysisManager) {
384 G4AnalysisManager::Instance()->FillH1(fH
385 if (r1 < fStepR) {
386 G4double w = edep;
387 if (r2 > fStepR) {
388 w *= (fStepR - r1) / (r2 - r1);
389 }
390 G4AnalysisManager::Instance()->FillH1(
391 }
392 }
393 }
394 else {
395 G4int bin;
396
397 if (bin2 < bin1) {
398 bin = bin2;
399 G4double z = z2;
400 bin2 = bin1;
401 z2 = z1;
402 bin1 = bin;
403 z1 = z;
404 }
405 G4double zz1 = z1;
406 G4double zz2 = (bin1 + 1) * fStepZ;
407 G4double rr1 = r1;
408 G4double dz = z2 - z1;
409 G4double dr = r2 - r1;
410 G4double rr2 = r1 + dr * (zz2 - zz1) / dz;
411 for (bin = bin1; bin <= bin2; bin++) {
412 if (fAnalysisManager) {
413 G4double de = edep * (zz2 - zz1) / dz;
414 G4double zf = (zz1 + zz2) * 0.5;
415 {
416 G4AnalysisManager::Instance()->FillH
417 }
418 if (rr1 < fStepR) {
419 G4double w = de;
420 if (rr2 > fStepR) w *= (fStepR - rr1
421 {
422 G4AnalysisManager::Instance()->Fil
423 }
424 }
425 }
426 zz1 = zz2;
427 zz2 = std::min(z2, zz1 + fStepZ);
428 rr1 = rr2;
429 rr2 = rr1 + dr * (zz2 - zz1) / dz;
430 }
431 }
432 }
433
434 //....oooOO0OOooo........oooOO0OOooo........oo
435