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