1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // $Id: G4Penelope08IonisationModel.cc,v 1.2 2010/12/15 07:39:12 gunter Exp $ |
---|
27 | // GEANT4 tag $Name: geant4-09-04-ref-00 $ |
---|
28 | // |
---|
29 | // Author: Luciano Pandola |
---|
30 | // |
---|
31 | // History: |
---|
32 | // -------- |
---|
33 | // 27 Jul 2010 L Pandola First complete implementation |
---|
34 | // |
---|
35 | |
---|
36 | #include "G4Penelope08IonisationModel.hh" |
---|
37 | #include "G4ParticleDefinition.hh" |
---|
38 | #include "G4MaterialCutsCouple.hh" |
---|
39 | #include "G4ProductionCutsTable.hh" |
---|
40 | #include "G4DynamicParticle.hh" |
---|
41 | #include "G4AtomicTransitionManager.hh" |
---|
42 | #include "G4AtomicDeexcitation.hh" |
---|
43 | #include "G4AtomicShell.hh" |
---|
44 | #include "G4Gamma.hh" |
---|
45 | #include "G4Electron.hh" |
---|
46 | #include "G4Positron.hh" |
---|
47 | #include "G4AtomicDeexcitation.hh" |
---|
48 | #include "G4PenelopeOscillatorManager.hh" |
---|
49 | #include "G4PenelopeOscillator.hh" |
---|
50 | #include "G4PenelopeCrossSection.hh" |
---|
51 | #include "G4PhysicsFreeVector.hh" |
---|
52 | #include "G4PhysicsLogVector.hh" |
---|
53 | |
---|
54 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
55 | |
---|
56 | |
---|
57 | G4Penelope08IonisationModel::G4Penelope08IonisationModel(const G4ParticleDefinition*, |
---|
58 | const G4String& nam) |
---|
59 | :G4VEmModel(nam),isInitialised(false), |
---|
60 | kineticEnergy1(0.*eV), |
---|
61 | cosThetaPrimary(1.0), |
---|
62 | energySecondary(0.*eV), |
---|
63 | cosThetaSecondary(0.0), |
---|
64 | targetOscillator(-1),XSTableElectron(0),XSTablePositron(0), |
---|
65 | theDeltaTable(0),energyGrid(0) |
---|
66 | { |
---|
67 | fIntrinsicLowEnergyLimit = 100.0*eV; |
---|
68 | fIntrinsicHighEnergyLimit = 100.0*GeV; |
---|
69 | // SetLowEnergyLimit(fIntrinsicLowEnergyLimit); |
---|
70 | SetHighEnergyLimit(fIntrinsicHighEnergyLimit); |
---|
71 | nBins = 200; |
---|
72 | // |
---|
73 | oscManager = G4PenelopeOscillatorManager::GetOscillatorManager(); |
---|
74 | // |
---|
75 | verboseLevel= 0; |
---|
76 | |
---|
77 | // Verbosity scale: |
---|
78 | // 0 = nothing |
---|
79 | // 1 = warning for energy non-conservation |
---|
80 | // 2 = details of energy budget |
---|
81 | // 3 = calculation of cross sections, file openings, sampling of atoms |
---|
82 | // 4 = entering in methods |
---|
83 | |
---|
84 | // Atomic deexcitation model activated by default |
---|
85 | SetDeexcitationFlag(true); |
---|
86 | ActivateAuger(false); |
---|
87 | } |
---|
88 | |
---|
89 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
90 | |
---|
91 | G4Penelope08IonisationModel::~G4Penelope08IonisationModel() |
---|
92 | { |
---|
93 | ClearTables(); |
---|
94 | } |
---|
95 | |
---|
96 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
97 | void G4Penelope08IonisationModel::Initialise(const G4ParticleDefinition*, |
---|
98 | const G4DataVector&) |
---|
99 | { |
---|
100 | if (verboseLevel > 3) |
---|
101 | G4cout << "Calling G4Penelope08IonisationModel::Initialise()" << G4endl; |
---|
102 | |
---|
103 | //Clear and re-build the tables |
---|
104 | ClearTables(); |
---|
105 | XSTableElectron = new |
---|
106 | std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>; |
---|
107 | XSTablePositron = new |
---|
108 | std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>; |
---|
109 | |
---|
110 | theDeltaTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; |
---|
111 | |
---|
112 | //Set the number of bins for the tables. 20 points per decade |
---|
113 | nBins = (size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit())); |
---|
114 | nBins = std::max(nBins,(size_t)100); |
---|
115 | |
---|
116 | energyGrid = new G4PhysicsLogVector(LowEnergyLimit(), |
---|
117 | HighEnergyLimit(), |
---|
118 | nBins-1); //one hidden bin is added |
---|
119 | |
---|
120 | if (verboseLevel > 2) { |
---|
121 | G4cout << "Penelope Ionisation model is initialized " << G4endl |
---|
122 | << "Energy range: " |
---|
123 | << LowEnergyLimit() / keV << " keV - " |
---|
124 | << HighEnergyLimit() / GeV << " GeV. Using " |
---|
125 | << nBins << " bins." |
---|
126 | << G4endl; |
---|
127 | } |
---|
128 | |
---|
129 | if(isInitialised) return; |
---|
130 | fParticleChange = GetParticleChangeForLoss(); |
---|
131 | isInitialised = true; |
---|
132 | } |
---|
133 | |
---|
134 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
135 | |
---|
136 | G4double G4Penelope08IonisationModel::CrossSectionPerVolume(const G4Material* material, |
---|
137 | const G4ParticleDefinition* theParticle, |
---|
138 | G4double energy, |
---|
139 | G4double cutEnergy, |
---|
140 | G4double) |
---|
141 | { |
---|
142 | // Penelope model to calculate the cross section for inelastic collisions above the |
---|
143 | // threshold. It makes use of the Generalised Oscillator Strength (GOS) model from |
---|
144 | // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567 |
---|
145 | // |
---|
146 | // The total cross section is calculated analytically by taking |
---|
147 | // into account the atomic oscillators coming into the play for a given threshold. |
---|
148 | // |
---|
149 | // For incident e- the maximum energy allowed for the delta-rays is energy/2. |
---|
150 | // because particles are undistinghishable. |
---|
151 | // |
---|
152 | // The contribution is splitted in three parts: distant longitudinal collisions, |
---|
153 | // distant transverse collisions and close collisions. Each term is described by |
---|
154 | // its own analytical function. |
---|
155 | // Fermi density correction is calculated analytically according to |
---|
156 | // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1 |
---|
157 | // |
---|
158 | if (verboseLevel > 3) |
---|
159 | G4cout << "Calling CrossSectionPerVolume() of G4Penelope08IonisationModel" << G4endl; |
---|
160 | |
---|
161 | SetupForMaterial(theParticle, material, energy); |
---|
162 | |
---|
163 | G4double totalCross = 0.0; |
---|
164 | G4double crossPerMolecule = 0.; |
---|
165 | |
---|
166 | G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material, |
---|
167 | cutEnergy); |
---|
168 | |
---|
169 | if (theXS) |
---|
170 | crossPerMolecule = theXS->GetHardCrossSection(energy); |
---|
171 | |
---|
172 | G4double atomDensity = material->GetTotNbOfAtomsPerVolume(); |
---|
173 | G4double atPerMol = oscManager->GetAtomsPerMolecule(material); |
---|
174 | |
---|
175 | if (verboseLevel > 3) |
---|
176 | G4cout << "Material " << material->GetName() << " has " << atPerMol << |
---|
177 | "atoms per molecule" << G4endl; |
---|
178 | |
---|
179 | G4double moleculeDensity = 0.; |
---|
180 | if (atPerMol) |
---|
181 | moleculeDensity = atomDensity/atPerMol; |
---|
182 | |
---|
183 | G4double crossPerVolume = crossPerMolecule*moleculeDensity; |
---|
184 | |
---|
185 | if (verboseLevel > 2) |
---|
186 | { |
---|
187 | G4cout << "G4Penelope08IonisationModel " << G4endl; |
---|
188 | G4cout << "Mean free path for delta emission > " << cutEnergy/keV << " keV at " << |
---|
189 | energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl; |
---|
190 | if (theXS) |
---|
191 | totalCross = (theXS->GetTotalCrossSection(energy))*moleculeDensity; |
---|
192 | G4cout << "Total free path for ionisation (no threshold) at " << |
---|
193 | energy/keV << " keV = " << (1./totalCross)/mm << " mm" << G4endl; |
---|
194 | } |
---|
195 | return crossPerVolume; |
---|
196 | } |
---|
197 | |
---|
198 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
199 | |
---|
200 | //This is a dummy method. Never inkoved by the tracking, it just issues |
---|
201 | //a warning if one tries to get Cross Sections per Atom via the |
---|
202 | //G4EmCalculator. |
---|
203 | G4double G4Penelope08IonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, |
---|
204 | G4double, |
---|
205 | G4double, |
---|
206 | G4double, |
---|
207 | G4double, |
---|
208 | G4double) |
---|
209 | { |
---|
210 | G4cout << "*** G4Penelope08IonisationModel -- WARNING ***" << G4endl; |
---|
211 | G4cout << "Penelope08 Ionisation model does not calculate cross section _per atom_ " << G4endl; |
---|
212 | G4cout << "so the result is always zero. For physics values, please invoke " << G4endl; |
---|
213 | G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl; |
---|
214 | return 0; |
---|
215 | } |
---|
216 | |
---|
217 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
218 | |
---|
219 | G4double G4Penelope08IonisationModel::ComputeDEDXPerVolume(const G4Material* material, |
---|
220 | const G4ParticleDefinition* theParticle, |
---|
221 | G4double kineticEnergy, |
---|
222 | G4double cutEnergy) |
---|
223 | { |
---|
224 | // Penelope model to calculate the stopping power for soft inelastic collisions |
---|
225 | // below the threshold. It makes use of the Generalised Oscillator Strength (GOS) |
---|
226 | // model from |
---|
227 | // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567 |
---|
228 | // |
---|
229 | // The stopping power is calculated analytically using the dsigma/dW cross |
---|
230 | // section from the GOS models, which includes separate contributions from |
---|
231 | // distant longitudinal collisions, distant transverse collisions and |
---|
232 | // close collisions. Only the atomic oscillators that come in the play |
---|
233 | // (according to the threshold) are considered for the calculation. |
---|
234 | // Differential cross sections have a different form for e+ and e-. |
---|
235 | // |
---|
236 | // Fermi density correction is calculated analytically according to |
---|
237 | // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1 |
---|
238 | |
---|
239 | if (verboseLevel > 3) |
---|
240 | G4cout << "Calling ComputeDEDX() of G4Penelope08IonisationModel" << G4endl; |
---|
241 | |
---|
242 | |
---|
243 | G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material, |
---|
244 | cutEnergy); |
---|
245 | G4double sPowerPerMolecule = 0.0; |
---|
246 | if (theXS) |
---|
247 | sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy); |
---|
248 | |
---|
249 | |
---|
250 | G4double atomDensity = material->GetTotNbOfAtomsPerVolume(); |
---|
251 | G4double atPerMol = oscManager->GetAtomsPerMolecule(material); |
---|
252 | |
---|
253 | G4double moleculeDensity = 0.; |
---|
254 | if (atPerMol) |
---|
255 | moleculeDensity = atomDensity/atPerMol; |
---|
256 | |
---|
257 | G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity; |
---|
258 | |
---|
259 | if (verboseLevel > 2) |
---|
260 | { |
---|
261 | G4cout << "G4Penelope08IonisationModel " << G4endl; |
---|
262 | G4cout << "Stopping power < " << cutEnergy/keV << " keV at " << |
---|
263 | kineticEnergy/keV << " keV = " << |
---|
264 | sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl; |
---|
265 | } |
---|
266 | return sPowerPerVolume; |
---|
267 | } |
---|
268 | |
---|
269 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
270 | |
---|
271 | G4double G4Penelope08IonisationModel::MinEnergyCut(const G4ParticleDefinition*, |
---|
272 | const G4MaterialCutsCouple*) |
---|
273 | { |
---|
274 | return fIntrinsicLowEnergyLimit; |
---|
275 | } |
---|
276 | |
---|
277 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
278 | |
---|
279 | void G4Penelope08IonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, |
---|
280 | const G4MaterialCutsCouple* couple, |
---|
281 | const G4DynamicParticle* aDynamicParticle, |
---|
282 | G4double cutE, G4double) |
---|
283 | { |
---|
284 | // Penelope model to sample the final state following an hard inelastic interaction. |
---|
285 | // It makes use of the Generalised Oscillator Strength (GOS) model from |
---|
286 | // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567 |
---|
287 | // |
---|
288 | // The GOS model is used to calculate the individual cross sections for all |
---|
289 | // the atomic oscillators coming in the play, taking into account the three |
---|
290 | // contributions (distant longitudinal collisions, distant transverse collisions and |
---|
291 | // close collisions). Then the target shell and the interaction channel are |
---|
292 | // sampled. Final state of the delta-ray (energy, angle) are generated according |
---|
293 | // to the analytical distributions (dSigma/dW) for the selected interaction |
---|
294 | // channels. |
---|
295 | // For e-, the maximum energy for the delta-ray is initialEnergy/2. (because |
---|
296 | // particles are indistinghusbable), while it is the full initialEnergy for |
---|
297 | // e+. |
---|
298 | // The efficiency of the random sampling algorithm (e.g. for close collisions) |
---|
299 | // decreases when initial and cutoff energy increase (e.g. 87% for 10-keV primary |
---|
300 | // and 1 keV threshold, 99% for 10-MeV primary and 10-keV threshold). |
---|
301 | // Differential cross sections have a different form for e+ and e-. |
---|
302 | // |
---|
303 | // WARNING: The model provides an _average_ description of inelastic collisions. |
---|
304 | // Anyway, the energy spectrum associated to distant excitations of a given |
---|
305 | // atomic shell is approximated as a single resonance. The simulated energy spectra |
---|
306 | // show _unphysical_ narrow peaks at energies that are multiple of the shell |
---|
307 | // resonance energies. The spurious speaks are automatically smoothed out after |
---|
308 | // multiple inelastic collisions. |
---|
309 | // |
---|
310 | // The model determines also the original shell from which the delta-ray is expelled, |
---|
311 | // in order to produce fluorescence de-excitation (from G4DeexcitationManager) |
---|
312 | // |
---|
313 | // Fermi density correction is calculated analytically according to |
---|
314 | // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1 |
---|
315 | |
---|
316 | if (verboseLevel > 3) |
---|
317 | G4cout << "Calling SamplingSecondaries() of G4Penelope08IonisationModel" << G4endl; |
---|
318 | |
---|
319 | G4double kineticEnergy0 = aDynamicParticle->GetKineticEnergy(); |
---|
320 | const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition(); |
---|
321 | |
---|
322 | if (kineticEnergy0 <= fIntrinsicLowEnergyLimit) |
---|
323 | { |
---|
324 | fParticleChange->SetProposedKineticEnergy(0.); |
---|
325 | fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy0); |
---|
326 | return ; |
---|
327 | } |
---|
328 | const G4Material* material = couple->GetMaterial(); |
---|
329 | G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(material); |
---|
330 | |
---|
331 | G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection(); |
---|
332 | |
---|
333 | //Initialise final-state variables. The proper values will be set by the methods |
---|
334 | // SampleFinalStateElectron() and SampleFinalStatePositron() |
---|
335 | kineticEnergy1=kineticEnergy0; |
---|
336 | cosThetaPrimary=1.0; |
---|
337 | energySecondary=0.0; |
---|
338 | cosThetaSecondary=1.0; |
---|
339 | targetOscillator = -1; |
---|
340 | |
---|
341 | if (theParticle == G4Electron::Electron()) |
---|
342 | SampleFinalStateElectron(material,cutE,kineticEnergy0); |
---|
343 | else if (theParticle == G4Positron::Positron()) |
---|
344 | SampleFinalStatePositron(material,cutE,kineticEnergy0); |
---|
345 | else |
---|
346 | { |
---|
347 | G4cout << "G4Penelope08IonisationModel::SamplingSecondaries() - " ; |
---|
348 | G4cout << "Invalid particle " << theParticle->GetParticleName() << G4endl; |
---|
349 | G4Exception(); |
---|
350 | } |
---|
351 | |
---|
352 | if (verboseLevel > 3) |
---|
353 | { |
---|
354 | G4cout << "G4Penelope08IonisationModel::SamplingSecondaries() for " << |
---|
355 | theParticle->GetParticleName() << G4endl; |
---|
356 | G4cout << "Final eKin = " << kineticEnergy1 << " keV" << G4endl; |
---|
357 | G4cout << "Final cosTheta = " << cosThetaPrimary << G4endl; |
---|
358 | G4cout << "Delta-ray eKin = " << energySecondary << " keV" << G4endl; |
---|
359 | G4cout << "Delta-ray cosTheta = " << cosThetaSecondary << G4endl; |
---|
360 | G4cout << "Oscillator: " << targetOscillator << G4endl; |
---|
361 | } |
---|
362 | |
---|
363 | //Update the primary particle |
---|
364 | G4double sint = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary); |
---|
365 | G4double phiPrimary = twopi * G4UniformRand(); |
---|
366 | G4double dirx = sint * std::cos(phiPrimary); |
---|
367 | G4double diry = sint * std::sin(phiPrimary); |
---|
368 | G4double dirz = cosThetaPrimary; |
---|
369 | |
---|
370 | G4ThreeVector electronDirection1(dirx,diry,dirz); |
---|
371 | electronDirection1.rotateUz(particleDirection0); |
---|
372 | |
---|
373 | if (kineticEnergy1 > 0) |
---|
374 | { |
---|
375 | fParticleChange->ProposeMomentumDirection(electronDirection1); |
---|
376 | fParticleChange->SetProposedKineticEnergy(kineticEnergy1); |
---|
377 | } |
---|
378 | else |
---|
379 | { |
---|
380 | fParticleChange->SetProposedKineticEnergy(0.); |
---|
381 | } |
---|
382 | |
---|
383 | //Generate the delta ray |
---|
384 | G4double ionEnergyInPenelopeDatabase = |
---|
385 | (*theTable)[targetOscillator]->GetIonisationEnergy(); |
---|
386 | |
---|
387 | //Now, try to handle fluorescence |
---|
388 | //Notice: merged levels are indicated with Z=0 and flag=30 |
---|
389 | G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag(); |
---|
390 | G4int Z = (G4int) (*theTable)[targetOscillator]->GetParentZ(); |
---|
391 | |
---|
392 | //initialize here, then check photons created by Atomic-Deexcitation, and the final state e- |
---|
393 | std::vector<G4DynamicParticle*>* photonVector=0; |
---|
394 | const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance(); |
---|
395 | G4double bindingEnergy = 0.*eV; |
---|
396 | G4int shellId = 0; |
---|
397 | |
---|
398 | //Real level |
---|
399 | if (Z > 0 && shFlag<30) |
---|
400 | { |
---|
401 | const G4AtomicShell* shell = transitionManager->Shell(Z,shFlag-1); |
---|
402 | bindingEnergy = shell->BindingEnergy(); |
---|
403 | shellId = shell->ShellId(); |
---|
404 | } |
---|
405 | |
---|
406 | //correct the energySecondary to account for the fact that the Penelope |
---|
407 | //database of ionisation energies is in general (slightly) different |
---|
408 | //from the fluorescence database used in Geant4. |
---|
409 | energySecondary += ionEnergyInPenelopeDatabase-bindingEnergy; |
---|
410 | |
---|
411 | G4double localEnergyDeposit = bindingEnergy; |
---|
412 | //testing purposes only |
---|
413 | G4double energyInFluorescence = 0; |
---|
414 | |
---|
415 | if (energySecondary < 0) |
---|
416 | { |
---|
417 | //It means that there was some problem/mismatch between the two databases. |
---|
418 | //In this case, the available energy is ok to excite the level according |
---|
419 | //to the Penelope database, but not according to the Geant4 database |
---|
420 | //Full residual energy is deposited locally |
---|
421 | localEnergyDeposit += energySecondary; |
---|
422 | energySecondary = 0.0; |
---|
423 | } |
---|
424 | |
---|
425 | //the local energy deposit is what remains: part of this may be spent for fluorescence. |
---|
426 | if(DeexcitationFlag() && Z > 5) |
---|
427 | { |
---|
428 | const G4ProductionCutsTable* theCoupleTable= |
---|
429 | G4ProductionCutsTable::GetProductionCutsTable(); |
---|
430 | |
---|
431 | size_t index = couple->GetIndex(); |
---|
432 | G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index]; |
---|
433 | G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index]; |
---|
434 | |
---|
435 | // Generation of fluorescence |
---|
436 | // Data in EADL are available only for Z > 5 |
---|
437 | // Protection to avoid generating photons in the unphysical case of |
---|
438 | // shell binding energy > photon energy |
---|
439 | if (localEnergyDeposit > cutg || localEnergyDeposit > cute) |
---|
440 | { |
---|
441 | G4DynamicParticle* aPhoton; |
---|
442 | deexcitationManager.SetCutForSecondaryPhotons(cutg); |
---|
443 | deexcitationManager.SetCutForAugerElectrons(cute); |
---|
444 | |
---|
445 | photonVector = deexcitationManager.GenerateParticles(Z,shellId); |
---|
446 | if(photonVector) |
---|
447 | { |
---|
448 | size_t nPhotons = photonVector->size(); |
---|
449 | for (size_t k=0; k<nPhotons; k++) |
---|
450 | { |
---|
451 | aPhoton = (*photonVector)[k]; |
---|
452 | if (aPhoton) |
---|
453 | { |
---|
454 | G4double itsEnergy = aPhoton->GetKineticEnergy(); |
---|
455 | if (itsEnergy <= localEnergyDeposit) |
---|
456 | { |
---|
457 | localEnergyDeposit -= itsEnergy; |
---|
458 | if (aPhoton->GetDefinition() == G4Gamma::Gamma()) |
---|
459 | energyInFluorescence += itsEnergy;; |
---|
460 | fvect->push_back(aPhoton); |
---|
461 | } |
---|
462 | else |
---|
463 | { |
---|
464 | delete aPhoton; |
---|
465 | (*photonVector)[k]=0; |
---|
466 | } |
---|
467 | } |
---|
468 | } |
---|
469 | delete photonVector; |
---|
470 | } |
---|
471 | } |
---|
472 | } |
---|
473 | |
---|
474 | //Always produce explicitely the delta electron |
---|
475 | G4DynamicParticle* electron = 0; |
---|
476 | G4double sinThetaE = std::sqrt(1.-cosThetaSecondary*cosThetaSecondary); |
---|
477 | G4double phiEl = phiPrimary+pi; //pi with respect to the primary electron/positron |
---|
478 | G4double xEl = sinThetaE * std::cos(phiEl); |
---|
479 | G4double yEl = sinThetaE * std::sin(phiEl); |
---|
480 | G4double zEl = cosThetaSecondary; |
---|
481 | G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction |
---|
482 | eDirection.rotateUz(particleDirection0); |
---|
483 | electron = new G4DynamicParticle (G4Electron::Electron(), |
---|
484 | eDirection,energySecondary) ; |
---|
485 | fvect->push_back(electron); |
---|
486 | |
---|
487 | |
---|
488 | if (localEnergyDeposit < 0) |
---|
489 | { |
---|
490 | G4cout << "WARNING-" |
---|
491 | << "G4Penelope08IonisationModel::SampleSecondaries - Negative energy deposit" |
---|
492 | << G4endl; |
---|
493 | localEnergyDeposit=0.; |
---|
494 | } |
---|
495 | fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit); |
---|
496 | |
---|
497 | if (verboseLevel > 1) |
---|
498 | { |
---|
499 | G4cout << "-----------------------------------------------------------" << G4endl; |
---|
500 | G4cout << "Energy balance from G4Penelope08Ionisation" << G4endl; |
---|
501 | G4cout << "Incoming primary energy: " << kineticEnergy0/keV << " keV" << G4endl; |
---|
502 | G4cout << "-----------------------------------------------------------" << G4endl; |
---|
503 | G4cout << "Outgoing primary energy: " << kineticEnergy1/keV << " keV" << G4endl; |
---|
504 | G4cout << "Delta ray " << energySecondary/keV << " keV" << G4endl; |
---|
505 | G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl; |
---|
506 | G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl; |
---|
507 | G4cout << "Total final state: " << (energySecondary+energyInFluorescence+kineticEnergy1+ |
---|
508 | localEnergyDeposit)/keV << |
---|
509 | " keV" << G4endl; |
---|
510 | G4cout << "-----------------------------------------------------------" << G4endl; |
---|
511 | } |
---|
512 | |
---|
513 | if (verboseLevel > 0) |
---|
514 | { |
---|
515 | G4double energyDiff = std::fabs(energySecondary+energyInFluorescence+kineticEnergy1+ |
---|
516 | localEnergyDeposit-kineticEnergy0); |
---|
517 | if (energyDiff > 0.05*keV) |
---|
518 | G4cout << "Warning from G4PenelopeIonisation: problem with energy conservation: " << |
---|
519 | (energySecondary+energyInFluorescence+kineticEnergy1+localEnergyDeposit)/keV << |
---|
520 | " keV (final) vs. " << |
---|
521 | kineticEnergy0/keV << " keV (initial)" << G4endl; |
---|
522 | } |
---|
523 | |
---|
524 | } |
---|
525 | |
---|
526 | |
---|
527 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
528 | |
---|
529 | void G4Penelope08IonisationModel::ActivateAuger(G4bool augerbool) |
---|
530 | { |
---|
531 | if (!DeexcitationFlag() && augerbool) |
---|
532 | { |
---|
533 | G4cout << "WARNING - G4Penelope08IonisationModel" << G4endl; |
---|
534 | G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl; |
---|
535 | G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl; |
---|
536 | } |
---|
537 | deexcitationManager.ActivateAugerElectronProduction(augerbool); |
---|
538 | if (verboseLevel > 1) |
---|
539 | G4cout << "Auger production set to " << augerbool << G4endl; |
---|
540 | } |
---|
541 | |
---|
542 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
543 | |
---|
544 | void G4Penelope08IonisationModel::ClearTables() |
---|
545 | { |
---|
546 | std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>::iterator i; |
---|
547 | if (XSTableElectron) |
---|
548 | { |
---|
549 | for (i=XSTableElectron->begin(); i != XSTableElectron->end(); i++) |
---|
550 | { |
---|
551 | G4PenelopeCrossSection* tab = i->second; |
---|
552 | delete tab; |
---|
553 | } |
---|
554 | delete XSTableElectron; |
---|
555 | XSTableElectron = 0; |
---|
556 | } |
---|
557 | |
---|
558 | if (XSTablePositron) |
---|
559 | { |
---|
560 | for (i=XSTablePositron->begin(); i != XSTablePositron->end(); i++) |
---|
561 | { |
---|
562 | G4PenelopeCrossSection* tab = i->second; |
---|
563 | delete tab; |
---|
564 | } |
---|
565 | delete XSTablePositron; |
---|
566 | XSTablePositron = 0; |
---|
567 | } |
---|
568 | |
---|
569 | std::map<const G4Material*,G4PhysicsFreeVector*>::iterator k; |
---|
570 | if (theDeltaTable) |
---|
571 | { |
---|
572 | for (k=theDeltaTable->begin();k!=theDeltaTable->end();k++) |
---|
573 | delete k->second; |
---|
574 | delete theDeltaTable; |
---|
575 | theDeltaTable = 0; |
---|
576 | } |
---|
577 | |
---|
578 | if (energyGrid) |
---|
579 | delete energyGrid; |
---|
580 | |
---|
581 | if (verboseLevel > 2) |
---|
582 | G4cout << "G4Penelope08IonisationModel: cleared tables" << G4endl; |
---|
583 | return; |
---|
584 | } |
---|
585 | |
---|
586 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
587 | |
---|
588 | G4PenelopeCrossSection* |
---|
589 | G4Penelope08IonisationModel::GetCrossSectionTableForCouple(const G4ParticleDefinition* part, |
---|
590 | const G4Material* mat, |
---|
591 | G4double cut) |
---|
592 | { |
---|
593 | if (part != G4Electron::Electron() && part != G4Positron::Positron()) |
---|
594 | { |
---|
595 | G4cout << "G4Penelope08IonisationModel::GetCrossSectionTableForCouple" << G4endl; |
---|
596 | G4cout << "Invalid particle: " << part->GetParticleName() << G4endl; |
---|
597 | G4Exception(); |
---|
598 | return NULL; |
---|
599 | } |
---|
600 | |
---|
601 | if (part == G4Electron::Electron()) |
---|
602 | { |
---|
603 | if (!XSTableElectron) |
---|
604 | { |
---|
605 | G4cout << "G4Penelope08IonisationModel::GetCrossSectionTableForCouple()" << G4endl; |
---|
606 | G4cout << "The Cross Section Table for e- was not initialized correctly!" << G4endl; |
---|
607 | G4Exception(); |
---|
608 | } |
---|
609 | std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut); |
---|
610 | if (XSTableElectron->count(theKey)) //table already built |
---|
611 | return XSTableElectron->find(theKey)->second; |
---|
612 | else |
---|
613 | { |
---|
614 | BuildXSTable(mat,cut,part); |
---|
615 | if (XSTableElectron->count(theKey)) //now it should be ok! |
---|
616 | return XSTableElectron->find(theKey)->second; |
---|
617 | else |
---|
618 | { |
---|
619 | G4cout << "G4Penelope08IonisationModel::GetCrossSectionTableForCouple()" << G4endl; |
---|
620 | G4cout << "Unable to build e- table for " << mat->GetName() << G4endl; |
---|
621 | G4Exception(); |
---|
622 | } |
---|
623 | } |
---|
624 | } |
---|
625 | |
---|
626 | if (part == G4Positron::Positron()) |
---|
627 | { |
---|
628 | if (!XSTablePositron) |
---|
629 | { |
---|
630 | G4cout << "G4Penelope08IonisationModel::GetCrossSectionTableForCouple()" << G4endl; |
---|
631 | G4cout << "The Cross Section Table for e+ was not initialized correctly!" << G4endl; |
---|
632 | G4Exception(); |
---|
633 | } |
---|
634 | std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut); |
---|
635 | if (XSTablePositron->count(theKey)) //table already built |
---|
636 | return XSTablePositron->find(theKey)->second; |
---|
637 | else |
---|
638 | { |
---|
639 | BuildXSTable(mat,cut,part); |
---|
640 | if (XSTablePositron->count(theKey)) //now it should be ok! |
---|
641 | return XSTablePositron->find(theKey)->second; |
---|
642 | else |
---|
643 | { |
---|
644 | G4cout << "G4Penelope08IonisationModel::GetCrossSectionTableForCouple()" << G4endl; |
---|
645 | G4cout << "Unable to build e+ table for " << mat->GetName() << G4endl; |
---|
646 | G4Exception(); |
---|
647 | } |
---|
648 | } |
---|
649 | } |
---|
650 | return NULL; |
---|
651 | } |
---|
652 | |
---|
653 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
654 | |
---|
655 | void G4Penelope08IonisationModel::BuildXSTable(const G4Material* mat,G4double cut, |
---|
656 | const G4ParticleDefinition* part) |
---|
657 | { |
---|
658 | // |
---|
659 | //This method fills the G4PenelopeCrossSection containers for electrons or positrons |
---|
660 | //and for the given material/cut couple. The calculation is done as sum over the |
---|
661 | //individual shells. |
---|
662 | //Equivalent of subroutines EINaT and PINaT of Penelope |
---|
663 | // |
---|
664 | if (verboseLevel > 2) |
---|
665 | { |
---|
666 | G4cout << "G4Penelope08IonisationModel: going to build cross section table " << G4endl; |
---|
667 | G4cout << "for " << part->GetParticleName() << " in " << mat->GetName() << G4endl; |
---|
668 | } |
---|
669 | |
---|
670 | //Tables have been already created (checked by GetCrossSectionTableForCouple) |
---|
671 | G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat); |
---|
672 | size_t numberOfOscillators = theTable->size(); |
---|
673 | |
---|
674 | if (energyGrid->GetVectorLength() != nBins) |
---|
675 | { |
---|
676 | G4cout << "G4Penelope08IonisationModel::BuildXSTable" << G4endl; |
---|
677 | G4cout << "Energy Grid looks not initialized" << G4endl; |
---|
678 | G4cout << nBins << " " << energyGrid->GetVectorLength() << G4endl; |
---|
679 | G4Exception(); |
---|
680 | } |
---|
681 | |
---|
682 | G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(nBins,numberOfOscillators); |
---|
683 | |
---|
684 | //loop on the energy grid |
---|
685 | for (size_t bin=0;bin<nBins;bin++) |
---|
686 | { |
---|
687 | G4double energy = energyGrid->GetLowEdgeEnergy(bin); |
---|
688 | G4double XH0=0, XH1=0, XH2=0; |
---|
689 | G4double XS0=0, XS1=0, XS2=0; |
---|
690 | |
---|
691 | //oscillator loop |
---|
692 | for (size_t iosc=0;iosc<numberOfOscillators;iosc++) |
---|
693 | { |
---|
694 | G4DataVector* tempStorage = 0; |
---|
695 | |
---|
696 | G4PenelopeOscillator* theOsc = (*theTable)[iosc]; |
---|
697 | G4double delta = GetDensityCorrection(mat,energy); |
---|
698 | if (part == G4Electron::Electron()) |
---|
699 | tempStorage = ComputeShellCrossSectionsElectron(theOsc,energy,cut,delta); |
---|
700 | else if (part == G4Positron::Positron()) |
---|
701 | tempStorage = ComputeShellCrossSectionsPositron(theOsc,energy,cut,delta); |
---|
702 | //check results are all right |
---|
703 | if (!tempStorage) |
---|
704 | { |
---|
705 | G4cout << "G4Penelope08IonisationModel::BuildXSTable" << G4endl; |
---|
706 | G4cout << "Problem in calculating the shell XS " << G4endl; |
---|
707 | G4Exception(); |
---|
708 | } |
---|
709 | if (tempStorage->size() != 6) |
---|
710 | { |
---|
711 | G4cout << "G4Penelope08IonisationModel::BuildXSTable" << G4endl; |
---|
712 | G4cout << "Problem in calculating the shell XS " << G4endl; |
---|
713 | G4cout << "Result has dimension " << tempStorage->size() << " instead of 6" << G4endl; |
---|
714 | G4Exception(); |
---|
715 | } |
---|
716 | G4double stre = theOsc->GetOscillatorStrength(); |
---|
717 | |
---|
718 | XH0 += stre*(*tempStorage)[0]; |
---|
719 | XH1 += stre*(*tempStorage)[1]; |
---|
720 | XH2 += stre*(*tempStorage)[2]; |
---|
721 | XS0 += stre*(*tempStorage)[3]; |
---|
722 | XS1 += stre*(*tempStorage)[4]; |
---|
723 | XS2 += stre*(*tempStorage)[5]; |
---|
724 | XSEntry->AddShellCrossSectionPoint(bin,iosc,energy,stre*(*tempStorage)[0]); |
---|
725 | if (tempStorage) |
---|
726 | { |
---|
727 | delete tempStorage; |
---|
728 | tempStorage = 0; |
---|
729 | } |
---|
730 | } |
---|
731 | XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2); |
---|
732 | } |
---|
733 | |
---|
734 | |
---|
735 | //Normalize shell cross sections before insertion in the table |
---|
736 | XSEntry->NormalizeShellCrossSections(); |
---|
737 | |
---|
738 | |
---|
739 | //Insert in the appropriate table |
---|
740 | std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut); |
---|
741 | if (part == G4Electron::Electron()) |
---|
742 | XSTableElectron->insert(std::make_pair(theKey,XSEntry)); |
---|
743 | else if (part == G4Positron::Positron()) |
---|
744 | XSTablePositron->insert(std::make_pair(theKey,XSEntry)); |
---|
745 | |
---|
746 | return; |
---|
747 | } |
---|
748 | |
---|
749 | |
---|
750 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
751 | |
---|
752 | G4double G4Penelope08IonisationModel::GetDensityCorrection(const G4Material* mat, |
---|
753 | G4double energy) |
---|
754 | { |
---|
755 | G4double result = 0; |
---|
756 | if (!theDeltaTable) |
---|
757 | { |
---|
758 | G4cout << "G4Penelope08IonisationModel::GetDensityCorrection()" << G4endl; |
---|
759 | G4cout << "Delta Table not initialized. Was Initialise() run?" << G4endl; |
---|
760 | G4Exception(); |
---|
761 | } |
---|
762 | if (energy <= 0*eV) |
---|
763 | { |
---|
764 | G4cout << "G4Penelope08IonisationModel::GetDensityCorrection()" << G4endl; |
---|
765 | G4cout << "Invalid energy " << energy/eV << " eV " << G4endl; |
---|
766 | return 0; |
---|
767 | } |
---|
768 | G4double logene = std::log(energy); |
---|
769 | |
---|
770 | //check if the material has been built |
---|
771 | if (!(theDeltaTable->count(mat))) |
---|
772 | BuildDeltaTable(mat); |
---|
773 | |
---|
774 | if (theDeltaTable->count(mat)) |
---|
775 | { |
---|
776 | G4PhysicsFreeVector* vec = theDeltaTable->find(mat)->second; |
---|
777 | result = vec->Value(logene); //the table has delta vs. ln(E) |
---|
778 | } |
---|
779 | else |
---|
780 | { |
---|
781 | G4cout << "G4Penelope08IonisationModel::GetDensityCorrection()" << G4endl; |
---|
782 | G4cout << "Unable to build table for " << mat->GetName() << G4endl; |
---|
783 | G4Exception(); |
---|
784 | } |
---|
785 | |
---|
786 | return result; |
---|
787 | } |
---|
788 | |
---|
789 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
790 | |
---|
791 | void G4Penelope08IonisationModel::BuildDeltaTable(const G4Material* mat) |
---|
792 | { |
---|
793 | G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat); |
---|
794 | G4double plasmaSq = oscManager->GetPlasmaEnergySquared(mat); |
---|
795 | G4double totalZ = oscManager->GetTotalZ(mat); |
---|
796 | size_t numberOfOscillators = theTable->size(); |
---|
797 | |
---|
798 | if (energyGrid->GetVectorLength() != nBins) |
---|
799 | { |
---|
800 | G4cout << "G4Penelope08IonisationModel::BuildDeltaTable" << G4endl; |
---|
801 | G4cout << "Energy Grid for Delta looks not initialized" << G4endl; |
---|
802 | G4cout << nBins << " " << energyGrid->GetVectorLength() << G4endl; |
---|
803 | G4Exception(); |
---|
804 | } |
---|
805 | |
---|
806 | G4PhysicsFreeVector* theVector = new G4PhysicsFreeVector(nBins); |
---|
807 | |
---|
808 | //loop on the energy grid |
---|
809 | for (size_t bin=0;bin<nBins;bin++) |
---|
810 | { |
---|
811 | G4double delta = 0.; |
---|
812 | G4double energy = energyGrid->GetLowEdgeEnergy(bin); |
---|
813 | |
---|
814 | //Here calculate delta |
---|
815 | G4double gam = 1.0+(energy/electron_mass_c2); |
---|
816 | G4double gamSq = gam*gam; |
---|
817 | |
---|
818 | G4double TST = totalZ/(gamSq*plasmaSq); |
---|
819 | G4double wl2 = 0; |
---|
820 | G4double fdel = 0; |
---|
821 | |
---|
822 | //loop on oscillators |
---|
823 | for (size_t i=0;i<numberOfOscillators;i++) |
---|
824 | { |
---|
825 | G4PenelopeOscillator* theOsc = (*theTable)[i]; |
---|
826 | G4double wri = theOsc->GetResonanceEnergy(); |
---|
827 | fdel += theOsc->GetOscillatorStrength()/(wri*wri+wl2); |
---|
828 | } |
---|
829 | if (fdel >= TST) //if fdel < TST, delta = 0 |
---|
830 | { |
---|
831 | //get last oscillator |
---|
832 | G4PenelopeOscillator* theOsc = (*theTable)[numberOfOscillators-1]; |
---|
833 | wl2 = theOsc->GetResonanceEnergy()*theOsc->GetResonanceEnergy(); |
---|
834 | |
---|
835 | //First iteration |
---|
836 | G4bool loopAgain = false; |
---|
837 | do |
---|
838 | { |
---|
839 | loopAgain = false; |
---|
840 | wl2 += wl2; |
---|
841 | fdel = 0.; |
---|
842 | for (size_t i=0;i<numberOfOscillators;i++) |
---|
843 | { |
---|
844 | G4PenelopeOscillator* theOsc = (*theTable)[i]; |
---|
845 | G4double wri = theOsc->GetResonanceEnergy(); |
---|
846 | fdel += theOsc->GetOscillatorStrength()/(wri*wri+wl2); |
---|
847 | } |
---|
848 | if (fdel > TST) |
---|
849 | loopAgain = true; |
---|
850 | }while(loopAgain); |
---|
851 | |
---|
852 | G4double wl2l = 0; |
---|
853 | G4double wl2u = wl2; |
---|
854 | //second iteration |
---|
855 | do |
---|
856 | { |
---|
857 | loopAgain = false; |
---|
858 | wl2 = 0.5*(wl2l+wl2u); |
---|
859 | fdel = 0; |
---|
860 | for (size_t i=0;i<numberOfOscillators;i++) |
---|
861 | { |
---|
862 | G4PenelopeOscillator* theOsc = (*theTable)[i]; |
---|
863 | G4double wri = theOsc->GetResonanceEnergy(); |
---|
864 | fdel += theOsc->GetOscillatorStrength()/(wri*wri+wl2); |
---|
865 | } |
---|
866 | if (fdel > TST) |
---|
867 | wl2l = wl2; |
---|
868 | else |
---|
869 | wl2u = wl2; |
---|
870 | if ((wl2u-wl2l)>1e-12*wl2) |
---|
871 | loopAgain = true; |
---|
872 | }while(loopAgain); |
---|
873 | |
---|
874 | //Eventually get density correction |
---|
875 | delta = 0.; |
---|
876 | for (size_t i=0;i<numberOfOscillators;i++) |
---|
877 | { |
---|
878 | G4PenelopeOscillator* theOsc = (*theTable)[i]; |
---|
879 | G4double wri = theOsc->GetResonanceEnergy(); |
---|
880 | delta += theOsc->GetOscillatorStrength()* |
---|
881 | std::log(1.0+(wl2/(wri*wri))); |
---|
882 | } |
---|
883 | delta = (delta/totalZ)-wl2/(gamSq*plasmaSq); |
---|
884 | } |
---|
885 | energy = std::max(1e-9*eV,energy); //prevents log(0) |
---|
886 | theVector->PutValue(bin,std::log(energy),delta); |
---|
887 | } |
---|
888 | theDeltaTable->insert(std::make_pair(mat,theVector)); |
---|
889 | return; |
---|
890 | } |
---|
891 | |
---|
892 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
893 | G4DataVector* G4Penelope08IonisationModel::ComputeShellCrossSectionsElectron(G4PenelopeOscillator* theOsc, |
---|
894 | G4double energy, |
---|
895 | G4double cut, |
---|
896 | G4double delta) |
---|
897 | { |
---|
898 | // |
---|
899 | //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for |
---|
900 | //the given oscillator/cut and at the given energy. |
---|
901 | //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2) |
---|
902 | //Equivalent of subroutines EINaT1 of Penelope |
---|
903 | // |
---|
904 | // Results are _per target electron_ |
---|
905 | // |
---|
906 | G4DataVector* result = new G4DataVector(); |
---|
907 | for (size_t i=0;i<6;i++) |
---|
908 | result->push_back(0.); |
---|
909 | G4double ionEnergy = theOsc->GetIonisationEnergy(); |
---|
910 | |
---|
911 | //return a set of zero's if the energy it too low to excite the current oscillator |
---|
912 | if (energy < ionEnergy) |
---|
913 | return result; |
---|
914 | |
---|
915 | G4double H0=0.,H1=0.,H2=0.; |
---|
916 | G4double S0=0.,S1=0.,S2=0.; |
---|
917 | |
---|
918 | //Define useful constants to be used in the calculation |
---|
919 | G4double gamma = 1.0+energy/electron_mass_c2; |
---|
920 | G4double gammaSq = gamma*gamma; |
---|
921 | G4double beta = (gammaSq-1.0)/gammaSq; |
---|
922 | G4double pielr2 = pi*classic_electr_radius*classic_electr_radius; //pi*re^2 |
---|
923 | G4double constant = pielr2*2.0*electron_mass_c2/beta; |
---|
924 | G4double XHDT0 = std::log(gammaSq)-beta; |
---|
925 | |
---|
926 | G4double cpSq = energy*(energy+2.0*electron_mass_c2); |
---|
927 | G4double cp = std::sqrt(cpSq); |
---|
928 | G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2)); |
---|
929 | |
---|
930 | // |
---|
931 | // Distant interactions |
---|
932 | // |
---|
933 | G4double resEne = theOsc->GetResonanceEnergy(); |
---|
934 | G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy(); |
---|
935 | if (energy > resEne) |
---|
936 | { |
---|
937 | G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2); |
---|
938 | G4double cp1 = std::sqrt(cp1Sq); |
---|
939 | |
---|
940 | //Distant longitudinal interactions |
---|
941 | G4double QM = 0; |
---|
942 | if (resEne > 1e-6*energy) |
---|
943 | QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2; |
---|
944 | else |
---|
945 | { |
---|
946 | QM = resEne*resEne/(beta*2.0*electron_mass_c2); |
---|
947 | QM = QM*(1.0-0.5*QM/electron_mass_c2); |
---|
948 | } |
---|
949 | G4double SDL1 = 0; |
---|
950 | if (QM < cutoffEne) |
---|
951 | SDL1 = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2))); |
---|
952 | |
---|
953 | //Distant transverse interactions |
---|
954 | if (SDL1) |
---|
955 | { |
---|
956 | G4double SDT1 = std::max(XHDT0-delta,0.0); |
---|
957 | G4double SD1 = SDL1+SDT1; |
---|
958 | if (cut > resEne) |
---|
959 | { |
---|
960 | S1 = SD1; //XS1 |
---|
961 | S0 = SD1/resEne; //XS0 |
---|
962 | S2 = SD1*resEne; //XS2 |
---|
963 | } |
---|
964 | else |
---|
965 | { |
---|
966 | H1 = SD1; //XH1 |
---|
967 | H0 = SD1/resEne; //XH0 |
---|
968 | H2 = SD1*resEne; //XH2 |
---|
969 | } |
---|
970 | } |
---|
971 | } |
---|
972 | // |
---|
973 | // Close collisions (Moller's cross section) |
---|
974 | // |
---|
975 | G4double wl = std::max(cut,cutoffEne); |
---|
976 | G4double ee = energy + ionEnergy; |
---|
977 | G4double wu = 0.5*ee; |
---|
978 | if (wl < wu-(1e-5*eV)) |
---|
979 | { |
---|
980 | H0 += (1.0/(ee-wu)) - (1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) + |
---|
981 | (1.0-amol)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee + |
---|
982 | amol*(wu-wl)/(ee*ee); |
---|
983 | H1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) + |
---|
984 | (2.0-amol)*std::log((ee-wu)/(ee-wl)) + |
---|
985 | amol*(wu*wu-wl*wl)/(2.0*ee*ee); |
---|
986 | H2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) - |
---|
987 | (wl*(2.0*ee-wl)/(ee-wl)) + |
---|
988 | (3.0-amol)*ee*std::log((ee-wu)/(ee-wl)) + |
---|
989 | amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee); |
---|
990 | wu = wl; |
---|
991 | } |
---|
992 | wl = cutoffEne; |
---|
993 | |
---|
994 | if (wl > wu-(1e-5*eV)) |
---|
995 | { |
---|
996 | (*result)[0] = constant*H0; |
---|
997 | (*result)[1] = constant*H1; |
---|
998 | (*result)[2] = constant*H2; |
---|
999 | (*result)[3] = constant*S0; |
---|
1000 | (*result)[4] = constant*S1; |
---|
1001 | (*result)[5] = constant*S2; |
---|
1002 | return result; |
---|
1003 | } |
---|
1004 | |
---|
1005 | S0 += (1.0/(ee-wu))-(1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) + |
---|
1006 | (1.0-amol)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee + |
---|
1007 | amol*(wu-wl)/(ee*ee); |
---|
1008 | S1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) + |
---|
1009 | (2.0-amol)*std::log((ee-wu)/(ee-wl)) + |
---|
1010 | amol*(wu*wu-wl*wl)/(2.0*ee*ee); |
---|
1011 | S2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) - |
---|
1012 | (wl*(2.0*ee-wl)/(ee-wl)) + |
---|
1013 | (3.0-amol)*ee*std::log((ee-wu)/(ee-wl)) + |
---|
1014 | amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee); |
---|
1015 | |
---|
1016 | (*result)[0] = constant*H0; |
---|
1017 | (*result)[1] = constant*H1; |
---|
1018 | (*result)[2] = constant*H2; |
---|
1019 | (*result)[3] = constant*S0; |
---|
1020 | (*result)[4] = constant*S1; |
---|
1021 | (*result)[5] = constant*S2; |
---|
1022 | return result; |
---|
1023 | } |
---|
1024 | |
---|
1025 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
1026 | G4DataVector* G4Penelope08IonisationModel::ComputeShellCrossSectionsPositron(G4PenelopeOscillator* theOsc, |
---|
1027 | G4double energy, |
---|
1028 | G4double cut, |
---|
1029 | G4double delta) |
---|
1030 | { |
---|
1031 | // |
---|
1032 | //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for |
---|
1033 | //the given oscillator/cut and at the given energy. |
---|
1034 | //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2) |
---|
1035 | //Equivalent of subroutines PINaT1 of Penelope |
---|
1036 | // |
---|
1037 | // Results are _per target electron_ |
---|
1038 | // |
---|
1039 | G4DataVector* result = new G4DataVector(); |
---|
1040 | for (size_t i=0;i<6;i++) |
---|
1041 | result->push_back(0.); |
---|
1042 | G4double ionEnergy = theOsc->GetIonisationEnergy(); |
---|
1043 | |
---|
1044 | //return a set of zero's if the energy it too low to excite the current oscillator |
---|
1045 | if (energy < ionEnergy) |
---|
1046 | return result; |
---|
1047 | |
---|
1048 | G4double H0=0.,H1=0.,H2=0.; |
---|
1049 | G4double S0=0.,S1=0.,S2=0.; |
---|
1050 | |
---|
1051 | //Define useful constants to be used in the calculation |
---|
1052 | G4double gamma = 1.0+energy/electron_mass_c2; |
---|
1053 | G4double gammaSq = gamma*gamma; |
---|
1054 | G4double beta = (gammaSq-1.0)/gammaSq; |
---|
1055 | G4double pielr2 = pi*classic_electr_radius*classic_electr_radius; //pi*re^2 |
---|
1056 | G4double constant = pielr2*2.0*electron_mass_c2/beta; |
---|
1057 | G4double XHDT0 = std::log(gammaSq)-beta; |
---|
1058 | |
---|
1059 | G4double cpSq = energy*(energy+2.0*electron_mass_c2); |
---|
1060 | G4double cp = std::sqrt(cpSq); |
---|
1061 | G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2)); |
---|
1062 | G4double g12 = (gamma+1.0)*(gamma+1.0); |
---|
1063 | //Bhabha coefficients |
---|
1064 | G4double bha1 = amol*(2.0*g12-1.0)/(gammaSq-1.0); |
---|
1065 | G4double bha2 = amol*(3.0+1.0/g12); |
---|
1066 | G4double bha3 = amol*2.0*gamma*(gamma-1.0)/g12; |
---|
1067 | G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/g12; |
---|
1068 | |
---|
1069 | // |
---|
1070 | // Distant interactions |
---|
1071 | // |
---|
1072 | G4double resEne = theOsc->GetResonanceEnergy(); |
---|
1073 | G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy(); |
---|
1074 | if (energy > resEne) |
---|
1075 | { |
---|
1076 | G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2); |
---|
1077 | G4double cp1 = std::sqrt(cp1Sq); |
---|
1078 | |
---|
1079 | //Distant longitudinal interactions |
---|
1080 | G4double QM = 0; |
---|
1081 | if (resEne > 1e-6*energy) |
---|
1082 | QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2; |
---|
1083 | else |
---|
1084 | { |
---|
1085 | QM = resEne*resEne/(beta*2.0*electron_mass_c2); |
---|
1086 | QM = QM*(1.0-0.5*QM/electron_mass_c2); |
---|
1087 | } |
---|
1088 | G4double SDL1 = 0; |
---|
1089 | if (QM < cutoffEne) |
---|
1090 | SDL1 = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2))); |
---|
1091 | |
---|
1092 | //Distant transverse interactions |
---|
1093 | if (SDL1) |
---|
1094 | { |
---|
1095 | G4double SDT1 = std::max(XHDT0-delta,0.0); |
---|
1096 | G4double SD1 = SDL1+SDT1; |
---|
1097 | if (cut > resEne) |
---|
1098 | { |
---|
1099 | S1 = SD1; //XS1 |
---|
1100 | S0 = SD1/resEne; //XS0 |
---|
1101 | S2 = SD1*resEne; //XS2 |
---|
1102 | } |
---|
1103 | else |
---|
1104 | { |
---|
1105 | H1 = SD1; //XH1 |
---|
1106 | H0 = SD1/resEne; //XH0 |
---|
1107 | H2 = SD1*resEne; //XH2 |
---|
1108 | } |
---|
1109 | } |
---|
1110 | } |
---|
1111 | |
---|
1112 | // |
---|
1113 | // Close collisions (Bhabha's cross section) |
---|
1114 | // |
---|
1115 | G4double wl = std::max(cut,cutoffEne); |
---|
1116 | G4double wu = energy; |
---|
1117 | G4double energySq = energy*energy; |
---|
1118 | if (wl < wu-(1e-5*eV)) |
---|
1119 | { |
---|
1120 | G4double wlSq = wl*wl; |
---|
1121 | G4double wuSq = wu*wu; |
---|
1122 | H0 += (1.0/wl) - (1.0/wu)- bha1*std::log(wu/wl)/energy |
---|
1123 | + bha2*(wu-wl)/energySq |
---|
1124 | - bha3*(wuSq-wlSq)/(2.0*energySq*energy) |
---|
1125 | + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq); |
---|
1126 | H1 += std::log(wu/wl) - bha1*(wu-wl)/energy |
---|
1127 | + bha2*(wuSq-wlSq)/(2.0*energySq) |
---|
1128 | - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energy) |
---|
1129 | + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq); |
---|
1130 | H2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy) |
---|
1131 | + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq) |
---|
1132 | - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy) |
---|
1133 | + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq); |
---|
1134 | wu = wl; |
---|
1135 | } |
---|
1136 | wl = cutoffEne; |
---|
1137 | |
---|
1138 | if (wl > wu-(1e-5*eV)) |
---|
1139 | { |
---|
1140 | (*result)[0] = constant*H0; |
---|
1141 | (*result)[1] = constant*H1; |
---|
1142 | (*result)[2] = constant*H2; |
---|
1143 | (*result)[3] = constant*S0; |
---|
1144 | (*result)[4] = constant*S1; |
---|
1145 | (*result)[5] = constant*S2; |
---|
1146 | return result; |
---|
1147 | } |
---|
1148 | |
---|
1149 | G4double wlSq = wl*wl; |
---|
1150 | G4double wuSq = wu*wu; |
---|
1151 | |
---|
1152 | S0 += (1.0/wl) - (1.0/wu) - bha1*std::log(wu/wl)/energy |
---|
1153 | + bha2*(wu-wl)/energySq |
---|
1154 | - bha3*(wuSq-wlSq)/(2.0*energySq*energy) |
---|
1155 | + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq); |
---|
1156 | |
---|
1157 | S1 += std::log(wu/wl) - bha1*(wu-wl)/energy |
---|
1158 | + bha2*(wuSq-wlSq)/(2.0*energySq) |
---|
1159 | - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energy) |
---|
1160 | + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq); |
---|
1161 | |
---|
1162 | S2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy) |
---|
1163 | + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq) |
---|
1164 | - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy) |
---|
1165 | + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq); |
---|
1166 | |
---|
1167 | (*result)[0] = constant*H0; |
---|
1168 | (*result)[1] = constant*H1; |
---|
1169 | (*result)[2] = constant*H2; |
---|
1170 | (*result)[3] = constant*S0; |
---|
1171 | (*result)[4] = constant*S1; |
---|
1172 | (*result)[5] = constant*S2; |
---|
1173 | |
---|
1174 | return result; |
---|
1175 | } |
---|
1176 | |
---|
1177 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
1178 | void G4Penelope08IonisationModel::SampleFinalStateElectron(const G4Material* mat, |
---|
1179 | G4double cutEnergy, |
---|
1180 | G4double kineticEnergy) |
---|
1181 | { |
---|
1182 | // This method sets the final ionisation parameters |
---|
1183 | // kineticEnergy1, cosThetaPrimary (= updates of the primary e-) |
---|
1184 | // energySecondary, cosThetaSecondary (= info of delta-ray) |
---|
1185 | // targetOscillator (= ionised oscillator) |
---|
1186 | // |
---|
1187 | // The method implements SUBROUTINE EINa of Penelope |
---|
1188 | // |
---|
1189 | |
---|
1190 | G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat); |
---|
1191 | size_t numberOfOscillators = theTable->size(); |
---|
1192 | G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(G4Electron::Electron(),mat, |
---|
1193 | cutEnergy); |
---|
1194 | G4double delta = GetDensityCorrection(mat,kineticEnergy); |
---|
1195 | |
---|
1196 | // Selection of the active oscillator |
---|
1197 | G4double TST = G4UniformRand(); |
---|
1198 | targetOscillator = numberOfOscillators-1; //initialization, last oscillator |
---|
1199 | G4double XSsum = 0.; |
---|
1200 | |
---|
1201 | for (size_t i=0;i<numberOfOscillators-1;i++) |
---|
1202 | { |
---|
1203 | XSsum += theXS->GetShellCrossSection(i,kineticEnergy); |
---|
1204 | |
---|
1205 | if (XSsum > TST) |
---|
1206 | { |
---|
1207 | targetOscillator = (G4int) i; |
---|
1208 | break; |
---|
1209 | } |
---|
1210 | } |
---|
1211 | |
---|
1212 | |
---|
1213 | if (verboseLevel > 3) |
---|
1214 | { |
---|
1215 | G4cout << "SampleFinalStateElectron: sampled oscillator #" << targetOscillator << "." << G4endl; |
---|
1216 | G4cout << "Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/eV << |
---|
1217 | " eV " << G4endl; |
---|
1218 | G4cout << "Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/eV << " eV " |
---|
1219 | << G4endl; |
---|
1220 | } |
---|
1221 | //Constants |
---|
1222 | G4double rb = kineticEnergy + 2.0*electron_mass_c2; |
---|
1223 | G4double gam = 1.0+kineticEnergy/electron_mass_c2; |
---|
1224 | G4double gam2 = gam*gam; |
---|
1225 | G4double beta2 = (gam2-1.0)/gam2; |
---|
1226 | G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam); |
---|
1227 | |
---|
1228 | //Partial cross section of the active oscillator |
---|
1229 | G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy(); |
---|
1230 | G4double invResEne = 1.0/resEne; |
---|
1231 | G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy(); |
---|
1232 | G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy(); |
---|
1233 | G4double XHDL = 0.; |
---|
1234 | G4double XHDT = 0.; |
---|
1235 | G4double QM = 0.; |
---|
1236 | G4double cps = 0.; |
---|
1237 | G4double cp = 0.; |
---|
1238 | |
---|
1239 | //Distant excitations |
---|
1240 | if (resEne > cutEnergy && resEne < kineticEnergy) |
---|
1241 | { |
---|
1242 | cps = kineticEnergy*rb; |
---|
1243 | cp = std::sqrt(cps); |
---|
1244 | G4double XHDT0 = std::max(std::log(gam2)-beta2-delta,0.); |
---|
1245 | if (resEne > 1.0e-6*kineticEnergy) |
---|
1246 | { |
---|
1247 | G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2)); |
---|
1248 | QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2; |
---|
1249 | } |
---|
1250 | else |
---|
1251 | { |
---|
1252 | QM = resEne*resEne/(beta2*2.0*electron_mass_c2); |
---|
1253 | QM *= (1.0-QM*0.5/electron_mass_c2); |
---|
1254 | } |
---|
1255 | if (QM < cutoffEne) |
---|
1256 | { |
---|
1257 | XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2))) |
---|
1258 | *invResEne; |
---|
1259 | XHDT = XHDT0*invResEne; |
---|
1260 | } |
---|
1261 | else |
---|
1262 | { |
---|
1263 | QM = cutoffEne; |
---|
1264 | XHDL = 0.; |
---|
1265 | XHDT = 0.; |
---|
1266 | } |
---|
1267 | } |
---|
1268 | else |
---|
1269 | { |
---|
1270 | QM = cutoffEne; |
---|
1271 | cps = 0.; |
---|
1272 | cp = 0.; |
---|
1273 | XHDL = 0.; |
---|
1274 | XHDT = 0.; |
---|
1275 | } |
---|
1276 | |
---|
1277 | //Close collisions |
---|
1278 | G4double EE = kineticEnergy + ionEne; |
---|
1279 | G4double wmaxc = 0.5*EE; |
---|
1280 | G4double wcl = std::max(cutEnergy,cutoffEne); |
---|
1281 | G4double rcl = wcl/EE; |
---|
1282 | G4double XHC = 0.; |
---|
1283 | if (wcl < wmaxc) |
---|
1284 | { |
---|
1285 | G4double rl1 = 1.0-rcl; |
---|
1286 | G4double rrl1 = 1.0/rl1; |
---|
1287 | XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+ |
---|
1288 | (1.0-amol)*std::log(rcl*rrl1))/EE; |
---|
1289 | } |
---|
1290 | |
---|
1291 | //Total cross section per molecule for the active shell, in cm2 |
---|
1292 | G4double XHTOT = XHC + XHDL + XHDT; |
---|
1293 | |
---|
1294 | //very small cross section, do nothing |
---|
1295 | if (XHTOT < 1.e-14*barn) |
---|
1296 | { |
---|
1297 | kineticEnergy1=kineticEnergy; |
---|
1298 | cosThetaPrimary=1.0; |
---|
1299 | energySecondary=0.0; |
---|
1300 | cosThetaSecondary=1.0; |
---|
1301 | targetOscillator = numberOfOscillators-1; |
---|
1302 | return; |
---|
1303 | } |
---|
1304 | |
---|
1305 | //decide which kind of interaction we'll have |
---|
1306 | TST = XHTOT*G4UniformRand(); |
---|
1307 | |
---|
1308 | |
---|
1309 | // Hard close collision |
---|
1310 | G4double TS1 = XHC; |
---|
1311 | |
---|
1312 | if (TST < TS1) |
---|
1313 | { |
---|
1314 | G4double A = 5.0*amol; |
---|
1315 | G4double ARCL = A*0.5*rcl; |
---|
1316 | G4double rk=0.; |
---|
1317 | G4bool loopAgain = false; |
---|
1318 | do |
---|
1319 | { |
---|
1320 | loopAgain = false; |
---|
1321 | G4double fb = (1.0+ARCL)*G4UniformRand(); |
---|
1322 | if (fb < 1) |
---|
1323 | rk = rcl/(1.0-fb*(1.0-(rcl+rcl))); |
---|
1324 | else |
---|
1325 | rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL; |
---|
1326 | G4double rk2 = rk*rk; |
---|
1327 | G4double rkf = rk/(1.0-rk); |
---|
1328 | G4double phi = 1.0+rkf*rkf-rkf+amol*(rk2+rkf); |
---|
1329 | if (G4UniformRand()*(1.0+A*rk2) > phi) |
---|
1330 | loopAgain = true; |
---|
1331 | }while(loopAgain); |
---|
1332 | //energy and scattering angle (primary electron) |
---|
1333 | G4double deltaE = rk*EE; |
---|
1334 | |
---|
1335 | kineticEnergy1 = kineticEnergy - deltaE; |
---|
1336 | cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE))); |
---|
1337 | //energy and scattering angle of the delta ray |
---|
1338 | energySecondary = deltaE - ionEne; //subtract ionisation energy |
---|
1339 | cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2))); |
---|
1340 | if (verboseLevel > 3) |
---|
1341 | G4cout << "SampleFinalStateElectron: sampled close collision " << G4endl; |
---|
1342 | return; |
---|
1343 | } |
---|
1344 | |
---|
1345 | //Hard distant longitudinal collisions |
---|
1346 | TS1 += XHDL; |
---|
1347 | G4double deltaE = resEne; |
---|
1348 | kineticEnergy1 = kineticEnergy - deltaE; |
---|
1349 | |
---|
1350 | if (TST < TS1) |
---|
1351 | { |
---|
1352 | G4double QS = QM/(1.0+QM*0.5/electron_mass_c2); |
---|
1353 | G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand()) |
---|
1354 | - (QS*0.5/electron_mass_c2)); |
---|
1355 | G4double QTREV = Q*(Q+2.0*electron_mass_c2); |
---|
1356 | G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2); |
---|
1357 | cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps)); |
---|
1358 | if (cosThetaPrimary > 1.) |
---|
1359 | cosThetaPrimary = 1.0; |
---|
1360 | //energy and emission angle of the delta ray |
---|
1361 | energySecondary = deltaE - ionEne; |
---|
1362 | cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV); |
---|
1363 | if (cosThetaSecondary > 1.0) |
---|
1364 | cosThetaSecondary = 1.0; |
---|
1365 | if (verboseLevel > 3) |
---|
1366 | G4cout << "SampleFinalStateElectron: sampled distant longitudinal collision " << G4endl; |
---|
1367 | return; |
---|
1368 | } |
---|
1369 | |
---|
1370 | //Hard distant transverse collisions |
---|
1371 | cosThetaPrimary = 1.0; |
---|
1372 | //energy and emission angle of the delta ray |
---|
1373 | energySecondary = deltaE - ionEne; |
---|
1374 | cosThetaSecondary = 0.5; |
---|
1375 | if (verboseLevel > 3) |
---|
1376 | G4cout << "SampleFinalStateElectron: sampled distant transverse collision " << G4endl; |
---|
1377 | |
---|
1378 | return; |
---|
1379 | } |
---|
1380 | |
---|
1381 | |
---|
1382 | |
---|
1383 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
1384 | void G4Penelope08IonisationModel::SampleFinalStatePositron(const G4Material* mat, |
---|
1385 | G4double cutEnergy, |
---|
1386 | G4double kineticEnergy) |
---|
1387 | { |
---|
1388 | // This method sets the final ionisation parameters |
---|
1389 | // kineticEnergy1, cosThetaPrimary (= updates of the primary e-) |
---|
1390 | // energySecondary, cosThetaSecondary (= info of delta-ray) |
---|
1391 | // targetOscillator (= ionised oscillator) |
---|
1392 | // |
---|
1393 | // The method implements SUBROUTINE PINa of Penelope |
---|
1394 | // |
---|
1395 | |
---|
1396 | G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat); |
---|
1397 | size_t numberOfOscillators = theTable->size(); |
---|
1398 | G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(G4Positron::Positron(),mat, |
---|
1399 | cutEnergy); |
---|
1400 | G4double delta = GetDensityCorrection(mat,kineticEnergy); |
---|
1401 | |
---|
1402 | // Selection of the active oscillator |
---|
1403 | G4double TST = G4UniformRand(); |
---|
1404 | targetOscillator = numberOfOscillators-1; //initialization, last oscillator |
---|
1405 | G4double XSsum = 0.; |
---|
1406 | for (size_t i=0;i<numberOfOscillators-1;i++) |
---|
1407 | { |
---|
1408 | XSsum += theXS->GetShellCrossSection(i,kineticEnergy); |
---|
1409 | if (XSsum > TST) |
---|
1410 | { |
---|
1411 | targetOscillator = (G4int) i; |
---|
1412 | break; |
---|
1413 | } |
---|
1414 | } |
---|
1415 | |
---|
1416 | if (verboseLevel > 3) |
---|
1417 | { |
---|
1418 | G4cout << "SampleFinalStatePositron: sampled oscillator #" << targetOscillator << "." << G4endl; |
---|
1419 | G4cout << "Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/eV |
---|
1420 | << " eV " << G4endl; |
---|
1421 | G4cout << "Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/eV |
---|
1422 | << " eV " << G4endl; |
---|
1423 | } |
---|
1424 | |
---|
1425 | //Constants |
---|
1426 | G4double rb = kineticEnergy + 2.0*electron_mass_c2; |
---|
1427 | G4double gam = 1.0+kineticEnergy/electron_mass_c2; |
---|
1428 | G4double gam2 = gam*gam; |
---|
1429 | G4double beta2 = (gam2-1.0)/gam2; |
---|
1430 | G4double g12 = (gam+1.0)*(gam+1.0); |
---|
1431 | G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam); |
---|
1432 | //Bhabha coefficients |
---|
1433 | G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0); |
---|
1434 | G4double bha2 = amol*(3.0+1.0/g12); |
---|
1435 | G4double bha3 = amol*2.0*gam*(gam-1.0)/g12; |
---|
1436 | G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12; |
---|
1437 | |
---|
1438 | // |
---|
1439 | //Partial cross section of the active oscillator |
---|
1440 | // |
---|
1441 | G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy(); |
---|
1442 | G4double invResEne = 1.0/resEne; |
---|
1443 | G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy(); |
---|
1444 | G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy(); |
---|
1445 | |
---|
1446 | G4double XHDL = 0.; |
---|
1447 | G4double XHDT = 0.; |
---|
1448 | G4double QM = 0.; |
---|
1449 | G4double cps = 0.; |
---|
1450 | G4double cp = 0.; |
---|
1451 | |
---|
1452 | //Distant excitations XS (same as for electrons) |
---|
1453 | if (resEne > cutEnergy && resEne < kineticEnergy) |
---|
1454 | { |
---|
1455 | cps = kineticEnergy*rb; |
---|
1456 | cp = std::sqrt(cps); |
---|
1457 | G4double XHDT0 = std::max(std::log(gam2)-beta2-delta,0.); |
---|
1458 | if (resEne > 1.0e-6*kineticEnergy) |
---|
1459 | { |
---|
1460 | G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2)); |
---|
1461 | QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2; |
---|
1462 | } |
---|
1463 | else |
---|
1464 | { |
---|
1465 | QM = resEne*resEne/(beta2*2.0*electron_mass_c2); |
---|
1466 | QM *= (1.0-QM*0.5/electron_mass_c2); |
---|
1467 | } |
---|
1468 | if (QM < cutoffEne) |
---|
1469 | { |
---|
1470 | XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2))) |
---|
1471 | *invResEne; |
---|
1472 | XHDT = XHDT0*invResEne; |
---|
1473 | } |
---|
1474 | else |
---|
1475 | { |
---|
1476 | QM = cutoffEne; |
---|
1477 | XHDL = 0.; |
---|
1478 | XHDT = 0.; |
---|
1479 | } |
---|
1480 | } |
---|
1481 | else |
---|
1482 | { |
---|
1483 | QM = cutoffEne; |
---|
1484 | cps = 0.; |
---|
1485 | cp = 0.; |
---|
1486 | XHDL = 0.; |
---|
1487 | XHDT = 0.; |
---|
1488 | } |
---|
1489 | //Close collisions (Bhabha) |
---|
1490 | G4double wmaxc = kineticEnergy; |
---|
1491 | G4double wcl = std::max(cutEnergy,cutoffEne); |
---|
1492 | G4double rcl = wcl/kineticEnergy; |
---|
1493 | G4double XHC = 0.; |
---|
1494 | if (wcl < wmaxc) |
---|
1495 | { |
---|
1496 | G4double rl1 = 1.0-rcl; |
---|
1497 | XHC = ((1.0/rcl-1.0)+bha1*std::log(rcl)+bha2*rl1 |
---|
1498 | + (bha3/2.0)*(rcl*rcl-1.0) |
---|
1499 | + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineticEnergy; |
---|
1500 | } |
---|
1501 | |
---|
1502 | //Total cross section per molecule for the active shell, in cm2 |
---|
1503 | G4double XHTOT = XHC + XHDL + XHDT; |
---|
1504 | |
---|
1505 | //very small cross section, do nothing |
---|
1506 | if (XHTOT < 1.e-14*barn) |
---|
1507 | { |
---|
1508 | kineticEnergy1=kineticEnergy; |
---|
1509 | cosThetaPrimary=1.0; |
---|
1510 | energySecondary=0.0; |
---|
1511 | cosThetaSecondary=1.0; |
---|
1512 | targetOscillator = numberOfOscillators-1; |
---|
1513 | return; |
---|
1514 | } |
---|
1515 | |
---|
1516 | //decide which kind of interaction we'll have |
---|
1517 | TST = XHTOT*G4UniformRand(); |
---|
1518 | |
---|
1519 | // Hard close collision |
---|
1520 | G4double TS1 = XHC; |
---|
1521 | if (TST < TS1) |
---|
1522 | { |
---|
1523 | G4double rl1 = 1.0-rcl; |
---|
1524 | G4double rk=0.; |
---|
1525 | G4bool loopAgain = false; |
---|
1526 | do |
---|
1527 | { |
---|
1528 | loopAgain = false; |
---|
1529 | rk = rcl/(1.0-G4UniformRand()*rl1); |
---|
1530 | G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk))); |
---|
1531 | if (G4UniformRand() > phi) |
---|
1532 | loopAgain = true; |
---|
1533 | }while(loopAgain); |
---|
1534 | //energy and scattering angle (primary electron) |
---|
1535 | G4double deltaE = rk*kineticEnergy; |
---|
1536 | kineticEnergy1 = kineticEnergy - deltaE; |
---|
1537 | cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE))); |
---|
1538 | //energy and scattering angle of the delta ray |
---|
1539 | energySecondary = deltaE - ionEne; //subtract ionisation energy |
---|
1540 | cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2))); |
---|
1541 | if (verboseLevel > 3) |
---|
1542 | G4cout << "SampleFinalStatePositron: sampled close collision " << G4endl; |
---|
1543 | return; |
---|
1544 | } |
---|
1545 | |
---|
1546 | //Hard distant longitudinal collisions |
---|
1547 | TS1 += XHDL; |
---|
1548 | G4double deltaE = resEne; |
---|
1549 | kineticEnergy1 = kineticEnergy - deltaE; |
---|
1550 | if (TST < TS1) |
---|
1551 | { |
---|
1552 | G4double QS = QM/(1.0+QM*0.5/electron_mass_c2); |
---|
1553 | G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand()) |
---|
1554 | - (QS*0.5/electron_mass_c2)); |
---|
1555 | G4double QTREV = Q*(Q+2.0*electron_mass_c2); |
---|
1556 | G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2); |
---|
1557 | cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps)); |
---|
1558 | if (cosThetaPrimary > 1.) |
---|
1559 | cosThetaPrimary = 1.0; |
---|
1560 | //energy and emission angle of the delta ray |
---|
1561 | energySecondary = deltaE - ionEne; |
---|
1562 | cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV); |
---|
1563 | if (cosThetaSecondary > 1.0) |
---|
1564 | cosThetaSecondary = 1.0; |
---|
1565 | if (verboseLevel > 3) |
---|
1566 | G4cout << "SampleFinalStatePositron: sampled distant longitudinal collision " << G4endl; |
---|
1567 | return; |
---|
1568 | } |
---|
1569 | |
---|
1570 | //Hard distant transverse collisions |
---|
1571 | cosThetaPrimary = 1.0; |
---|
1572 | //energy and emission angle of the delta ray |
---|
1573 | energySecondary = deltaE - ionEne; |
---|
1574 | cosThetaSecondary = 0.5; |
---|
1575 | |
---|
1576 | if (verboseLevel > 3) |
---|
1577 | G4cout << "SampleFinalStatePositron: sampled distant transverse collision " << G4endl; |
---|
1578 | |
---|
1579 | return; |
---|
1580 | } |
---|
1581 | |
---|