| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 /// \file medical/GammaTherapy/src/Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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........oooOO0OOooo........oooOO0OOooo......
48
49 Run::Run(DetectorConstruction* det, HistoManager* histoMgr)
50 : G4Run(), fDetector(det), fHistoMgr(histoMgr)
51 {
52 fSumR = 0;
53 fNevt = fNelec = fNposit = fNgam = fNstep = fNgamPh = fNgamTar = fNeTar = fNePh = fNstepTarget =
54 0;
55
56 fAnalysisManager = G4AnalysisManager::Instance();
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 stepR= " << fStepR / mm << " mm " << G4endl;
90 G4cout << " " << fNBinsZ << " bins Z stepZ= " << fStepZ / mm << " mm " << G4endl;
91 G4cout << " " << fNBinsE << " bins E stepE= " << fStepE / MeV << " MeV " << G4endl;
92 G4cout << " " << fScoreBin << "-th bin in Z is used for R distribution" << G4endl;
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 * r2 - r1 * r1));
106 r1 = r2;
107 r2 += fStepR;
108 }
109
110 if (fAnalysisManager->GetFileName() != "") fHistoMgr->Update(det, true);
111 }
112
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114
115 Run::~Run() {}
116
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
118
119 void Run::Merge(const G4Run* run)
120 {
121 const Run* localRun = static_cast<const Run*>(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........oooOO0OOooo........oooOO0OOooo......
148
149 void Run::EndOfRun()
150 {
151 G4cout << "Histo: End of run actions are started" << G4endl;
152
153 // average
154 G4cout << "========================================================" << G4endl;
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 " << std::setprecision(8) << fNevt
170 << G4endl;
171 G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl;
172 G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl;
173 G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl;
174 G4cout << std::setprecision(4) << "Average number of steps in the phantom " << xs << G4endl;
175 G4cout << std::setprecision(4) << "Average number of e- steps in the target " << xes
176 << G4endl;
177 G4cout << std::setprecision(4) << "Average number of g produced in the target " << xgt
178 << G4endl;
179 G4cout << std::setprecision(4) << "Average number of e- produced in the target " << xet
180 << G4endl;
181 G4cout << std::setprecision(4) << "Average number of g produced in the phantom " << xph
182 << G4endl;
183 G4cout << std::setprecision(4) << "Average number of e- produced in the phantom " << xphe
184 << G4endl;
185 G4cout << std::setprecision(4) << "Total fGamma fluence in front of the phantom "
186 << x * fSumR / MeV << " MeV " << G4endl;
187 G4cout << "========================================================" << G4endl;
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])->scale(x);
202 }
203 G4double nr = fEdep[0];
204 if (nr > 0.0) {
205 nr = 1. / nr;
206 }
207 fAnalysisManager->GetH1(fHistoId[0])->scale(nr);
208
209 nr = sum * fStepR;
210 if (nr > 0.0) {
211 nr = 1. / nr;
212 }
213 fAnalysisManager->GetH1(fHistoId[1])->scale(nr);
214
215 fAnalysisManager->GetH1(fHistoId[3])
216 ->scale(1000.0 * cm3 / (CLHEP::pi * fAbsorberR * fAbsorberR * fStepZ));
217 fAnalysisManager->GetH1(fHistoId[4])->scale(1000.0 * cm3 * fVolumeR[0] / fStepZ);
218
219 // Write histogram file
220 if (!fAnalysisManager->Write()) {
221 G4Exception("Histo::Save()", "hist01", FatalException, "Cannot write ROOT file.");
222 }
223 G4cout << "### Histo::Save: Histograms are saved" << G4endl;
224 if (fAnalysisManager->CloseFile() && fVerbose) {
225 G4cout << " File is closed" << G4endl;
226 }
227 }
228 }
229 }
230
231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
232
233 void Run::AddTargetPhoton(const G4DynamicParticle* p)
234 {
235 ++fNgamTar;
236 if (fAnalysisManager) {
237 fAnalysisManager->FillH1(fHistoId[5], p->GetKineticEnergy() / MeV, 1.0);
238 }
239 }
240
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242
243 void Run::AddPhantomPhoton(const G4DynamicParticle*)
244 {
245 ++fNgamPh;
246 }
247
248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249
250 void Run::AddTargetElectron(const G4DynamicParticle* p)
251 {
252 ++fNeTar;
253 if (fAnalysisManager) {
254 fAnalysisManager->FillH1(fHistoId[8], p->GetKineticEnergy() / MeV, 1.0);
255 }
256 }
257
258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
259
260 void Run::AddPhantomElectron(const G4DynamicParticle* p)
261 {
262 ++fNePh;
263 if (fAnalysisManager) {
264 fAnalysisManager->FillH1(fHistoId[7], p->GetKineticEnergy() / MeV, 1.0);
265 }
266 }
267
268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
269
270 void Run::ScoreNewTrack(const G4Track* aTrack)
271 {
272 // Save primary parameters
273 const G4ParticleDefinition* particle = aTrack->GetParticleDefinition();
274 G4int pid = aTrack->GetParentID();
275 G4VPhysicalVolume* pv = aTrack->GetVolume();
276 const G4DynamicParticle* dp = aTrack->GetDynamicParticle();
277
278 // primary particle
279 if (0 == pid) {
280 ++fNevt;
281 if (fVerbose) {
282 G4ThreeVector pos = aTrack->GetVertexPosition();
283 G4ThreeVector dir = aTrack->GetMomentumDirection();
284 G4cout << "TrackingAction: Primary " << particle->GetParticleName()
285 << " Ekin(MeV)= " << aTrack->GetKineticEnergy() / MeV << "; pos= " << pos
286 << "; dir= " << dir << G4endl;
287 }
288 // secondary electron
289 }
290 else if (0 < pid && particle == fElectron) {
291 if (fVerbose) {
292 G4cout << "TrackingAction: Secondary electron " << G4endl;
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 positron " << G4endl;
307 }
308 AddPositron();
309
310 // secondary gamma
311 }
312 else if (0 < pid && particle == fGamma) {
313 if (fVerbose) {
314 G4cout << "TrackingAction: Secondary gamma; parentID= " << pid
315 << " E= " << aTrack->GetKineticEnergy() << G4endl;
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........oooOO0OOooo........oooOO0OOooo......
328
329 void Run::AddPhantomGamma(G4double e, G4double r)
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(fHistoId[6], e, 1.0);
344 G4AnalysisManager::Instance()->FillH1(fHistoId[9], r, e * fVolumeR[bin1]);
345 }
346 }
347
348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
349
350 void Run::AddPhantomStep(G4double edep, G4double r1, G4double z1, G4double r2, G4double z2,
351 G4double r0, G4double z0)
352 {
353 ++fNstep;
354 G4int nzbin = (G4int)(z0 / fStepZ);
355 if (fVerbose) {
356 G4cout << "Histo: edep(MeV)= " << edep / MeV << " at binz= " << nzbin << " r1= " << r1
357 << " z1= " << z1 << " r2= " << r2 << " z2= " << z2 << " r0= " << r0 << " z0= " << z0
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(fHistoId[0], r0, w);
370 G4AnalysisManager::Instance()->FillH1(fHistoId[1], r0, w);
371 G4AnalysisManager::Instance()->FillH1(fHistoId[2], r0, w);
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(fHistoId[3], z0, edep);
385 if (r1 < fStepR) {
386 G4double w = edep;
387 if (r2 > fStepR) {
388 w *= (fStepR - r1) / (r2 - r1);
389 }
390 G4AnalysisManager::Instance()->FillH1(fHistoId[4], z0, w);
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()->FillH1(fHistoId[3], zf, de);
417 }
418 if (rr1 < fStepR) {
419 G4double w = de;
420 if (rr2 > fStepR) w *= (fStepR - rr1) / (rr2 - rr1);
421 {
422 G4AnalysisManager::Instance()->FillH1(fHistoId[4], zf, w);
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........oooOO0OOooo........oooOO0OOooo......
435