| 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 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class class file
30 //
31 //
32 // File name: G4VAtomDeexcitation
33 //
34 // Author: Alfonso Mantero & Vladimir I
35 //
36 // Creation date: 21.04.2010
37 //
38 // Modifications:
39 //
40 // Class Description:
41 //
42 // Abstract interface to energy loss models
43
44 // -------------------------------------------
45 //
46
47 #include "G4VAtomDeexcitation.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4EmParameters.hh"
50 #include "G4ParticleDefinition.hh"
51 #include "G4DynamicParticle.hh"
52 #include "G4Step.hh"
53 #include "G4Region.hh"
54 #include "G4RegionStore.hh"
55 #include "G4MaterialCutsCouple.hh"
56 #include "G4MaterialCutsCouple.hh"
57 #include "G4Material.hh"
58 #include "G4Element.hh"
59 #include "G4ElementVector.hh"
60 #include "Randomize.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
70 //....oooOO0OOooo........oooOO0OOooo........oo
71
72 G4VAtomDeexcitation::G4VAtomDeexcitation(const
73 : name(modname)
74 {
75 vdyn.reserve(5);
76 theCoupleTable = nullptr;
77 gamma = G4Gamma::Gamma();
78 }
79
80 //....oooOO0OOooo........oooOO0OOooo........oo
81
82 G4VAtomDeexcitation::~G4VAtomDeexcitation() =
83
84 //....oooOO0OOooo........oooOO0OOooo........oo
85
86 void G4VAtomDeexcitation::InitialiseAtomicDeex
87 {
88 G4EmParameters* theParameters = G4EmParamete
89 theParameters->DefineRegParamForDeex(this);
90
91 // Define list of couples
92 theCoupleTable = G4ProductionCutsTable::GetP
93 nCouples = (G4int)theCoupleTable->GetTableSi
94
95 // needed for unit tests
96 std::size_t nn = std::max(nCouples, 1);
97 if(activeDeexcitationMedia.size() != nn) {
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
104 // initialisation of flags and options
105 // normally there is no locksed flags
106 if(!isActiveLocked) { isActive = theParamet
107 if(!isAugerLocked) { flagAuger = theParamet
108 if(!isPIXELocked) { flagPIXE = theParamet
109 ignoreCuts = theParameters->DeexcitationIgno
110
111 // Define list of regions
112 std::size_t nRegions = deRegions.size();
113 // check if deexcitation is active for the g
114 if(!isActive && 0 == nRegions) { return; }
115
116 // if no active regions add a world
117 if(0 == nRegions) {
118 SetDeexcitationActiveRegion("World",isActi
119 nRegions = deRegions.size();
120 }
121
122 if(0 < verbose) {
123 G4cout << G4endl;
124 G4cout << "### === Deexcitation model " <
125 << " is activated for " << nRegions
126 if(1 == nRegions) { G4cout << " region:" <
127 else { G4cout << " regions:"
128 }
129
130 // Identify active media
131 const G4RegionStore* regionStore = G4RegionS
132 for(std::size_t j=0; j<nRegions; ++j) {
133 const G4Region* reg = regionStore->GetRegi
134 if(nullptr != reg && 0 < nCouples) {
135 const G4ProductionCuts* rpcuts = reg->Ge
136 if(0 < verbose) {
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 }
151 }
152 std::size_t nelm = G4Element::GetNumberOfEle
153 //G4cout << nelm << G4endl;
154 for(std::size_t k=0; k<nelm; ++k) {
155 G4int Z = (*(G4Element::GetElementTable())
156 if(Z > 5 && Z < 93) {
157 activeZ[Z] = true;
158 //G4cout << "!!! Active de-excitation Z=
159 }
160 }
161
162 // Initialise derived class
163 InitialiseForNewRun();
164
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) {
174 G4cout << "### === PIXE model for hadrons
175 << theParameters->PIXECrossSectionM
176 << G4endl;
177 G4cout << "### === PIXE model for e+-:
178 << theParameters->PIXEElectronCross
179 << G4endl;
180 }
181 }
182
183 //....oooOO0OOooo........oooOO0OOooo........oo
184
185 void
186 G4VAtomDeexcitation::SetDeexcitationActiveRegi
187
188
189
190 {
191 // no PIXE in parallel world
192 if(rname == "DefaultRegionForParallelWorld")
193
194 G4String ss = rname;
195 /*
196 G4cout << "### G4VAtomDeexcitation::SetDeexc
197 << " " << valDeexcitation << " " <<
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
206 // Region already exist
207 if(ss == activeRegions[i]) {
208 deRegions[i] = valDeexcitation;
209 AugerRegions[i] = valAuger;
210 PIXERegions[i] = valPIXE;
211 return;
212 }
213 }
214 // New region
215 activeRegions.push_back(ss);
216 deRegions.push_back(valDeexcitation);
217 AugerRegions.push_back(valAuger);
218 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 }
256
257 //....oooOO0OOooo........oooOO0OOooo........oo
258
259 void
260 G4VAtomDeexcitation::AlongStepDeexcitation(std
261 con
262 G4d
263 G4i
264 {
265 G4double truelength = step.GetStepLength();
266 if(!flagPIXE && !activePIXEMedia[coupleIndex
267 if(eLossMax <= 0.0 || truelength <= 0.0)
268
269 // step parameters
270 const G4StepPoint* preStep = step.GetPreStep
271 const G4ThreeVector prePos = preStep->GetPos
272 const G4ThreeVector delta = step.GetPostStep
273 const G4double preTime = preStep->GetGlobalT
274 const G4double dt = step.GetPostStepPoint()-
275
276 // particle parameters
277 const G4Track* track = step.GetTrack();
278 const G4ParticleDefinition* part = track->Ge
279 G4double ekin = preStep->GetKineticEnergy();
280
281 // media parameters
282 G4double gCut = (*theCoupleTable->GetEnergyC
283 if(ignoreCuts) { gCut = 0.0; }
284 G4double eCut = DBL_MAX;
285 if(CheckAugerActiveRegion(coupleIndex)) {
286 eCut = (*theCoupleTable->GetEnergyCutsVect
287 if(ignoreCuts) { eCut = 0.0; }
288 }
289
290 //G4cout<<"!Sample PIXE gCut(MeV)= "<<gCut<<
291 // <<" Ekin(MeV)= " << ekin/MeV << G4
292
293 const G4Material* material = preStep->GetMat
294 const G4ElementVector* theElementVector = ma
295 const G4double* theAtomNumDensityVector =
296 material->GetVecNbOfAtomsPerVolume();
297 const std::size_t nelm = material->GetNumber
298
299 // loop over deexcitations
300 for(std::size_t i=0; i<nelm; ++i) {
301 G4int Z = (*theElementVector)[i]->GetZasIn
302 if(activeZ[Z] && Z < 93) {
303 G4int nshells =
304 std::min(9,(*theElementVector)[i]->Get
305 G4double rho = truelength*theAtomNumDens
306 //G4cout<<" Z "<< Z <<" is active x(m
307 for(G4int ii=0; ii<nshells; ++ii) {
308 auto as = (G4AtomicShellEnumerator)(ii
309 const G4AtomicShell* shell = GetAtomic
310 const G4double bindingEnergy = shell->
311
312 if(gCut > bindingEnergy) { break; }
313
314 if(eLossMax > bindingEnergy) {
315 G4double sig = rho*
316 GetShellIonisationCrossSectionPerA
317
318 // mfp is mean free path in units of
319 if(sig > 0.0) {
320 G4double mfp = 1.0/sig;
321 G4double stot = 0.0;
322 //G4cout << " Shell " << ii << " m
323 // sample ionisation points
324 do {
325 stot -= mfp*G4Log(G4UniformRand(
326 if( stot > 1.0 || eLossMax < bin
327 // sample deexcitation
328 vdyn.clear();
329 GenerateParticles(&vdyn, shell,
330 std::size_t nsec = vdyn.size();
331 if(nsec > 0) {
332 G4ThreeVector r = prePos + st
333 G4double time = preTime + st
334 for(std::size_t j=0; j<nsec; +
335 G4DynamicParticle* dp = vdyn
336 G4double e = dp->GetKineticE
337
338 // save new secondary if the
339 if(eLossMax >= e) {
340 eLossMax -= e;
341 G4Track* t = new G4Track(d
342
343 // defined secondary type
344 if(dp->GetDefinition() ==
345 t->SetCreatorModelID(_Ga
346 } else {
347 t->SetCreatorModelID(_eP
348 }
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 }
360 }
361 }
362 return;
363 }
364
365 //....oooOO0OOooo........oooOO0OOooo........oo
366