| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /// \file medical/GammaTherapy/src/Run.cc 26 /// \file medical/GammaTherapy/src/Run.cc
27 /// \brief Implementation of the Run class 27 /// \brief Implementation of the Run class
28 // 28 //
29 // 29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 32
33 #include "Run.hh" 33 #include "Run.hh"
34 <<
35 #include "DetectorConstruction.hh" 34 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh" <<
37 #include "PrimaryGeneratorAction.hh" 35 #include "PrimaryGeneratorAction.hh"
>> 36 #include "HistoManager.hh"
>> 37 #include "G4Track.hh"
>> 38 #include "G4VPhysicalVolume.hh"
38 39
39 #include "G4EmCalculator.hh" 40 #include "G4EmCalculator.hh"
40 #include "G4SystemOfUnits.hh" 41 #include "G4SystemOfUnits.hh"
41 #include "G4Track.hh" <<
42 #include "G4UnitsTable.hh" 42 #include "G4UnitsTable.hh"
43 #include "G4VPhysicalVolume.hh" <<
44 43
45 #include <iomanip> 44 #include <iomanip>
46 45
47 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 47
49 Run::Run(DetectorConstruction* det, HistoManag 48 Run::Run(DetectorConstruction* det, HistoManager* histoMgr)
50 : G4Run(), fDetector(det), fHistoMgr(histoMg << 49 : G4Run()
>> 50 , fDetector(det)
>> 51 , fHistoMgr(histoMgr)
51 { 52 {
52 fSumR = 0; 53 fSumR = 0;
53 fNevt = fNelec = fNposit = fNgam = fNstep = << 54 fNevt = fNelec = fNposit = fNgam = fNstep = fNgamPh = fNgamTar = fNeTar =
54 0; << 55 fNePh = fNstepTarget = 0;
55 56
56 fAnalysisManager = G4AnalysisManager::Instan 57 fAnalysisManager = G4AnalysisManager::Instance();
57 58
58 fHistoId = fHistoMgr->GetHistoIdentifiers(); 59 fHistoId = fHistoMgr->GetHistoIdentifiers();
59 fNHisto = fHistoId.size(); << 60 fNHisto = fHistoId.size();
60 61
61 fCheckVolume = fDetector->GetCheckVolume(); 62 fCheckVolume = fDetector->GetCheckVolume();
62 fGasVolume = fDetector->GetGasVolume(); << 63 fGasVolume = fDetector->GetGasVolume();
63 fPhantom = fDetector->GetPhantom(); << 64 fPhantom = fDetector->GetPhantom();
64 fTarget1 = fDetector->GetTarget1(); << 65 fTarget1 = fDetector->GetTarget1();
65 fTarget2 = fDetector->GetTarget2(); << 66 fTarget2 = fDetector->GetTarget2();
66 67
67 fNBinsR = fDetector->GetNumberDivR(); 68 fNBinsR = fDetector->GetNumberDivR();
68 fNBinsZ = fDetector->GetNumberDivZ(); 69 fNBinsZ = fDetector->GetNumberDivZ();
69 70
70 fScoreZ = fDetector->GetScoreZ(); << 71 fScoreZ = fDetector->GetScoreZ();
71 fAbsorberZ = fDetector->GetAbsorberZ(); 72 fAbsorberZ = fDetector->GetAbsorberZ();
72 fAbsorberR = fDetector->GetAbsorberR(); 73 fAbsorberR = fDetector->GetAbsorberR();
73 fMaxEnergy = fDetector->GetMaxEnergy(); 74 fMaxEnergy = fDetector->GetMaxEnergy();
74 75
75 fNBinsE = fDetector->GetNumberDivE(); << 76 fNBinsE = fDetector->GetNumberDivE();
76 fMaxEnergy = fDetector->GetMaxEnergy(); 77 fMaxEnergy = fDetector->GetMaxEnergy();
77 78
78 fStepZ = fAbsorberZ / (G4double)fNBinsZ; << 79 fStepZ = fAbsorberZ / (G4double) fNBinsZ;
79 fStepR = fAbsorberR / (G4double)fNBinsR; << 80 fStepR = fAbsorberR / (G4double) fNBinsR;
80 fStepE = fMaxEnergy / (G4double)fNBinsE; << 81 fStepE = fMaxEnergy / (G4double) fNBinsE;
81 fScoreBin = (G4int)(fScoreZ / fStepZ + 0.5); 82 fScoreBin = (G4int)(fScoreZ / fStepZ + 0.5);
82 83
83 fVerbose = fDetector->GetVerbose(); 84 fVerbose = fDetector->GetVerbose();
84 85
85 fGamma = G4Gamma::Gamma(); << 86 fGamma = G4Gamma::Gamma();
86 fElectron = G4Electron::Electron(); 87 fElectron = G4Electron::Electron();
87 fPositron = G4Positron::Positron(); 88 fPositron = G4Positron::Positron();
88 89
89 G4cout << " " << fNBinsR << " bins R ste << 90 G4cout << " " << fNBinsR << " bins R stepR= " << fStepR / mm << " mm "
90 G4cout << " " << fNBinsZ << " bins Z ste << 91 << G4endl;
91 G4cout << " " << fNBinsE << " bins E ste << 92 G4cout << " " << fNBinsZ << " bins Z stepZ= " << fStepZ / mm << " mm "
92 G4cout << " " << fScoreBin << "-th bin in << 93 << G4endl;
>> 94 G4cout << " " << fNBinsE << " bins E stepE= " << fStepE / MeV << " MeV "
>> 95 << G4endl;
>> 96 G4cout << " " << fScoreBin << "-th bin in Z is used for R distribution"
>> 97 << G4endl;
93 98
94 fVolumeR.clear(); 99 fVolumeR.clear();
95 fEdep.clear(); 100 fEdep.clear();
96 fGammaE.clear(); 101 fGammaE.clear();
97 102
98 fVolumeR.resize(fNBinsR, 0.0); 103 fVolumeR.resize(fNBinsR, 0.0);
99 fEdep.resize(fNBinsR, 0.0); 104 fEdep.resize(fNBinsR, 0.0);
100 fGammaE.resize(fNBinsE, 0.0); 105 fGammaE.resize(fNBinsE, 0.0);
101 106
102 G4double r1 = 0.0; 107 G4double r1 = 0.0;
103 G4double r2 = fStepR; 108 G4double r2 = fStepR;
104 for (G4int i = 0; i < fNBinsR; ++i) { << 109 for(G4int i = 0; i < fNBinsR; ++i)
>> 110 {
105 fVolumeR[i] = cm * cm / (CLHEP::pi * (r2 * 111 fVolumeR[i] = cm * cm / (CLHEP::pi * (r2 * r2 - r1 * r1));
106 r1 = r2; << 112 r1 = r2;
107 r2 += fStepR; 113 r2 += fStepR;
108 } 114 }
109 115
110 if (fAnalysisManager->GetFileName() != "") f << 116 if(fAnalysisManager->GetFileName() != "")
>> 117 fHistoMgr->Update(det, true);
111 } 118 }
112 119
113 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 121
115 Run::~Run() {} 122 Run::~Run() {}
116 123
117 //....oooOO0OOooo........oooOO0OOooo........oo 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
118 125
119 void Run::Merge(const G4Run* run) 126 void Run::Merge(const G4Run* run)
120 { 127 {
121 const Run* localRun = static_cast<const Run* 128 const Run* localRun = static_cast<const Run*>(run);
122 129
123 fNevt += localRun->fNevt; 130 fNevt += localRun->fNevt;
124 131
125 fNelec += localRun->fNelec; 132 fNelec += localRun->fNelec;
126 fNgam += localRun->fNgam; 133 fNgam += localRun->fNgam;
127 fNposit += localRun->fNposit; 134 fNposit += localRun->fNposit;
128 135
129 fNgamTar += localRun->fNgamTar; 136 fNgamTar += localRun->fNgamTar;
130 fNgamPh += localRun->fNgamPh; 137 fNgamPh += localRun->fNgamPh;
131 fNeTar += localRun->fNeTar; 138 fNeTar += localRun->fNeTar;
132 fNePh += localRun->fNePh; 139 fNePh += localRun->fNePh;
133 140
134 fNstep += localRun->fNstep; 141 fNstep += localRun->fNstep;
135 fNstepTarget += localRun->fNstepTarget; 142 fNstepTarget += localRun->fNstepTarget;
136 143
137 fSumR += localRun->fSumR; 144 fSumR += localRun->fSumR;
138 145
139 for (int i = 0; i < (int)fEdep.size(); i++) << 146 for(int i = 0; i < (int) fEdep.size(); i++)
140 fEdep[i] += localRun->fEdep[i]; 147 fEdep[i] += localRun->fEdep[i];
141 for (int i = 0; i < (int)fGammaE.size(); i++ << 148 for(int i = 0; i < (int) fGammaE.size(); i++)
142 fGammaE[i] += localRun->fGammaE[i]; 149 fGammaE[i] += localRun->fGammaE[i];
143 150
144 G4Run::Merge(run); 151 G4Run::Merge(run);
145 } 152 }
146 153
147 //....oooOO0OOooo........oooOO0OOooo........oo 154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
148 155
149 void Run::EndOfRun() 156 void Run::EndOfRun()
150 { 157 {
151 G4cout << "Histo: End of run actions are sta 158 G4cout << "Histo: End of run actions are started" << G4endl;
152 159
153 // average 160 // average
154 G4cout << "================================= << 161 G4cout << "========================================================"
155 G4double x = (G4double)fNevt; << 162 << G4endl;
156 if (fNevt > 0) { << 163 G4double x = (G4double) fNevt;
>> 164 if(fNevt > 0)
>> 165 {
157 x = 1.0 / x; 166 x = 1.0 / x;
158 } 167 }
159 G4double xe = x * (G4double)fNelec; << 168 G4double xe = x * (G4double) fNelec;
160 G4double xg = x * (G4double)fNgam; << 169 G4double xg = x * (G4double) fNgam;
161 G4double xp = x * (G4double)fNposit; << 170 G4double xp = x * (G4double) fNposit;
162 G4double xs = x * (G4double)fNstep; << 171 G4double xs = x * (G4double) fNstep;
163 G4double xph = x * (G4double)fNgamPh; << 172 G4double xph = x * (G4double) fNgamPh;
164 G4double xes = x * (G4double)fNstepTarget; << 173 G4double xes = x * (G4double) fNstepTarget;
165 G4double xgt = x * (G4double)fNgamTar; << 174 G4double xgt = x * (G4double) fNgamTar;
166 G4double xet = x * (G4double)fNeTar; << 175 G4double xet = x * (G4double) fNeTar;
167 G4double xphe = x * (G4double)fNePh; << 176 G4double xphe = x * (G4double) fNePh;
168 << 177
169 G4cout << "Number of events << 178 G4cout << "Number of events "
170 << G4endl; << 179 << std::setprecision(8) << fNevt << G4endl;
171 G4cout << std::setprecision(4) << "Average n << 180 G4cout << std::setprecision(4)
172 G4cout << std::setprecision(4) << "Average n << 181 << "Average number of e- " << xe << G4endl;
173 G4cout << std::setprecision(4) << "Average n << 182 G4cout << std::setprecision(4)
174 G4cout << std::setprecision(4) << "Average n << 183 << "Average number of gamma " << xg << G4endl;
175 G4cout << std::setprecision(4) << "Average n << 184 G4cout << std::setprecision(4)
176 << G4endl; << 185 << "Average number of e+ " << xp << G4endl;
177 G4cout << std::setprecision(4) << "Average n << 186 G4cout << std::setprecision(4)
178 << G4endl; << 187 << "Average number of steps in the phantom " << xs << G4endl;
179 G4cout << std::setprecision(4) << "Average n << 188 G4cout << std::setprecision(4)
180 << G4endl; << 189 << "Average number of e- steps in the target " << xes << G4endl;
181 G4cout << std::setprecision(4) << "Average n << 190 G4cout << std::setprecision(4)
>> 191 << "Average number of g produced in the target " << xgt << G4endl;
>> 192 G4cout << std::setprecision(4)
>> 193 << "Average number of e- produced in the target " << xet << G4endl;
>> 194 G4cout << std::setprecision(4)
>> 195 << "Average number of g produced in the phantom " << xph << G4endl;
>> 196 G4cout << std::setprecision(4)
>> 197 << "Average number of e- produced in the phantom " << xphe << G4endl;
>> 198 G4cout << std::setprecision(4)
>> 199 << "Total fGamma fluence in front of the phantom " << x * fSumR / MeV
>> 200 << " MeV " << G4endl;
>> 201 G4cout << "========================================================"
182 << G4endl; 202 << 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; 203 G4cout << G4endl;
189 G4cout << G4endl; 204 G4cout << G4endl;
190 205
191 G4double sum = 0.0; 206 G4double sum = 0.0;
192 for (G4int i = 0; i < fNBinsR; i++) { << 207 for(G4int i = 0; i < fNBinsR; i++)
>> 208 {
193 fEdep[i] *= x; 209 fEdep[i] *= x;
194 sum += fEdep[i]; 210 sum += fEdep[i];
195 } 211 }
196 212
197 if (fAnalysisManager) { << 213 if(fAnalysisManager)
198 if (fAnalysisManager->IsActive()) { << 214 {
>> 215 if(fAnalysisManager->IsActive())
>> 216 {
199 // normalise histograms 217 // normalise histograms
200 for (G4int i = 0; i < fNHisto; i++) { << 218 for(G4int i = 0; i < fNHisto; i++)
>> 219 {
201 fAnalysisManager->GetH1(fHistoId[i])-> 220 fAnalysisManager->GetH1(fHistoId[i])->scale(x);
202 } 221 }
203 G4double nr = fEdep[0]; 222 G4double nr = fEdep[0];
204 if (nr > 0.0) { << 223 if(nr > 0.0)
>> 224 {
205 nr = 1. / nr; 225 nr = 1. / nr;
206 } 226 }
207 fAnalysisManager->GetH1(fHistoId[0])->sc 227 fAnalysisManager->GetH1(fHistoId[0])->scale(nr);
208 228
209 nr = sum * fStepR; 229 nr = sum * fStepR;
210 if (nr > 0.0) { << 230 if(nr > 0.0)
>> 231 {
211 nr = 1. / nr; 232 nr = 1. / nr;
212 } 233 }
213 fAnalysisManager->GetH1(fHistoId[1])->sc 234 fAnalysisManager->GetH1(fHistoId[1])->scale(nr);
214 235
215 fAnalysisManager->GetH1(fHistoId[3]) 236 fAnalysisManager->GetH1(fHistoId[3])
216 ->scale(1000.0 * cm3 / (CLHEP::pi * fA 237 ->scale(1000.0 * cm3 / (CLHEP::pi * fAbsorberR * fAbsorberR * fStepZ));
217 fAnalysisManager->GetH1(fHistoId[4])->sc << 238 fAnalysisManager->GetH1(fHistoId[4])
>> 239 ->scale(1000.0 * cm3 * fVolumeR[0] / fStepZ);
218 240
219 // Write histogram file 241 // Write histogram file
220 if (!fAnalysisManager->Write()) { << 242 if(!fAnalysisManager->Write())
221 G4Exception("Histo::Save()", "hist01", << 243 {
>> 244 G4Exception("Histo::Save()", "hist01", FatalException,
>> 245 "Cannot write ROOT file.");
222 } 246 }
223 G4cout << "### Histo::Save: Histograms a 247 G4cout << "### Histo::Save: Histograms are saved" << G4endl;
224 if (fAnalysisManager->CloseFile() && fVe << 248 if(fAnalysisManager->CloseFile() && fVerbose)
>> 249 {
225 G4cout << " File is cl 250 G4cout << " File is closed" << G4endl;
226 } 251 }
227 } 252 }
228 } 253 }
229 } 254 }
230 255
231 //....oooOO0OOooo........oooOO0OOooo........oo 256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
232 257
233 void Run::AddTargetPhoton(const G4DynamicParti 258 void Run::AddTargetPhoton(const G4DynamicParticle* p)
234 { 259 {
235 ++fNgamTar; 260 ++fNgamTar;
236 if (fAnalysisManager) { << 261 if(fAnalysisManager)
>> 262 {
237 fAnalysisManager->FillH1(fHistoId[5], p->G 263 fAnalysisManager->FillH1(fHistoId[5], p->GetKineticEnergy() / MeV, 1.0);
238 } 264 }
239 } 265 }
240 266
241 //....oooOO0OOooo........oooOO0OOooo........oo 267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242 268
243 void Run::AddPhantomPhoton(const G4DynamicPart << 269 void Run::AddPhantomPhoton(const G4DynamicParticle*) { ++fNgamPh; }
244 { <<
245 ++fNgamPh; <<
246 } <<
247 270
248 //....oooOO0OOooo........oooOO0OOooo........oo 271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249 272
250 void Run::AddTargetElectron(const G4DynamicPar 273 void Run::AddTargetElectron(const G4DynamicParticle* p)
251 { 274 {
252 ++fNeTar; 275 ++fNeTar;
253 if (fAnalysisManager) { << 276 if(fAnalysisManager)
>> 277 {
254 fAnalysisManager->FillH1(fHistoId[8], p->G 278 fAnalysisManager->FillH1(fHistoId[8], p->GetKineticEnergy() / MeV, 1.0);
255 } 279 }
256 } 280 }
257 281
258 //....oooOO0OOooo........oooOO0OOooo........oo 282 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
259 283
260 void Run::AddPhantomElectron(const G4DynamicPa 284 void Run::AddPhantomElectron(const G4DynamicParticle* p)
261 { 285 {
262 ++fNePh; 286 ++fNePh;
263 if (fAnalysisManager) { << 287 if(fAnalysisManager)
>> 288 {
264 fAnalysisManager->FillH1(fHistoId[7], p->G 289 fAnalysisManager->FillH1(fHistoId[7], p->GetKineticEnergy() / MeV, 1.0);
265 } 290 }
266 } 291 }
267 292
268 //....oooOO0OOooo........oooOO0OOooo........oo 293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
269 294
270 void Run::ScoreNewTrack(const G4Track* aTrack) 295 void Run::ScoreNewTrack(const G4Track* aTrack)
271 { 296 {
272 // Save primary parameters 297 // Save primary parameters
273 const G4ParticleDefinition* particle = aTrac 298 const G4ParticleDefinition* particle = aTrack->GetParticleDefinition();
274 G4int pid = aTrack->GetParentID(); << 299 G4int pid = aTrack->GetParentID();
275 G4VPhysicalVolume* pv = aTrack->GetVolume(); << 300 G4VPhysicalVolume* pv = aTrack->GetVolume();
276 const G4DynamicParticle* dp = aTrack->GetDyn << 301 const G4DynamicParticle* dp = aTrack->GetDynamicParticle();
277 302
278 // primary particle 303 // primary particle
279 if (0 == pid) { << 304 if(0 == pid)
>> 305 {
280 ++fNevt; 306 ++fNevt;
281 if (fVerbose) { << 307 if(fVerbose)
>> 308 {
282 G4ThreeVector pos = aTrack->GetVertexPos 309 G4ThreeVector pos = aTrack->GetVertexPosition();
283 G4ThreeVector dir = aTrack->GetMomentumD 310 G4ThreeVector dir = aTrack->GetMomentumDirection();
284 G4cout << "TrackingAction: Primary " << 311 G4cout << "TrackingAction: Primary " << particle->GetParticleName()
285 << " Ekin(MeV)= " << aTrack->GetK << 312 << " Ekin(MeV)= " << aTrack->GetKineticEnergy() / MeV
286 << "; dir= " << dir << G4endl; << 313 << "; pos= " << pos << "; dir= " << dir << G4endl;
287 } 314 }
288 // secondary electron 315 // secondary electron
289 } 316 }
290 else if (0 < pid && particle == fElectron) { << 317 else if(0 < pid && particle == fElectron)
291 if (fVerbose) { << 318 {
>> 319 if(fVerbose)
>> 320 {
292 G4cout << "TrackingAction: Secondary ele 321 G4cout << "TrackingAction: Secondary electron " << G4endl;
293 } 322 }
294 AddElectron(); 323 AddElectron();
295 if (pv == fPhantom) { << 324 if(pv == fPhantom)
>> 325 {
296 AddPhantomElectron(dp); 326 AddPhantomElectron(dp);
297 } 327 }
298 else if (pv == fTarget1 || pv == fTarget2) << 328 else if(pv == fTarget1 || pv == fTarget2)
>> 329 {
299 AddTargetElectron(dp); 330 AddTargetElectron(dp);
300 } 331 }
301 332
302 // secondary positron 333 // secondary positron
303 } 334 }
304 else if (0 < pid && particle == fPositron) { << 335 else if(0 < pid && particle == fPositron)
305 if (fVerbose) { << 336 {
>> 337 if(fVerbose)
>> 338 {
306 G4cout << "TrackingAction: Secondary pos 339 G4cout << "TrackingAction: Secondary positron " << G4endl;
307 } 340 }
308 AddPositron(); 341 AddPositron();
309 342
310 // secondary gamma 343 // secondary gamma
311 } 344 }
312 else if (0 < pid && particle == fGamma) { << 345 else if(0 < pid && particle == fGamma)
313 if (fVerbose) { << 346 {
>> 347 if(fVerbose)
>> 348 {
314 G4cout << "TrackingAction: Secondary gam 349 G4cout << "TrackingAction: Secondary gamma; parentID= " << pid
315 << " E= " << aTrack->GetKineticEn 350 << " E= " << aTrack->GetKineticEnergy() << G4endl;
316 } 351 }
317 AddPhoton(); 352 AddPhoton();
318 if (pv == fPhantom) { << 353 if(pv == fPhantom)
>> 354 {
319 AddPhantomPhoton(dp); 355 AddPhantomPhoton(dp);
320 } 356 }
321 else if (pv == fTarget1 || pv == fTarget2) << 357 else if(pv == fTarget1 || pv == fTarget2)
>> 358 {
322 AddTargetPhoton(dp); 359 AddTargetPhoton(dp);
323 } 360 }
324 } 361 }
325 } 362 }
326 363
327 //....oooOO0OOooo........oooOO0OOooo........oo 364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
328 365
329 void Run::AddPhantomGamma(G4double e, G4double 366 void Run::AddPhantomGamma(G4double e, G4double r)
330 { 367 {
331 e /= MeV; 368 e /= MeV;
332 fSumR += e; 369 fSumR += e;
333 G4int bin = (G4int)(e / fStepE); 370 G4int bin = (G4int)(e / fStepE);
334 if (bin >= fNBinsE) { << 371 if(bin >= fNBinsE)
>> 372 {
335 bin = fNBinsE - 1; 373 bin = fNBinsE - 1;
336 } 374 }
337 fGammaE[bin] += e; 375 fGammaE[bin] += e;
338 G4int bin1 = (G4int)(r / fStepR); 376 G4int bin1 = (G4int)(r / fStepR);
339 if (bin1 >= fNBinsR) { << 377 if(bin1 >= fNBinsR)
>> 378 {
340 bin1 = fNBinsR - 1; 379 bin1 = fNBinsR - 1;
341 } 380 }
342 if (fAnalysisManager) { << 381 if(fAnalysisManager)
>> 382 {
343 G4AnalysisManager::Instance()->FillH1(fHis 383 G4AnalysisManager::Instance()->FillH1(fHistoId[6], e, 1.0);
344 G4AnalysisManager::Instance()->FillH1(fHis 384 G4AnalysisManager::Instance()->FillH1(fHistoId[9], r, e * fVolumeR[bin1]);
345 } 385 }
346 } 386 }
347 387
348 //....oooOO0OOooo........oooOO0OOooo........oo 388 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
349 389
350 void Run::AddPhantomStep(G4double edep, G4doub << 390 void Run::AddPhantomStep(G4double edep, G4double r1, G4double z1, G4double r2,
351 G4double r0, G4double << 391 G4double z2, G4double r0, G4double z0)
352 { 392 {
353 ++fNstep; 393 ++fNstep;
354 G4int nzbin = (G4int)(z0 / fStepZ); 394 G4int nzbin = (G4int)(z0 / fStepZ);
355 if (fVerbose) { << 395 if(fVerbose)
356 G4cout << "Histo: edep(MeV)= " << edep / M << 396 {
357 << " z1= " << z1 << " r2= " << r2 < << 397 G4cout << "Histo: edep(MeV)= " << edep / MeV << " at binz= " << nzbin
358 << G4endl; << 398 << " r1= " << r1 << " z1= " << z1 << " r2= " << r2 << " z2= " << z2
>> 399 << " r0= " << r0 << " z0= " << z0 << G4endl;
359 } 400 }
360 if (nzbin == fScoreBin) { << 401 if(nzbin == fScoreBin)
>> 402 {
361 G4int bin = (G4int)(r0 / fStepR); 403 G4int bin = (G4int)(r0 / fStepR);
362 if (bin >= fNBinsR) { << 404 if(bin >= fNBinsR)
>> 405 {
363 bin = fNBinsR - 1; 406 bin = fNBinsR - 1;
364 } 407 }
365 double w = edep * fVolumeR[bin]; 408 double w = edep * fVolumeR[bin];
366 fEdep[bin] += w; 409 fEdep[bin] += w;
367 410
368 if (fAnalysisManager) { << 411 if(fAnalysisManager)
>> 412 {
369 G4AnalysisManager::Instance()->FillH1(fH 413 G4AnalysisManager::Instance()->FillH1(fHistoId[0], r0, w);
370 G4AnalysisManager::Instance()->FillH1(fH 414 G4AnalysisManager::Instance()->FillH1(fHistoId[1], r0, w);
371 G4AnalysisManager::Instance()->FillH1(fH 415 G4AnalysisManager::Instance()->FillH1(fHistoId[2], r0, w);
372 } 416 }
373 } 417 }
374 G4int bin1 = (G4int)(z1 / fStepZ); 418 G4int bin1 = (G4int)(z1 / fStepZ);
375 if (bin1 >= fNBinsZ) { << 419 if(bin1 >= fNBinsZ)
>> 420 {
376 bin1 = fNBinsZ - 1; 421 bin1 = fNBinsZ - 1;
377 } 422 }
378 G4int bin2 = (G4int)(z2 / fStepZ); 423 G4int bin2 = (G4int)(z2 / fStepZ);
379 if (bin2 >= fNBinsZ) { << 424 if(bin2 >= fNBinsZ)
>> 425 {
380 bin2 = fNBinsZ - 1; 426 bin2 = fNBinsZ - 1;
381 } 427 }
382 if (bin1 == bin2) { << 428 if(bin1 == bin2)
383 if (fAnalysisManager) { << 429 {
>> 430 if(fAnalysisManager)
>> 431 {
384 G4AnalysisManager::Instance()->FillH1(fH 432 G4AnalysisManager::Instance()->FillH1(fHistoId[3], z0, edep);
385 if (r1 < fStepR) { << 433 if(r1 < fStepR)
>> 434 {
386 G4double w = edep; 435 G4double w = edep;
387 if (r2 > fStepR) { << 436 if(r2 > fStepR)
>> 437 {
388 w *= (fStepR - r1) / (r2 - r1); 438 w *= (fStepR - r1) / (r2 - r1);
389 } 439 }
390 G4AnalysisManager::Instance()->FillH1( 440 G4AnalysisManager::Instance()->FillH1(fHistoId[4], z0, w);
391 } 441 }
392 } 442 }
393 } 443 }
394 else { << 444 else
>> 445 {
395 G4int bin; 446 G4int bin;
396 447
397 if (bin2 < bin1) { << 448 if(bin2 < bin1)
398 bin = bin2; << 449 {
>> 450 bin = bin2;
399 G4double z = z2; 451 G4double z = z2;
400 bin2 = bin1; << 452 bin2 = bin1;
401 z2 = z1; << 453 z2 = z1;
402 bin1 = bin; << 454 bin1 = bin;
403 z1 = z; << 455 z1 = z;
404 } 456 }
405 G4double zz1 = z1; 457 G4double zz1 = z1;
406 G4double zz2 = (bin1 + 1) * fStepZ; 458 G4double zz2 = (bin1 + 1) * fStepZ;
407 G4double rr1 = r1; 459 G4double rr1 = r1;
408 G4double dz = z2 - z1; << 460 G4double dz = z2 - z1;
409 G4double dr = r2 - r1; << 461 G4double dr = r2 - r1;
410 G4double rr2 = r1 + dr * (zz2 - zz1) / dz; 462 G4double rr2 = r1 + dr * (zz2 - zz1) / dz;
411 for (bin = bin1; bin <= bin2; bin++) { << 463 for(bin = bin1; bin <= bin2; bin++)
412 if (fAnalysisManager) { << 464 {
>> 465 if(fAnalysisManager)
>> 466 {
413 G4double de = edep * (zz2 - zz1) / dz; 467 G4double de = edep * (zz2 - zz1) / dz;
414 G4double zf = (zz1 + zz2) * 0.5; 468 G4double zf = (zz1 + zz2) * 0.5;
415 { 469 {
416 G4AnalysisManager::Instance()->FillH 470 G4AnalysisManager::Instance()->FillH1(fHistoId[3], zf, de);
417 } 471 }
418 if (rr1 < fStepR) { << 472 if(rr1 < fStepR)
>> 473 {
419 G4double w = de; 474 G4double w = de;
420 if (rr2 > fStepR) w *= (fStepR - rr1 << 475 if(rr2 > fStepR)
>> 476 w *= (fStepR - rr1) / (rr2 - rr1);
421 { 477 {
422 G4AnalysisManager::Instance()->Fil 478 G4AnalysisManager::Instance()->FillH1(fHistoId[4], zf, w);
423 } 479 }
424 } 480 }
425 } 481 }
426 zz1 = zz2; 482 zz1 = zz2;
427 zz2 = std::min(z2, zz1 + fStepZ); 483 zz2 = std::min(z2, zz1 + fStepZ);
428 rr1 = rr2; 484 rr1 = rr2;
429 rr2 = rr1 + dr * (zz2 - zz1) / dz; 485 rr2 = rr1 + dr * (zz2 - zz1) / dz;
430 } 486 }
431 } 487 }
432 } 488 }
433 489
434 //....oooOO0OOooo........oooOO0OOooo........oo 490 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
435 491