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