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 // 26 // 27 // ------------------------------------------- 27 // ------------------------------------------------------------------- 28 // 28 // 29 // GEANT4 Class class file 29 // GEANT4 Class class file 30 // 30 // 31 // 31 // 32 // File name: G4VAtomDeexcitation 32 // File name: G4VAtomDeexcitation 33 // 33 // 34 // Author: Alfonso Mantero & Vladimir I 34 // Author: Alfonso Mantero & Vladimir Ivanchenko 35 // 35 // 36 // Creation date: 21.04.2010 36 // Creation date: 21.04.2010 37 // 37 // 38 // Modifications: 38 // Modifications: 39 // 39 // 40 // Class Description: 40 // Class Description: 41 // 41 // 42 // Abstract interface to energy loss models 42 // Abstract interface to energy loss models 43 43 44 // ------------------------------------------- 44 // ------------------------------------------------------------------- 45 // 45 // 46 46 47 #include "G4VAtomDeexcitation.hh" 47 #include "G4VAtomDeexcitation.hh" 48 #include "G4SystemOfUnits.hh" 48 #include "G4SystemOfUnits.hh" 49 #include "G4EmParameters.hh" 49 #include "G4EmParameters.hh" 50 #include "G4ParticleDefinition.hh" 50 #include "G4ParticleDefinition.hh" 51 #include "G4DynamicParticle.hh" 51 #include "G4DynamicParticle.hh" 52 #include "G4Step.hh" 52 #include "G4Step.hh" 53 #include "G4Region.hh" 53 #include "G4Region.hh" 54 #include "G4RegionStore.hh" 54 #include "G4RegionStore.hh" 55 #include "G4MaterialCutsCouple.hh" 55 #include "G4MaterialCutsCouple.hh" 56 #include "G4MaterialCutsCouple.hh" 56 #include "G4MaterialCutsCouple.hh" 57 #include "G4Material.hh" 57 #include "G4Material.hh" 58 #include "G4Element.hh" 58 #include "G4Element.hh" 59 #include "G4ElementVector.hh" 59 #include "G4ElementVector.hh" 60 #include "Randomize.hh" 60 #include "Randomize.hh" 61 #include "G4VParticleChange.hh" 61 #include "G4VParticleChange.hh" 62 #include "G4EmSecondaryParticleType.hh" 62 #include "G4EmSecondaryParticleType.hh" 63 #include "G4Gamma.hh" 63 #include "G4Gamma.hh" 64 #include "G4Log.hh" 64 #include "G4Log.hh" 65 65 66 #ifdef G4MULTITHREADED 66 #ifdef G4MULTITHREADED 67 G4Mutex G4VAtomDeexcitation::atomDeexcitatio 67 G4Mutex G4VAtomDeexcitation::atomDeexcitationMutex = G4MUTEX_INITIALIZER; 68 #endif 68 #endif 69 69 70 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 71 71 72 G4VAtomDeexcitation::G4VAtomDeexcitation(const 72 G4VAtomDeexcitation::G4VAtomDeexcitation(const G4String& modname) 73 : name(modname) 73 : name(modname) 74 { 74 { 75 vdyn.reserve(5); 75 vdyn.reserve(5); 76 theCoupleTable = nullptr; 76 theCoupleTable = nullptr; 77 gamma = G4Gamma::Gamma(); 77 gamma = G4Gamma::Gamma(); 78 } 78 } 79 79 80 //....oooOO0OOooo........oooOO0OOooo........oo 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 81 81 82 G4VAtomDeexcitation::~G4VAtomDeexcitation() = 82 G4VAtomDeexcitation::~G4VAtomDeexcitation() = default; 83 83 84 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 85 85 86 void G4VAtomDeexcitation::InitialiseAtomicDeex 86 void G4VAtomDeexcitation::InitialiseAtomicDeexcitation() 87 { 87 { 88 G4EmParameters* theParameters = G4EmParamete 88 G4EmParameters* theParameters = G4EmParameters::Instance(); 89 theParameters->DefineRegParamForDeex(this); 89 theParameters->DefineRegParamForDeex(this); 90 90 91 // Define list of couples 91 // Define list of couples 92 theCoupleTable = G4ProductionCutsTable::GetP 92 theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable(); 93 nCouples = (G4int)theCoupleTable->GetTableSi 93 nCouples = (G4int)theCoupleTable->GetTableSize(); 94 94 95 // needed for unit tests 95 // needed for unit tests 96 std::size_t nn = std::max(nCouples, 1); 96 std::size_t nn = std::max(nCouples, 1); 97 if(activeDeexcitationMedia.size() != nn) { 97 if(activeDeexcitationMedia.size() != nn) { 98 activeDeexcitationMedia.resize(nn, false); 98 activeDeexcitationMedia.resize(nn, false); 99 activeAugerMedia.resize(nn, false); 99 activeAugerMedia.resize(nn, false); 100 activePIXEMedia.resize(nn, false); 100 activePIXEMedia.resize(nn, false); 101 } 101 } 102 if(activeZ.size() != 93) { activeZ.resize(93 102 if(activeZ.size() != 93) { activeZ.resize(93, false); } 103 103 104 // initialisation of flags and options 104 // initialisation of flags and options 105 // normally there is no locksed flags 105 // normally there is no locksed flags 106 if(!isActiveLocked) { isActive = theParamet 106 if(!isActiveLocked) { isActive = theParameters->Fluo(); } 107 if(!isAugerLocked) { flagAuger = theParamet 107 if(!isAugerLocked) { flagAuger = theParameters->Auger(); } 108 if(!isPIXELocked) { flagPIXE = theParamet 108 if(!isPIXELocked) { flagPIXE = theParameters->Pixe(); } 109 ignoreCuts = theParameters->DeexcitationIgno 109 ignoreCuts = theParameters->DeexcitationIgnoreCut(); 110 110 111 // Define list of regions 111 // Define list of regions 112 std::size_t nRegions = deRegions.size(); 112 std::size_t nRegions = deRegions.size(); 113 // check if deexcitation is active for the g 113 // check if deexcitation is active for the given run 114 if(!isActive && 0 == nRegions) { return; } 114 if(!isActive && 0 == nRegions) { return; } 115 115 116 // if no active regions add a world 116 // if no active regions add a world 117 if(0 == nRegions) { 117 if(0 == nRegions) { 118 SetDeexcitationActiveRegion("World",isActi 118 SetDeexcitationActiveRegion("World",isActive,flagAuger,flagPIXE); 119 nRegions = deRegions.size(); 119 nRegions = deRegions.size(); 120 } 120 } 121 121 122 if(0 < verbose) { 122 if(0 < verbose) { 123 G4cout << G4endl; 123 G4cout << G4endl; 124 G4cout << "### === Deexcitation model " < 124 G4cout << "### === Deexcitation model " << name 125 << " is activated for " << nRegions 125 << " is activated for " << nRegions; 126 if(1 == nRegions) { G4cout << " region:" < 126 if(1 == nRegions) { G4cout << " region:" << G4endl; } 127 else { G4cout << " regions:" 127 else { G4cout << " regions:" << G4endl;} 128 } 128 } 129 129 130 // Identify active media 130 // Identify active media 131 const G4RegionStore* regionStore = G4RegionS 131 const G4RegionStore* regionStore = G4RegionStore::GetInstance(); 132 for(std::size_t j=0; j<nRegions; ++j) { 132 for(std::size_t j=0; j<nRegions; ++j) { 133 const G4Region* reg = regionStore->GetRegi 133 const G4Region* reg = regionStore->GetRegion(activeRegions[j], false); 134 if(nullptr != reg && 0 < nCouples) { 134 if(nullptr != reg && 0 < nCouples) { 135 const G4ProductionCuts* rpcuts = reg->Ge 135 const G4ProductionCuts* rpcuts = reg->GetProductionCuts(); 136 if(0 < verbose) { 136 if(0 < verbose) { 137 G4cout << " " << activeRegion 137 G4cout << " " << activeRegions[j] 138 << " " << deRegions[j] << " 138 << " " << deRegions[j] << " " << AugerRegions[j] 139 << " " << PIXERegions[j] << G4 139 << " " << PIXERegions[j] << G4endl; 140 } 140 } 141 for(G4int i=0; i<nCouples; ++i) { 141 for(G4int i=0; i<nCouples; ++i) { 142 const G4MaterialCutsCouple* couple = 142 const G4MaterialCutsCouple* couple = 143 theCoupleTable->GetMaterialCutsCoupl 143 theCoupleTable->GetMaterialCutsCouple(i); 144 if (couple->GetProductionCuts() == rpc 144 if (couple->GetProductionCuts() == rpcuts) { 145 activeDeexcitationMedia[i] = deRegio 145 activeDeexcitationMedia[i] = deRegions[j]; 146 activeAugerMedia[i] = AugerRegions[j 146 activeAugerMedia[i] = AugerRegions[j]; 147 activePIXEMedia[i] = PIXERegions[j]; 147 activePIXEMedia[i] = PIXERegions[j]; 148 } 148 } 149 } 149 } 150 } 150 } 151 } 151 } 152 std::size_t nelm = G4Element::GetNumberOfEle 152 std::size_t nelm = G4Element::GetNumberOfElements(); 153 //G4cout << nelm << G4endl; 153 //G4cout << nelm << G4endl; 154 for(std::size_t k=0; k<nelm; ++k) { 154 for(std::size_t k=0; k<nelm; ++k) { 155 G4int Z = (*(G4Element::GetElementTable()) 155 G4int Z = (*(G4Element::GetElementTable()))[k]->GetZasInt(); 156 if(Z > 5 && Z < 93) { 156 if(Z > 5 && Z < 93) { 157 activeZ[Z] = true; 157 activeZ[Z] = true; 158 //G4cout << "!!! Active de-excitation Z= 158 //G4cout << "!!! Active de-excitation Z= " << Z << G4endl; 159 } 159 } 160 } 160 } 161 161 162 // Initialise derived class 162 // Initialise derived class 163 InitialiseForNewRun(); 163 InitialiseForNewRun(); 164 164 165 if(0 < verbose && flagAuger) { 165 if(0 < verbose && flagAuger) { 166 G4cout << "### === Auger flag: " << flagA 166 G4cout << "### === Auger flag: " << flagAuger 167 << G4endl; 167 << G4endl; 168 } 168 } 169 if(0 < verbose) { 169 if(0 < verbose) { 170 G4cout << "### === Ignore cuts flag: " 170 G4cout << "### === Ignore cuts flag: " << ignoreCuts 171 << G4endl; 171 << G4endl; 172 } 172 } 173 if(0 < verbose && flagPIXE) { 173 if(0 < verbose && flagPIXE) { 174 G4cout << "### === PIXE model for hadrons 174 G4cout << "### === PIXE model for hadrons: " 175 << theParameters->PIXECrossSectionM 175 << theParameters->PIXECrossSectionModel() 176 << G4endl; 176 << G4endl; 177 G4cout << "### === PIXE model for e+-: 177 G4cout << "### === PIXE model for e+-: " 178 << theParameters->PIXEElectronCross 178 << theParameters->PIXEElectronCrossSectionModel() 179 << G4endl; 179 << G4endl; 180 } 180 } 181 } 181 } 182 182 183 //....oooOO0OOooo........oooOO0OOooo........oo 183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 184 184 185 void 185 void 186 G4VAtomDeexcitation::SetDeexcitationActiveRegi 186 G4VAtomDeexcitation::SetDeexcitationActiveRegion(const G4String& rname, 187 187 G4bool valDeexcitation, 188 188 G4bool valAuger, 189 189 G4bool valPIXE) 190 { 190 { 191 // no PIXE in parallel world 191 // no PIXE in parallel world 192 if(rname == "DefaultRegionForParallelWorld") 192 if(rname == "DefaultRegionForParallelWorld") { return; } 193 193 194 G4String ss = rname; 194 G4String ss = rname; 195 /* 195 /* 196 G4cout << "### G4VAtomDeexcitation::SetDeexc 196 G4cout << "### G4VAtomDeexcitation::SetDeexcitationActiveRegion " << ss 197 << " " << valDeexcitation << " " << 197 << " " << valDeexcitation << " " << valAuger 198 << " " << valPIXE << G4endl; 198 << " " << valPIXE << G4endl; 199 */ 199 */ 200 if(ss == "world" || ss == "World" || ss == " 200 if(ss == "world" || ss == "World" || ss == "WORLD") { 201 ss = "DefaultRegionForTheWorld"; 201 ss = "DefaultRegionForTheWorld"; 202 } 202 } 203 std::size_t n = deRegions.size(); 203 std::size_t n = deRegions.size(); 204 for(std::size_t i=0; i<n; ++i) { 204 for(std::size_t i=0; i<n; ++i) { 205 205 206 // Region already exist 206 // Region already exist 207 if(ss == activeRegions[i]) { 207 if(ss == activeRegions[i]) { 208 deRegions[i] = valDeexcitation; 208 deRegions[i] = valDeexcitation; 209 AugerRegions[i] = valAuger; 209 AugerRegions[i] = valAuger; 210 PIXERegions[i] = valPIXE; 210 PIXERegions[i] = valPIXE; 211 return; 211 return; 212 } 212 } 213 } 213 } 214 // New region 214 // New region 215 activeRegions.push_back(ss); 215 activeRegions.push_back(ss); 216 deRegions.push_back(valDeexcitation); 216 deRegions.push_back(valDeexcitation); 217 AugerRegions.push_back(valAuger); 217 AugerRegions.push_back(valAuger); 218 PIXERegions.push_back(valPIXE); 218 PIXERegions.push_back(valPIXE); 219 219 220 // if de-excitation defined for the world vo 220 // if de-excitation defined for the world volume 221 // it should be active for all G4Regions 221 // it should be active for all G4Regions 222 if(ss == "DefaultRegionForTheWorld") { 222 if(ss == "DefaultRegionForTheWorld") { 223 G4RegionStore* regions = G4RegionStore::Ge 223 G4RegionStore* regions = G4RegionStore::GetInstance(); 224 std::size_t nn = regions->size(); 224 std::size_t nn = regions->size(); 225 for(std::size_t i=0; i<nn; ++i) { 225 for(std::size_t i=0; i<nn; ++i) { 226 if(ss == (*regions)[i]->GetName()) { con 226 if(ss == (*regions)[i]->GetName()) { continue; } 227 SetDeexcitationActiveRegion((*regions)[i 227 SetDeexcitationActiveRegion((*regions)[i]->GetName(), valDeexcitation, 228 valAuger, va 228 valAuger, valPIXE); 229 229 230 } 230 } 231 } 231 } 232 } 232 } 233 233 234 void G4VAtomDeexcitation::GenerateParticles(st 234 void G4VAtomDeexcitation::GenerateParticles(std::vector<G4DynamicParticle*>* v, 235 co 235 const G4AtomicShell* as, 236 G4 236 G4int Z, G4int idx) 237 { 237 { 238 G4double gCut = DBL_MAX; 238 G4double gCut = DBL_MAX; 239 if(ignoreCuts) { 239 if(ignoreCuts) { 240 gCut = 0.0; 240 gCut = 0.0; 241 } else if (nullptr != theCoupleTable) { 241 } else if (nullptr != theCoupleTable) { 242 gCut = (*(theCoupleTable->GetEnergyCutsVec 242 gCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[idx]; 243 } 243 } 244 if(gCut < as->BindingEnergy()) { 244 if(gCut < as->BindingEnergy()) { 245 G4double eCut = DBL_MAX; 245 G4double eCut = DBL_MAX; 246 if(CheckAugerActiveRegion(idx)) { 246 if(CheckAugerActiveRegion(idx)) { 247 if(ignoreCuts) { 247 if(ignoreCuts) { 248 eCut = 0.0; 248 eCut = 0.0; 249 } else if (nullptr != theCoupleTable) { 249 } else if (nullptr != theCoupleTable) { 250 eCut = (*(theCoupleTable->GetEnergyCut 250 eCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[idx]; 251 } 251 } 252 } 252 } 253 GenerateParticles(v, as, Z, gCut, eCut); 253 GenerateParticles(v, as, Z, gCut, eCut); 254 } 254 } 255 } 255 } 256 256 257 //....oooOO0OOooo........oooOO0OOooo........oo 257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 258 258 259 void 259 void 260 G4VAtomDeexcitation::AlongStepDeexcitation(std 260 G4VAtomDeexcitation::AlongStepDeexcitation(std::vector<G4Track*>& tracks, 261 con 261 const G4Step& step, 262 G4d 262 G4double& eLossMax, 263 G4i 263 G4int coupleIndex) 264 { 264 { 265 G4double truelength = step.GetStepLength(); 265 G4double truelength = step.GetStepLength(); 266 if(!flagPIXE && !activePIXEMedia[coupleIndex 266 if(!flagPIXE && !activePIXEMedia[coupleIndex]) { return; } 267 if(eLossMax <= 0.0 || truelength <= 0.0) 267 if(eLossMax <= 0.0 || truelength <= 0.0) { return; } 268 268 269 // step parameters 269 // step parameters 270 const G4StepPoint* preStep = step.GetPreStep 270 const G4StepPoint* preStep = step.GetPreStepPoint(); 271 const G4ThreeVector prePos = preStep->GetPos 271 const G4ThreeVector prePos = preStep->GetPosition(); 272 const G4ThreeVector delta = step.GetPostStep 272 const G4ThreeVector delta = step.GetPostStepPoint()->GetPosition() - prePos; 273 const G4double preTime = preStep->GetGlobalT 273 const G4double preTime = preStep->GetGlobalTime(); 274 const G4double dt = step.GetPostStepPoint()- 274 const G4double dt = step.GetPostStepPoint()->GetGlobalTime() - preTime; 275 275 276 // particle parameters 276 // particle parameters 277 const G4Track* track = step.GetTrack(); 277 const G4Track* track = step.GetTrack(); 278 const G4ParticleDefinition* part = track->Ge 278 const G4ParticleDefinition* part = track->GetDefinition(); 279 G4double ekin = preStep->GetKineticEnergy(); 279 G4double ekin = preStep->GetKineticEnergy(); 280 280 281 // media parameters 281 // media parameters 282 G4double gCut = (*theCoupleTable->GetEnergyC 282 G4double gCut = (*theCoupleTable->GetEnergyCutsVector(0))[coupleIndex]; 283 if(ignoreCuts) { gCut = 0.0; } 283 if(ignoreCuts) { gCut = 0.0; } 284 G4double eCut = DBL_MAX; 284 G4double eCut = DBL_MAX; 285 if(CheckAugerActiveRegion(coupleIndex)) { 285 if(CheckAugerActiveRegion(coupleIndex)) { 286 eCut = (*theCoupleTable->GetEnergyCutsVect 286 eCut = (*theCoupleTable->GetEnergyCutsVector(1))[coupleIndex]; 287 if(ignoreCuts) { eCut = 0.0; } 287 if(ignoreCuts) { eCut = 0.0; } 288 } 288 } 289 289 290 //G4cout<<"!Sample PIXE gCut(MeV)= "<<gCut<< 290 //G4cout<<"!Sample PIXE gCut(MeV)= "<<gCut<<" eCut(MeV)= "<<eCut 291 // <<" Ekin(MeV)= " << ekin/MeV << G4 291 // <<" Ekin(MeV)= " << ekin/MeV << G4endl; 292 292 293 const G4Material* material = preStep->GetMat 293 const G4Material* material = preStep->GetMaterial(); 294 const G4ElementVector* theElementVector = ma 294 const G4ElementVector* theElementVector = material->GetElementVector(); 295 const G4double* theAtomNumDensityVector = 295 const G4double* theAtomNumDensityVector = 296 material->GetVecNbOfAtomsPerVolume(); 296 material->GetVecNbOfAtomsPerVolume(); 297 const std::size_t nelm = material->GetNumber 297 const std::size_t nelm = material->GetNumberOfElements(); 298 298 299 // loop over deexcitations 299 // loop over deexcitations 300 for(std::size_t i=0; i<nelm; ++i) { 300 for(std::size_t i=0; i<nelm; ++i) { 301 G4int Z = (*theElementVector)[i]->GetZasIn 301 G4int Z = (*theElementVector)[i]->GetZasInt(); 302 if(activeZ[Z] && Z < 93) { 302 if(activeZ[Z] && Z < 93) { 303 G4int nshells = 303 G4int nshells = 304 std::min(9,(*theElementVector)[i]->Get 304 std::min(9,(*theElementVector)[i]->GetNbOfAtomicShells()); 305 G4double rho = truelength*theAtomNumDens 305 G4double rho = truelength*theAtomNumDensityVector[i]; 306 //G4cout<<" Z "<< Z <<" is active x(m 306 //G4cout<<" Z "<< Z <<" is active x(mm)= " << truelength/mm << G4endl; 307 for(G4int ii=0; ii<nshells; ++ii) { 307 for(G4int ii=0; ii<nshells; ++ii) { 308 auto as = (G4AtomicShellEnumerator)(ii 308 auto as = (G4AtomicShellEnumerator)(ii); 309 const G4AtomicShell* shell = GetAtomic 309 const G4AtomicShell* shell = GetAtomicShell(Z, as); 310 const G4double bindingEnergy = shell-> 310 const G4double bindingEnergy = shell->BindingEnergy(); 311 311 312 if(gCut > bindingEnergy) { break; } 312 if(gCut > bindingEnergy) { break; } 313 313 314 if(eLossMax > bindingEnergy) { 314 if(eLossMax > bindingEnergy) { 315 G4double sig = rho* 315 G4double sig = rho* 316 GetShellIonisationCrossSectionPerA 316 GetShellIonisationCrossSectionPerAtom(part, Z, as, ekin, material); 317 317 318 // mfp is mean free path in units of 318 // mfp is mean free path in units of step size 319 if(sig > 0.0) { 319 if(sig > 0.0) { 320 G4double mfp = 1.0/sig; 320 G4double mfp = 1.0/sig; 321 G4double stot = 0.0; 321 G4double stot = 0.0; 322 //G4cout << " Shell " << ii << " m 322 //G4cout << " Shell " << ii << " mfp(mm)= " << mfp/mm << G4endl; 323 // sample ionisation points 323 // sample ionisation points 324 do { 324 do { 325 stot -= mfp*G4Log(G4UniformRand( 325 stot -= mfp*G4Log(G4UniformRand()); 326 if( stot > 1.0 || eLossMax < bin 326 if( stot > 1.0 || eLossMax < bindingEnergy) { break; } 327 // sample deexcitation 327 // sample deexcitation 328 vdyn.clear(); 328 vdyn.clear(); 329 GenerateParticles(&vdyn, shell, 329 GenerateParticles(&vdyn, shell, Z, gCut, eCut); 330 std::size_t nsec = vdyn.size(); 330 std::size_t nsec = vdyn.size(); 331 if(nsec > 0) { 331 if(nsec > 0) { 332 G4ThreeVector r = prePos + st 332 G4ThreeVector r = prePos + stot*delta; 333 G4double time = preTime + st 333 G4double time = preTime + stot*dt; 334 for(std::size_t j=0; j<nsec; + 334 for(std::size_t j=0; j<nsec; ++j) { 335 G4DynamicParticle* dp = vdyn 335 G4DynamicParticle* dp = vdyn[j]; 336 G4double e = dp->GetKineticE 336 G4double e = dp->GetKineticEnergy(); 337 337 338 // save new secondary if the 338 // save new secondary if there is enough energy 339 if(eLossMax >= e) { 339 if(eLossMax >= e) { 340 eLossMax -= e; 340 eLossMax -= e; 341 G4Track* t = new G4Track(d 341 G4Track* t = new G4Track(dp, time, r); 342 342 343 // defined secondary type 343 // defined secondary type 344 if(dp->GetDefinition() == 344 if(dp->GetDefinition() == gamma) { 345 t->SetCreatorModelID(_Ga 345 t->SetCreatorModelID(_GammaPIXE); 346 } else { 346 } else { 347 t->SetCreatorModelID(_eP 347 t->SetCreatorModelID(_ePIXE); 348 } 348 } 349 tracks.push_back(t); 349 tracks.push_back(t); 350 } else { 350 } else { 351 delete dp; 351 delete dp; 352 } 352 } 353 } 353 } 354 } 354 } 355 // Loop checking, 03-Aug-2015, V 355 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 356 } while (stot < 1.0); 356 } while (stot < 1.0); 357 } 357 } 358 } 358 } 359 } 359 } 360 } 360 } 361 } 361 } 362 return; 362 return; 363 } 363 } 364 364 365 //....oooOO0OOooo........oooOO0OOooo........oo 365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 366 366