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