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