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) 182 << G4endl; << 191 << "Average number of g produced in the target " << xgt << G4endl; 183 G4cout << std::setprecision(4) << "Average n << 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 << "========================================================" 184 << G4endl; 202 << 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 } >> 253 >> 254 delete fAnalysisManager; >> 255 fAnalysisManager = 0; 228 } 256 } 229 } 257 } 230 258 231 //....oooOO0OOooo........oooOO0OOooo........oo 259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 232 260 233 void Run::AddTargetPhoton(const G4DynamicParti 261 void Run::AddTargetPhoton(const G4DynamicParticle* p) 234 { 262 { 235 ++fNgamTar; 263 ++fNgamTar; 236 if (fAnalysisManager) { << 264 if(fAnalysisManager) >> 265 { 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*) { ++fNgamPh; } 244 { << 245 ++fNgamPh; << 246 } << 247 273 248 //....oooOO0OOooo........oooOO0OOooo........oo 274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 249 275 250 void Run::AddTargetElectron(const G4DynamicPar 276 void Run::AddTargetElectron(const G4DynamicParticle* p) 251 { 277 { 252 ++fNeTar; 278 ++fNeTar; 253 if (fAnalysisManager) { << 279 if(fAnalysisManager) >> 280 { 254 fAnalysisManager->FillH1(fHistoId[8], p->G 281 fAnalysisManager->FillH1(fHistoId[8], p->GetKineticEnergy() / MeV, 1.0); 255 } 282 } 256 } 283 } 257 284 258 //....oooOO0OOooo........oooOO0OOooo........oo 285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 259 286 260 void Run::AddPhantomElectron(const G4DynamicPa 287 void Run::AddPhantomElectron(const G4DynamicParticle* p) 261 { 288 { 262 ++fNePh; 289 ++fNePh; 263 if (fAnalysisManager) { << 290 if(fAnalysisManager) >> 291 { 264 fAnalysisManager->FillH1(fHistoId[7], p->G 292 fAnalysisManager->FillH1(fHistoId[7], p->GetKineticEnergy() / MeV, 1.0); 265 } 293 } 266 } 294 } 267 295 268 //....oooOO0OOooo........oooOO0OOooo........oo 296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 269 297 270 void Run::ScoreNewTrack(const G4Track* aTrack) 298 void Run::ScoreNewTrack(const G4Track* aTrack) 271 { 299 { 272 // Save primary parameters 300 // Save primary parameters 273 const G4ParticleDefinition* particle = aTrac 301 const G4ParticleDefinition* particle = aTrack->GetParticleDefinition(); 274 G4int pid = aTrack->GetParentID(); << 302 G4int pid = aTrack->GetParentID(); 275 G4VPhysicalVolume* pv = aTrack->GetVolume(); << 303 G4VPhysicalVolume* pv = aTrack->GetVolume(); 276 const G4DynamicParticle* dp = aTrack->GetDyn << 304 const G4DynamicParticle* dp = aTrack->GetDynamicParticle(); 277 305 278 // primary particle 306 // primary particle 279 if (0 == pid) { << 307 if(0 == pid) >> 308 { 280 ++fNevt; 309 ++fNevt; 281 if (fVerbose) { << 310 if(fVerbose) >> 311 { 282 G4ThreeVector pos = aTrack->GetVertexPos 312 G4ThreeVector pos = aTrack->GetVertexPosition(); 283 G4ThreeVector dir = aTrack->GetMomentumD 313 G4ThreeVector dir = aTrack->GetMomentumDirection(); 284 G4cout << "TrackingAction: Primary " << 314 G4cout << "TrackingAction: Primary " << particle->GetParticleName() 285 << " Ekin(MeV)= " << aTrack->GetK << 315 << " Ekin(MeV)= " << aTrack->GetKineticEnergy() / MeV 286 << "; dir= " << dir << G4endl; << 316 << "; pos= " << pos << "; dir= " << dir << G4endl; 287 } 317 } 288 // secondary electron 318 // secondary electron 289 } 319 } 290 else if (0 < pid && particle == fElectron) { << 320 else if(0 < pid && particle == fElectron) 291 if (fVerbose) { << 321 { >> 322 if(fVerbose) >> 323 { 292 G4cout << "TrackingAction: Secondary ele 324 G4cout << "TrackingAction: Secondary electron " << G4endl; 293 } 325 } 294 AddElectron(); 326 AddElectron(); 295 if (pv == fPhantom) { << 327 if(pv == fPhantom) >> 328 { 296 AddPhantomElectron(dp); 329 AddPhantomElectron(dp); 297 } 330 } 298 else if (pv == fTarget1 || pv == fTarget2) << 331 else if(pv == fTarget1 || pv == fTarget2) >> 332 { 299 AddTargetElectron(dp); 333 AddTargetElectron(dp); 300 } 334 } 301 335 302 // secondary positron 336 // secondary positron 303 } 337 } 304 else if (0 < pid && particle == fPositron) { << 338 else if(0 < pid && particle == fPositron) 305 if (fVerbose) { << 339 { >> 340 if(fVerbose) >> 341 { 306 G4cout << "TrackingAction: Secondary pos 342 G4cout << "TrackingAction: Secondary positron " << G4endl; 307 } 343 } 308 AddPositron(); 344 AddPositron(); 309 345 310 // secondary gamma 346 // secondary gamma 311 } 347 } 312 else if (0 < pid && particle == fGamma) { << 348 else if(0 < pid && particle == fGamma) 313 if (fVerbose) { << 349 { >> 350 if(fVerbose) >> 351 { 314 G4cout << "TrackingAction: Secondary gam 352 G4cout << "TrackingAction: Secondary gamma; parentID= " << pid 315 << " E= " << aTrack->GetKineticEn 353 << " E= " << aTrack->GetKineticEnergy() << G4endl; 316 } 354 } 317 AddPhoton(); 355 AddPhoton(); 318 if (pv == fPhantom) { << 356 if(pv == fPhantom) >> 357 { 319 AddPhantomPhoton(dp); 358 AddPhantomPhoton(dp); 320 } 359 } 321 else if (pv == fTarget1 || pv == fTarget2) << 360 else if(pv == fTarget1 || pv == fTarget2) >> 361 { 322 AddTargetPhoton(dp); 362 AddTargetPhoton(dp); 323 } 363 } 324 } 364 } 325 } 365 } 326 366 327 //....oooOO0OOooo........oooOO0OOooo........oo 367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 328 368 329 void Run::AddPhantomGamma(G4double e, G4double 369 void Run::AddPhantomGamma(G4double e, G4double r) 330 { 370 { 331 e /= MeV; 371 e /= MeV; 332 fSumR += e; 372 fSumR += e; 333 G4int bin = (G4int)(e / fStepE); 373 G4int bin = (G4int)(e / fStepE); 334 if (bin >= fNBinsE) { << 374 if(bin >= fNBinsE) >> 375 { 335 bin = fNBinsE - 1; 376 bin = fNBinsE - 1; 336 } 377 } 337 fGammaE[bin] += e; 378 fGammaE[bin] += e; 338 G4int bin1 = (G4int)(r / fStepR); 379 G4int bin1 = (G4int)(r / fStepR); 339 if (bin1 >= fNBinsR) { << 380 if(bin1 >= fNBinsR) >> 381 { 340 bin1 = fNBinsR - 1; 382 bin1 = fNBinsR - 1; 341 } 383 } 342 if (fAnalysisManager) { << 384 if(fAnalysisManager) >> 385 { 343 G4AnalysisManager::Instance()->FillH1(fHis 386 G4AnalysisManager::Instance()->FillH1(fHistoId[6], e, 1.0); 344 G4AnalysisManager::Instance()->FillH1(fHis 387 G4AnalysisManager::Instance()->FillH1(fHistoId[9], r, e * fVolumeR[bin1]); 345 } 388 } 346 } 389 } 347 390 348 //....oooOO0OOooo........oooOO0OOooo........oo 391 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 349 392 350 void Run::AddPhantomStep(G4double edep, G4doub << 393 void Run::AddPhantomStep(G4double edep, G4double r1, G4double z1, G4double r2, 351 G4double r0, G4double << 394 G4double z2, G4double r0, G4double z0) 352 { 395 { 353 ++fNstep; 396 ++fNstep; 354 G4int nzbin = (G4int)(z0 / fStepZ); 397 G4int nzbin = (G4int)(z0 / fStepZ); 355 if (fVerbose) { << 398 if(fVerbose) 356 G4cout << "Histo: edep(MeV)= " << edep / M << 399 { 357 << " z1= " << z1 << " r2= " << r2 < << 400 G4cout << "Histo: edep(MeV)= " << edep / MeV << " at binz= " << nzbin 358 << G4endl; << 401 << " r1= " << r1 << " z1= " << z1 << " r2= " << r2 << " z2= " << z2 >> 402 << " r0= " << r0 << " z0= " << z0 << G4endl; 359 } 403 } 360 if (nzbin == fScoreBin) { << 404 if(nzbin == fScoreBin) >> 405 { 361 G4int bin = (G4int)(r0 / fStepR); 406 G4int bin = (G4int)(r0 / fStepR); 362 if (bin >= fNBinsR) { << 407 if(bin >= fNBinsR) >> 408 { 363 bin = fNBinsR - 1; 409 bin = fNBinsR - 1; 364 } 410 } 365 double w = edep * fVolumeR[bin]; 411 double w = edep * fVolumeR[bin]; 366 fEdep[bin] += w; 412 fEdep[bin] += w; 367 413 368 if (fAnalysisManager) { << 414 if(fAnalysisManager) >> 415 { 369 G4AnalysisManager::Instance()->FillH1(fH 416 G4AnalysisManager::Instance()->FillH1(fHistoId[0], r0, w); 370 G4AnalysisManager::Instance()->FillH1(fH 417 G4AnalysisManager::Instance()->FillH1(fHistoId[1], r0, w); 371 G4AnalysisManager::Instance()->FillH1(fH 418 G4AnalysisManager::Instance()->FillH1(fHistoId[2], r0, w); 372 } 419 } 373 } 420 } 374 G4int bin1 = (G4int)(z1 / fStepZ); 421 G4int bin1 = (G4int)(z1 / fStepZ); 375 if (bin1 >= fNBinsZ) { << 422 if(bin1 >= fNBinsZ) >> 423 { 376 bin1 = fNBinsZ - 1; 424 bin1 = fNBinsZ - 1; 377 } 425 } 378 G4int bin2 = (G4int)(z2 / fStepZ); 426 G4int bin2 = (G4int)(z2 / fStepZ); 379 if (bin2 >= fNBinsZ) { << 427 if(bin2 >= fNBinsZ) >> 428 { 380 bin2 = fNBinsZ - 1; 429 bin2 = fNBinsZ - 1; 381 } 430 } 382 if (bin1 == bin2) { << 431 if(bin1 == bin2) 383 if (fAnalysisManager) { << 432 { >> 433 if(fAnalysisManager) >> 434 { 384 G4AnalysisManager::Instance()->FillH1(fH 435 G4AnalysisManager::Instance()->FillH1(fHistoId[3], z0, edep); 385 if (r1 < fStepR) { << 436 if(r1 < fStepR) >> 437 { 386 G4double w = edep; 438 G4double w = edep; 387 if (r2 > fStepR) { << 439 if(r2 > fStepR) >> 440 { 388 w *= (fStepR - r1) / (r2 - r1); 441 w *= (fStepR - r1) / (r2 - r1); 389 } 442 } 390 G4AnalysisManager::Instance()->FillH1( 443 G4AnalysisManager::Instance()->FillH1(fHistoId[4], z0, w); 391 } 444 } 392 } 445 } 393 } 446 } 394 else { << 447 else >> 448 { 395 G4int bin; 449 G4int bin; 396 450 397 if (bin2 < bin1) { << 451 if(bin2 < bin1) 398 bin = bin2; << 452 { >> 453 bin = bin2; 399 G4double z = z2; 454 G4double z = z2; 400 bin2 = bin1; << 455 bin2 = bin1; 401 z2 = z1; << 456 z2 = z1; 402 bin1 = bin; << 457 bin1 = bin; 403 z1 = z; << 458 z1 = z; 404 } 459 } 405 G4double zz1 = z1; 460 G4double zz1 = z1; 406 G4double zz2 = (bin1 + 1) * fStepZ; 461 G4double zz2 = (bin1 + 1) * fStepZ; 407 G4double rr1 = r1; 462 G4double rr1 = r1; 408 G4double dz = z2 - z1; << 463 G4double dz = z2 - z1; 409 G4double dr = r2 - r1; << 464 G4double dr = r2 - r1; 410 G4double rr2 = r1 + dr * (zz2 - zz1) / dz; 465 G4double rr2 = r1 + dr * (zz2 - zz1) / dz; 411 for (bin = bin1; bin <= bin2; bin++) { << 466 for(bin = bin1; bin <= bin2; bin++) 412 if (fAnalysisManager) { << 467 { >> 468 if(fAnalysisManager) >> 469 { 413 G4double de = edep * (zz2 - zz1) / dz; 470 G4double de = edep * (zz2 - zz1) / dz; 414 G4double zf = (zz1 + zz2) * 0.5; 471 G4double zf = (zz1 + zz2) * 0.5; 415 { 472 { 416 G4AnalysisManager::Instance()->FillH 473 G4AnalysisManager::Instance()->FillH1(fHistoId[3], zf, de); 417 } 474 } 418 if (rr1 < fStepR) { << 475 if(rr1 < fStepR) >> 476 { 419 G4double w = de; 477 G4double w = de; 420 if (rr2 > fStepR) w *= (fStepR - rr1 << 478 if(rr2 > fStepR) >> 479 w *= (fStepR - rr1) / (rr2 - rr1); 421 { 480 { 422 G4AnalysisManager::Instance()->Fil 481 G4AnalysisManager::Instance()->FillH1(fHistoId[4], zf, w); 423 } 482 } 424 } 483 } 425 } 484 } 426 zz1 = zz2; 485 zz1 = zz2; 427 zz2 = std::min(z2, zz1 + fStepZ); 486 zz2 = std::min(z2, zz1 + fStepZ); 428 rr1 = rr2; 487 rr1 = rr2; 429 rr2 = rr1 + dr * (zz2 - zz1) / dz; 488 rr2 = rr1 + dr * (zz2 - zz1) / dz; 430 } 489 } 431 } 490 } 432 } 491 } 433 492 434 //....oooOO0OOooo........oooOO0OOooo........oo 493 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 435 494