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: G4Penelope08ComptonModel.cc,v 1.7 2010/07/28 07:09:16 pandola Exp $ |
---|
27 | // GEANT4 tag $Name: geant4-09-04-ref-00 $ |
---|
28 | // |
---|
29 | // Author: Luciano Pandola |
---|
30 | // |
---|
31 | // History: |
---|
32 | // -------- |
---|
33 | // 15 Feb 2010 L Pandola Implementation |
---|
34 | // 18 Mar 2010 L. Pandola Removed GetAtomsPerMolecule(), now demanded |
---|
35 | // to G4PenelopeOscillatorManager |
---|
36 | // |
---|
37 | #include "G4Penelope08ComptonModel.hh" |
---|
38 | #include "G4ParticleDefinition.hh" |
---|
39 | #include "G4MaterialCutsCouple.hh" |
---|
40 | #include "G4ProductionCutsTable.hh" |
---|
41 | #include "G4DynamicParticle.hh" |
---|
42 | #include "G4VEMDataSet.hh" |
---|
43 | #include "G4PhysicsTable.hh" |
---|
44 | #include "G4PhysicsLogVector.hh" |
---|
45 | #include "G4AtomicTransitionManager.hh" |
---|
46 | #include "G4AtomicShell.hh" |
---|
47 | #include "G4Gamma.hh" |
---|
48 | #include "G4Electron.hh" |
---|
49 | #include "G4PenelopeOscillatorManager.hh" |
---|
50 | #include "G4PenelopeOscillator.hh" |
---|
51 | |
---|
52 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
53 | |
---|
54 | |
---|
55 | G4Penelope08ComptonModel::G4Penelope08ComptonModel(const G4ParticleDefinition*, |
---|
56 | const G4String& nam) |
---|
57 | :G4VEmModel(nam),isInitialised(false),oscManager(0) |
---|
58 | { |
---|
59 | fIntrinsicLowEnergyLimit = 100.0*eV; |
---|
60 | fIntrinsicHighEnergyLimit = 100.0*GeV; |
---|
61 | // SetLowEnergyLimit(fIntrinsicLowEnergyLimit); |
---|
62 | SetHighEnergyLimit(fIntrinsicHighEnergyLimit); |
---|
63 | // |
---|
64 | oscManager = G4PenelopeOscillatorManager::GetOscillatorManager(); |
---|
65 | |
---|
66 | verboseLevel= 0; |
---|
67 | // Verbosity scale: |
---|
68 | // 0 = nothing |
---|
69 | // 1 = warning for energy non-conservation |
---|
70 | // 2 = details of energy budget |
---|
71 | // 3 = calculation of cross sections, file openings, sampling of atoms |
---|
72 | // 4 = entering in methods |
---|
73 | |
---|
74 | //by default, the model will use atomic deexcitation |
---|
75 | SetDeexcitationFlag(true); |
---|
76 | ActivateAuger(false); |
---|
77 | |
---|
78 | } |
---|
79 | |
---|
80 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
81 | |
---|
82 | G4Penelope08ComptonModel::~G4Penelope08ComptonModel() |
---|
83 | {;} |
---|
84 | |
---|
85 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
86 | |
---|
87 | void G4Penelope08ComptonModel::Initialise(const G4ParticleDefinition*, |
---|
88 | const G4DataVector&) |
---|
89 | { |
---|
90 | if (verboseLevel > 3) |
---|
91 | G4cout << "Calling G4Penelope08ComptonModel::Initialise()" << G4endl; |
---|
92 | |
---|
93 | if (verboseLevel > 0) { |
---|
94 | G4cout << "Penelope Compton model is initialized " << G4endl |
---|
95 | << "Energy range: " |
---|
96 | << LowEnergyLimit() / keV << " keV - " |
---|
97 | << HighEnergyLimit() / GeV << " GeV" |
---|
98 | << G4endl; |
---|
99 | } |
---|
100 | |
---|
101 | if(isInitialised) return; |
---|
102 | fParticleChange = GetParticleChangeForGamma(); |
---|
103 | isInitialised = true; |
---|
104 | } |
---|
105 | |
---|
106 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
107 | |
---|
108 | G4double G4Penelope08ComptonModel::CrossSectionPerVolume(const G4Material* material, |
---|
109 | const G4ParticleDefinition* p, |
---|
110 | G4double energy, |
---|
111 | G4double, |
---|
112 | G4double) |
---|
113 | { |
---|
114 | // Penelope model to calculate the Compton scattering cross section: |
---|
115 | // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167 |
---|
116 | // |
---|
117 | // The cross section for Compton scattering is calculated according to the Klein-Nishina |
---|
118 | // formula for energy > 5 MeV. |
---|
119 | // For E < 5 MeV it is used a parametrization for the differential cross-section dSigma/dOmega, |
---|
120 | // which is integrated numerically in cos(theta), G4Penelope08ComptonModel::DifferentialCrossSection(). |
---|
121 | // The parametrization includes the J(p) |
---|
122 | // distribution profiles for the atomic shells, that are tabulated from Hartree-Fock calculations |
---|
123 | // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201 |
---|
124 | // |
---|
125 | if (verboseLevel > 3) |
---|
126 | G4cout << "Calling CrossSectionPerVolume() of G4Penelope08ComptonModel" << G4endl; |
---|
127 | SetupForMaterial(p, material, energy); |
---|
128 | |
---|
129 | //Retrieve the oscillator table for this material |
---|
130 | G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material); |
---|
131 | |
---|
132 | G4double cs = 0; |
---|
133 | |
---|
134 | if (energy < 5*MeV) //explicit calculation for E < 5 MeV |
---|
135 | { |
---|
136 | size_t numberOfOscillators = theTable->size(); |
---|
137 | for (size_t i=0;i<numberOfOscillators;i++) |
---|
138 | { |
---|
139 | G4PenelopeOscillator* theOsc = (*theTable)[i]; |
---|
140 | //sum contributions coming from each oscillator |
---|
141 | cs += OscillatorTotalCrossSection(energy,theOsc); |
---|
142 | } |
---|
143 | } |
---|
144 | else //use Klein-Nishina for E>5 MeV |
---|
145 | cs = KleinNishinaCrossSection(energy,material); |
---|
146 | |
---|
147 | //cross sections are in units of pi*classic_electr_radius^2 |
---|
148 | cs *= pi*classic_electr_radius*classic_electr_radius; |
---|
149 | |
---|
150 | //Now, cs is the cross section *per molecule*, let's calculate the |
---|
151 | //cross section per volume |
---|
152 | |
---|
153 | G4double atomDensity = material->GetTotNbOfAtomsPerVolume(); |
---|
154 | G4double atPerMol = oscManager->GetAtomsPerMolecule(material); |
---|
155 | |
---|
156 | if (verboseLevel > 3) |
---|
157 | G4cout << "Material " << material->GetName() << " has " << atPerMol << |
---|
158 | "atoms per molecule" << G4endl; |
---|
159 | |
---|
160 | G4double moleculeDensity = 0.; |
---|
161 | |
---|
162 | if (atPerMol) |
---|
163 | moleculeDensity = atomDensity/atPerMol; |
---|
164 | |
---|
165 | G4double csvolume = cs*moleculeDensity; |
---|
166 | |
---|
167 | if (verboseLevel > 2) |
---|
168 | G4cout << "Compton mean free path at " << energy/keV << " keV for material " << |
---|
169 | material->GetName() << " = " << (1./csvolume)/mm << " mm" << G4endl; |
---|
170 | return csvolume; |
---|
171 | } |
---|
172 | |
---|
173 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
174 | |
---|
175 | //This is a dummy method. Never inkoved by the tracking, it just issues |
---|
176 | //a warning if one tries to get Cross Sections per Atom via the |
---|
177 | //G4EmCalculator. |
---|
178 | G4double G4Penelope08ComptonModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, |
---|
179 | G4double, |
---|
180 | G4double, |
---|
181 | G4double, |
---|
182 | G4double, |
---|
183 | G4double) |
---|
184 | { |
---|
185 | G4cout << "*** G4Penelope08ComptonModel -- WARNING ***" << G4endl; |
---|
186 | G4cout << "Penelope Compton model does not calculate cross section _per atom_ " << G4endl; |
---|
187 | G4cout << "so the result is always zero. For physics values, please invoke " << G4endl; |
---|
188 | G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl; |
---|
189 | return 0; |
---|
190 | } |
---|
191 | |
---|
192 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
193 | |
---|
194 | void G4Penelope08ComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, |
---|
195 | const G4MaterialCutsCouple* couple, |
---|
196 | const G4DynamicParticle* aDynamicGamma, |
---|
197 | G4double, |
---|
198 | G4double) |
---|
199 | { |
---|
200 | |
---|
201 | // Penelope model to sample the Compton scattering final state. |
---|
202 | // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167 |
---|
203 | // The model determines also the original shell from which the electron is expelled, |
---|
204 | // in order to produce fluorescence de-excitation (from G4DeexcitationManager) |
---|
205 | // |
---|
206 | // The final state for Compton scattering is calculated according to the Klein-Nishina |
---|
207 | // formula for energy > 5 MeV. In this case, the Doppler broadening is negligible and |
---|
208 | // one can assume that the target electron is at rest. |
---|
209 | // For E < 5 MeV it is used the parametrization for the differential cross-section dSigma/dOmega, |
---|
210 | // to sample the scattering angle and the energy of the emerging electron, which is |
---|
211 | // G4Penelope08ComptonModel::DifferentialCrossSection(). The rejection method is |
---|
212 | // used to sample cos(theta). The efficiency increases monotonically with photon energy and is |
---|
213 | // nearly independent on the Z; typical values are 35%, 80% and 95% for 1 keV, 1 MeV and 10 MeV, |
---|
214 | // respectively. |
---|
215 | // The parametrization includes the J(p) distribution profiles for the atomic shells, that are |
---|
216 | // tabulated |
---|
217 | // from Hartree-Fock calculations from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201. |
---|
218 | // Doppler broadening is included. |
---|
219 | // |
---|
220 | |
---|
221 | if (verboseLevel > 3) |
---|
222 | G4cout << "Calling SampleSecondaries() of G4Penelope08ComptonModel" << G4endl; |
---|
223 | |
---|
224 | G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy(); |
---|
225 | |
---|
226 | if (photonEnergy0 <= fIntrinsicLowEnergyLimit) |
---|
227 | { |
---|
228 | fParticleChange->ProposeTrackStatus(fStopAndKill); |
---|
229 | fParticleChange->SetProposedKineticEnergy(0.); |
---|
230 | fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0); |
---|
231 | return ; |
---|
232 | } |
---|
233 | |
---|
234 | G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection(); |
---|
235 | const G4Material* material = couple->GetMaterial(); |
---|
236 | |
---|
237 | G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material); |
---|
238 | |
---|
239 | const G4int nmax = 64; |
---|
240 | G4double rn[nmax],pac[nmax]; |
---|
241 | |
---|
242 | G4double S=0.0; |
---|
243 | G4double epsilon = 0.0; |
---|
244 | G4double cosTheta = 1.0; |
---|
245 | G4double hartreeFunc = 0.0; |
---|
246 | G4double oscStren = 0.0; |
---|
247 | size_t numberOfOscillators = theTable->size(); |
---|
248 | size_t targetOscillator = 0; |
---|
249 | G4double ionEnergy = 0.0*eV; |
---|
250 | |
---|
251 | G4double ek = photonEnergy0/electron_mass_c2; |
---|
252 | G4double ek2 = 2.*ek+1.0; |
---|
253 | G4double eks = ek*ek; |
---|
254 | G4double ek1 = eks-ek2-1.0; |
---|
255 | |
---|
256 | G4double taumin = 1.0/ek2; |
---|
257 | G4double a1 = std::log(ek2); |
---|
258 | G4double a2 = a1+2.0*ek*(1.0+ek)/(ek2*ek2); |
---|
259 | |
---|
260 | G4double TST = 0; |
---|
261 | G4double tau = 0.; |
---|
262 | |
---|
263 | //If the incoming photon is above 5 MeV, the quicker approach based on the |
---|
264 | //pure Klein-Nishina formula is used |
---|
265 | if (photonEnergy0 > 5*MeV) |
---|
266 | { |
---|
267 | do{ |
---|
268 | do{ |
---|
269 | if ((a2*G4UniformRand()) < a1) |
---|
270 | tau = std::pow(taumin,G4UniformRand()); |
---|
271 | else |
---|
272 | tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0)); |
---|
273 | //rejection function |
---|
274 | TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau)); |
---|
275 | }while (G4UniformRand()> TST); |
---|
276 | epsilon=tau; |
---|
277 | cosTheta = 1.0 - (1.0-tau)/(ek*tau); |
---|
278 | |
---|
279 | //Target shell electrons |
---|
280 | TST = oscManager->GetTotalZ(material)*G4UniformRand(); |
---|
281 | targetOscillator = numberOfOscillators-1; //last level |
---|
282 | S=0.0; |
---|
283 | G4bool levelFound = false; |
---|
284 | for (size_t j=0;j<numberOfOscillators && !levelFound; j++) |
---|
285 | { |
---|
286 | S += (*theTable)[j]->GetOscillatorStrength(); |
---|
287 | if (S > TST) |
---|
288 | { |
---|
289 | targetOscillator = j; |
---|
290 | levelFound = true; |
---|
291 | } |
---|
292 | } |
---|
293 | //check whether the level is valid |
---|
294 | ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy(); |
---|
295 | }while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0); |
---|
296 | } |
---|
297 | else //photonEnergy0 < 5 MeV |
---|
298 | { |
---|
299 | //Incoherent scattering function for theta=PI |
---|
300 | G4double s0=0.0; |
---|
301 | G4double pzomc=0.0; |
---|
302 | G4double rni=0.0; |
---|
303 | G4double aux=0.0; |
---|
304 | for (size_t i=0;i<numberOfOscillators;i++) |
---|
305 | { |
---|
306 | ionEnergy = (*theTable)[i]->GetIonisationEnergy(); |
---|
307 | if (photonEnergy0 > ionEnergy) |
---|
308 | { |
---|
309 | G4double aux = photonEnergy0*(photonEnergy0-ionEnergy)*2.0; |
---|
310 | hartreeFunc = (*theTable)[i]->GetHartreeFactor(); |
---|
311 | oscStren = (*theTable)[i]->GetOscillatorStrength(); |
---|
312 | pzomc = hartreeFunc*(aux-electron_mass_c2*ionEnergy)/ |
---|
313 | (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy)); |
---|
314 | if (pzomc > 0) |
---|
315 | rni = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)* |
---|
316 | (std::sqrt(0.5)+std::sqrt(2.0)*pzomc)); |
---|
317 | else |
---|
318 | rni = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)* |
---|
319 | (std::sqrt(0.5)-std::sqrt(2.0)*pzomc)); |
---|
320 | s0 += oscStren*rni; |
---|
321 | } |
---|
322 | } |
---|
323 | //Sampling tau |
---|
324 | G4double cdt1 = 0.; |
---|
325 | do |
---|
326 | { |
---|
327 | if ((G4UniformRand()*a2) < a1) |
---|
328 | tau = std::pow(taumin,G4UniformRand()); |
---|
329 | else |
---|
330 | tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0)); |
---|
331 | cdt1 = (1.0-tau)/(ek*tau); |
---|
332 | //Incoherent scattering function |
---|
333 | S = 0.; |
---|
334 | for (size_t i=0;i<numberOfOscillators;i++) |
---|
335 | { |
---|
336 | ionEnergy = (*theTable)[i]->GetIonisationEnergy(); |
---|
337 | if (photonEnergy0 > ionEnergy) //sum only on excitable levels |
---|
338 | { |
---|
339 | aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1; |
---|
340 | hartreeFunc = (*theTable)[i]->GetHartreeFactor(); |
---|
341 | oscStren = (*theTable)[i]->GetOscillatorStrength(); |
---|
342 | pzomc = hartreeFunc*(aux-electron_mass_c2*ionEnergy)/ |
---|
343 | (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy)); |
---|
344 | if (pzomc > 0) |
---|
345 | rn[i] = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)* |
---|
346 | (std::sqrt(0.5)+std::sqrt(2.0)*pzomc)); |
---|
347 | else |
---|
348 | rn[i] = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)* |
---|
349 | (std::sqrt(0.5)-std::sqrt(2.0)*pzomc)); |
---|
350 | S += oscStren*rn[i]; |
---|
351 | pac[i] = S; |
---|
352 | } |
---|
353 | else |
---|
354 | pac[i] = S-1e-6; |
---|
355 | } |
---|
356 | //Rejection function |
---|
357 | TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau)); |
---|
358 | }while ((G4UniformRand()*s0) > TST); |
---|
359 | |
---|
360 | cosTheta = 1.0 - cdt1; |
---|
361 | G4double fpzmax=0.0,fpz=0.0; |
---|
362 | G4double A=0.0; |
---|
363 | //Target electron shell |
---|
364 | do |
---|
365 | { |
---|
366 | do |
---|
367 | { |
---|
368 | TST = S*G4UniformRand(); |
---|
369 | targetOscillator = numberOfOscillators-1; //last level |
---|
370 | G4bool levelFound = false; |
---|
371 | for (size_t i=0;i<numberOfOscillators && !levelFound;i++) |
---|
372 | { |
---|
373 | if (pac[i]>TST) |
---|
374 | { |
---|
375 | targetOscillator = i; |
---|
376 | levelFound = true; |
---|
377 | } |
---|
378 | } |
---|
379 | A = G4UniformRand()*rn[targetOscillator]; |
---|
380 | hartreeFunc = (*theTable)[targetOscillator]->GetHartreeFactor(); |
---|
381 | oscStren = (*theTable)[targetOscillator]->GetOscillatorStrength(); |
---|
382 | if (A < 0.5) |
---|
383 | pzomc = (std::sqrt(0.5)-std::sqrt(0.5-std::log(2.0*A)))/ |
---|
384 | (std::sqrt(2.0)*hartreeFunc); |
---|
385 | else |
---|
386 | pzomc = (std::sqrt(0.5-std::log(2.0-2.0*A))-std::sqrt(0.5))/ |
---|
387 | (std::sqrt(2.0)*hartreeFunc); |
---|
388 | } while (pzomc < -1); |
---|
389 | |
---|
390 | // F(EP) rejection |
---|
391 | G4double XQC = 1.0+tau*(tau-2.0*cosTheta); |
---|
392 | G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC); |
---|
393 | if (AF > 0) |
---|
394 | fpzmax = 1.0+AF*0.2; |
---|
395 | else |
---|
396 | fpzmax = 1.0-AF*0.2; |
---|
397 | fpz = 1.0+AF*std::max(std::min(pzomc,0.2),-0.2); |
---|
398 | }while ((fpzmax*G4UniformRand())>fpz); |
---|
399 | |
---|
400 | //Energy of the scattered photon |
---|
401 | G4double T = pzomc*pzomc; |
---|
402 | G4double b1 = 1.0-T*tau*tau; |
---|
403 | G4double b2 = 1.0-T*tau*cosTheta; |
---|
404 | if (pzomc > 0.0) |
---|
405 | epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T)))); |
---|
406 | else |
---|
407 | epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T)))); |
---|
408 | } //energy < 5 MeV |
---|
409 | |
---|
410 | //Ok, the kinematics has been calculated. |
---|
411 | G4double sinTheta = std::sqrt(1-cosTheta*cosTheta); |
---|
412 | G4double phi = twopi * G4UniformRand() ; |
---|
413 | G4double dirx = sinTheta * std::cos(phi); |
---|
414 | G4double diry = sinTheta * std::sin(phi); |
---|
415 | G4double dirz = cosTheta ; |
---|
416 | |
---|
417 | // Update G4VParticleChange for the scattered photon |
---|
418 | G4ThreeVector photonDirection1(dirx,diry,dirz); |
---|
419 | photonDirection1.rotateUz(photonDirection0); |
---|
420 | fParticleChange->ProposeMomentumDirection(photonDirection1) ; |
---|
421 | |
---|
422 | G4double photonEnergy1 = epsilon * photonEnergy0; |
---|
423 | |
---|
424 | if (photonEnergy1 > 0.) |
---|
425 | fParticleChange->SetProposedKineticEnergy(photonEnergy1) ; |
---|
426 | else |
---|
427 | { |
---|
428 | fParticleChange->SetProposedKineticEnergy(0.) ; |
---|
429 | fParticleChange->ProposeTrackStatus(fStopAndKill); |
---|
430 | } |
---|
431 | |
---|
432 | //Create scattered electron |
---|
433 | G4double diffEnergy = photonEnergy0*(1-epsilon); |
---|
434 | ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy(); |
---|
435 | |
---|
436 | G4double Q2 = |
---|
437 | photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta); |
---|
438 | G4double cosThetaE = 0.; //scattering angle for the electron |
---|
439 | |
---|
440 | if (Q2 > 1.0e-12) |
---|
441 | cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2); |
---|
442 | else |
---|
443 | cosThetaE = 1.0; |
---|
444 | G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE); |
---|
445 | |
---|
446 | //Now, try to handle fluorescence |
---|
447 | //Notice: merged levels are indicated with Z=0 and flag=30 |
---|
448 | G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag(); |
---|
449 | G4int Z = (G4int) (*theTable)[targetOscillator]->GetParentZ(); |
---|
450 | |
---|
451 | //initialize here, then check photons created by Atomic-Deexcitation, and the final state e- |
---|
452 | std::vector<G4DynamicParticle*>* photonVector=0; |
---|
453 | const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance(); |
---|
454 | G4double bindingEnergy = 0.*eV; |
---|
455 | G4int shellId = 0; |
---|
456 | |
---|
457 | //Real level |
---|
458 | if (Z > 0 && shFlag<30) |
---|
459 | { |
---|
460 | const G4AtomicShell* shell = transitionManager->Shell(Z,shFlag-1); |
---|
461 | bindingEnergy = shell->BindingEnergy(); |
---|
462 | shellId = shell->ShellId(); |
---|
463 | } |
---|
464 | |
---|
465 | G4double ionEnergyInPenelopeDatabase = ionEnergy; |
---|
466 | //protection against energy non-conservation |
---|
467 | ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase); |
---|
468 | |
---|
469 | //subtract the excitation energy. If not emitted by fluorescence |
---|
470 | //the ionization energy is deposited as local energy deposition |
---|
471 | G4double eKineticEnergy = diffEnergy - ionEnergy; |
---|
472 | G4double localEnergyDeposit = ionEnergy; |
---|
473 | G4double energyInFluorescence = 0.; //testing purposes only |
---|
474 | |
---|
475 | if (eKineticEnergy < 0) |
---|
476 | { |
---|
477 | //It means that there was some problem/mismatch between the two databases. |
---|
478 | //Try to make it work |
---|
479 | //In this case available Energy (diffEnergy) < ionEnergy |
---|
480 | //Full residual energy is deposited locally |
---|
481 | localEnergyDeposit = diffEnergy; |
---|
482 | eKineticEnergy = 0.0; |
---|
483 | } |
---|
484 | |
---|
485 | //the local energy deposit is what remains: part of this may be spent for fluorescence. |
---|
486 | if(DeexcitationFlag() && Z > 5) { |
---|
487 | |
---|
488 | const G4ProductionCutsTable* theCoupleTable= |
---|
489 | G4ProductionCutsTable::GetProductionCutsTable(); |
---|
490 | |
---|
491 | size_t index = couple->GetIndex(); |
---|
492 | G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index]; |
---|
493 | G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index]; |
---|
494 | |
---|
495 | // Generation of fluorescence |
---|
496 | // Data in EADL are available only for Z > 5 |
---|
497 | // Protection to avoid generating photons in the unphysical case of |
---|
498 | // shell binding energy > photon energy |
---|
499 | if (localEnergyDeposit > cutg || localEnergyDeposit > cute) |
---|
500 | { |
---|
501 | G4DynamicParticle* aPhoton; |
---|
502 | deexcitationManager.SetCutForSecondaryPhotons(cutg); |
---|
503 | deexcitationManager.SetCutForAugerElectrons(cute); |
---|
504 | |
---|
505 | photonVector = deexcitationManager.GenerateParticles(Z,shellId); |
---|
506 | if(photonVector) |
---|
507 | { |
---|
508 | size_t nPhotons = photonVector->size(); |
---|
509 | for (size_t k=0; k<nPhotons; k++) |
---|
510 | { |
---|
511 | aPhoton = (*photonVector)[k]; |
---|
512 | if (aPhoton) |
---|
513 | { |
---|
514 | G4double itsEnergy = aPhoton->GetKineticEnergy(); |
---|
515 | if (itsEnergy <= localEnergyDeposit) |
---|
516 | { |
---|
517 | localEnergyDeposit -= itsEnergy; |
---|
518 | if (aPhoton->GetDefinition() == G4Gamma::Gamma()) |
---|
519 | energyInFluorescence += itsEnergy;; |
---|
520 | fvect->push_back(aPhoton); |
---|
521 | } |
---|
522 | else |
---|
523 | { |
---|
524 | delete aPhoton; |
---|
525 | (*photonVector)[k]=0; |
---|
526 | } |
---|
527 | } |
---|
528 | } |
---|
529 | delete photonVector; |
---|
530 | } |
---|
531 | } |
---|
532 | } |
---|
533 | |
---|
534 | //Always produce explicitely the electron |
---|
535 | G4DynamicParticle* electron = 0; |
---|
536 | |
---|
537 | G4double xEl = sinThetaE * std::cos(phi+pi); |
---|
538 | G4double yEl = sinThetaE * std::sin(phi+pi); |
---|
539 | G4double zEl = cosThetaE; |
---|
540 | G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction |
---|
541 | eDirection.rotateUz(photonDirection0); |
---|
542 | electron = new G4DynamicParticle (G4Electron::Electron(), |
---|
543 | eDirection,eKineticEnergy) ; |
---|
544 | fvect->push_back(electron); |
---|
545 | |
---|
546 | |
---|
547 | if (localEnergyDeposit < 0) |
---|
548 | { |
---|
549 | G4cout << "WARNING-" |
---|
550 | << "G4Penelope08ComptonModel::SampleSecondaries - Negative energy deposit" |
---|
551 | << G4endl; |
---|
552 | localEnergyDeposit=0.; |
---|
553 | } |
---|
554 | fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit); |
---|
555 | |
---|
556 | G4double electronEnergy = 0.; |
---|
557 | if (verboseLevel > 1) |
---|
558 | { |
---|
559 | G4cout << "-----------------------------------------------------------" << G4endl; |
---|
560 | G4cout << "Energy balance from G4Penelope08Compton" << G4endl; |
---|
561 | G4cout << "Incoming photon energy: " << photonEnergy0/keV << " keV" << G4endl; |
---|
562 | G4cout << "-----------------------------------------------------------" << G4endl; |
---|
563 | G4cout << "Scattered photon: " << photonEnergy1/keV << " keV" << G4endl; |
---|
564 | if (electron) |
---|
565 | electronEnergy = eKineticEnergy; |
---|
566 | G4cout << "Scattered electron " << electronEnergy/keV << " keV" << G4endl; |
---|
567 | G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl; |
---|
568 | G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl; |
---|
569 | G4cout << "Total final state: " << (photonEnergy1+electronEnergy+energyInFluorescence+ |
---|
570 | localEnergyDeposit)/keV << |
---|
571 | " keV" << G4endl; |
---|
572 | G4cout << "-----------------------------------------------------------" << G4endl; |
---|
573 | } |
---|
574 | if (verboseLevel > 0) |
---|
575 | { |
---|
576 | G4double energyDiff = std::fabs(photonEnergy1+ |
---|
577 | electronEnergy+energyInFluorescence+ |
---|
578 | localEnergyDeposit-photonEnergy0); |
---|
579 | if (energyDiff > 0.05*keV) |
---|
580 | G4cout << "Warning from G4Penelope08Compton: problem with energy conservation: " << |
---|
581 | (photonEnergy1+electronEnergy+energyInFluorescence+localEnergyDeposit)/keV << |
---|
582 | " keV (final) vs. " << |
---|
583 | photonEnergy0/keV << " keV (initial)" << G4endl; |
---|
584 | } |
---|
585 | } |
---|
586 | |
---|
587 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
588 | |
---|
589 | G4double G4Penelope08ComptonModel::DifferentialCrossSection(G4double cosTheta,G4double energy, |
---|
590 | G4PenelopeOscillator* osc) |
---|
591 | { |
---|
592 | // |
---|
593 | // Penelope model. Single differential cross section *per electron* |
---|
594 | // for photon Compton scattering by |
---|
595 | // electrons in the given atomic oscillator, differential in the direction of the |
---|
596 | // scattering photon. This is in units of pi*classic_electr_radius**2 |
---|
597 | // |
---|
598 | // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167 |
---|
599 | // The parametrization includes the J(p) distribution profiles for the atomic shells, |
---|
600 | // that are tabulated from Hartree-Fock calculations |
---|
601 | // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201 |
---|
602 | // |
---|
603 | G4double ionEnergy = osc->GetIonisationEnergy(); |
---|
604 | G4double harFunc = osc->GetHartreeFactor(); |
---|
605 | |
---|
606 | const G4double k2 = std::sqrt(2.); |
---|
607 | const G4double k1 = 1./k2; |
---|
608 | |
---|
609 | if (energy < ionEnergy) |
---|
610 | return 0; |
---|
611 | |
---|
612 | //energy of the Compton line |
---|
613 | G4double cdt1 = 1.0-cosTheta; |
---|
614 | G4double EOEC = 1.0+(energy/electron_mass_c2)*cdt1; |
---|
615 | G4double ECOE = 1.0/EOEC; |
---|
616 | |
---|
617 | //Incoherent scattering function (analytical profile) |
---|
618 | G4double aux = energy*(energy-ionEnergy)*cdt1; |
---|
619 | G4double Pzimax = |
---|
620 | (aux - electron_mass_c2*ionEnergy)/(electron_mass_c2*std::sqrt(2*aux+ionEnergy*ionEnergy)); |
---|
621 | G4double sia = 0.0; |
---|
622 | G4double x = harFunc*Pzimax; |
---|
623 | if (x > 0) |
---|
624 | sia = 1.0-0.5*std::exp(0.5-(k1+k2*x)*(k1+k2*x)); |
---|
625 | else |
---|
626 | sia = 0.5*std::exp(0.5-(k1-k2*x)*(k1-k2*x)); |
---|
627 | |
---|
628 | //1st order correction, integral of Pz times the Compton profile. |
---|
629 | //Calculated approximately using a free-electron gas profile |
---|
630 | G4double pf = 3.0/(4.0*harFunc); |
---|
631 | if (std::fabs(Pzimax) < pf) |
---|
632 | { |
---|
633 | G4double QCOE2 = 1.0+ECOE*ECOE-2.0*ECOE*cosTheta; |
---|
634 | G4double p2 = Pzimax*Pzimax; |
---|
635 | G4double dspz = std::sqrt(QCOE2)* |
---|
636 | (1.0+ECOE*(ECOE-cosTheta)/QCOE2)*harFunc |
---|
637 | *0.25*(2*p2-(p2*p2)/(pf*pf)-(pf*pf)); |
---|
638 | sia += std::max(dspz,-1.0*sia); |
---|
639 | } |
---|
640 | |
---|
641 | G4double XKN = EOEC+ECOE-1.0+cosTheta*cosTheta; |
---|
642 | |
---|
643 | //Differential cross section (per electron, in units of pi*classic_electr_radius**2) |
---|
644 | G4double diffCS = ECOE*ECOE*XKN*sia; |
---|
645 | |
---|
646 | return diffCS; |
---|
647 | } |
---|
648 | |
---|
649 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
650 | |
---|
651 | void G4Penelope08ComptonModel::ActivateAuger(G4bool augerbool) |
---|
652 | { |
---|
653 | if (!DeexcitationFlag() && augerbool) |
---|
654 | { |
---|
655 | G4cout << "WARNING - G4Penelope08ComptonModel" << G4endl; |
---|
656 | G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl; |
---|
657 | G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl; |
---|
658 | } |
---|
659 | deexcitationManager.ActivateAugerElectronProduction(augerbool); |
---|
660 | if (verboseLevel > 1) |
---|
661 | G4cout << "Auger production set to " << augerbool << G4endl; |
---|
662 | } |
---|
663 | |
---|
664 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
665 | |
---|
666 | G4double G4Penelope08ComptonModel::OscillatorTotalCrossSection(G4double energy,G4PenelopeOscillator* osc) |
---|
667 | { |
---|
668 | //Total cross section (integrated) for the given oscillator in units of |
---|
669 | //pi*classic_electr_radius^2 |
---|
670 | |
---|
671 | //Integrate differential cross section for each oscillator |
---|
672 | G4double stre = osc->GetOscillatorStrength(); |
---|
673 | |
---|
674 | // here one uses the using the 20-point |
---|
675 | // Gauss quadrature method with an adaptive bipartition scheme |
---|
676 | const G4int npoints=10; |
---|
677 | const G4int ncallsmax=20000; |
---|
678 | const G4int nst=256; |
---|
679 | static G4double Abscissas[10] = {7.652651133497334e-02,2.2778585114164508e-01,3.7370608871541956e-01, |
---|
680 | 5.1086700195082710e-01,6.3605368072651503e-01,7.4633190646015079e-01, |
---|
681 | 8.3911697182221882e-01,9.1223442825132591e-01,9.6397192727791379e-01, |
---|
682 | 9.9312859918509492e-01}; |
---|
683 | static G4double Weights[10] = {1.5275338713072585e-01,1.4917298647260375e-01,1.4209610931838205e-01, |
---|
684 | 1.3168863844917663e-01,1.1819453196151842e-01,1.0193011981724044e-01, |
---|
685 | 8.3276741576704749e-02,6.2672048334109064e-02,4.0601429800386941e-02, |
---|
686 | 1.7614007139152118e-02}; |
---|
687 | |
---|
688 | G4double MaxError = 1e-5; |
---|
689 | //Error control |
---|
690 | G4double Ctol = std::min(std::max(MaxError,1e-13),1e-02); |
---|
691 | G4double Ptol = 0.01*Ctol; |
---|
692 | G4double Err=1e35; |
---|
693 | |
---|
694 | //Gauss integration from -1 to 1 |
---|
695 | G4double LowPoint = -1.0; |
---|
696 | G4double HighPoint = 1.0; |
---|
697 | |
---|
698 | G4double h=HighPoint-LowPoint; |
---|
699 | G4double sumga=0.0; |
---|
700 | G4double a=0.5*(HighPoint-LowPoint); |
---|
701 | G4double b=0.5*(HighPoint+LowPoint); |
---|
702 | G4double c=a*Abscissas[0]; |
---|
703 | G4double d= Weights[0]* |
---|
704 | (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc)); |
---|
705 | for (G4int i=2;i<=npoints;i++) |
---|
706 | { |
---|
707 | c=a*Abscissas[i-1]; |
---|
708 | d += Weights[i-1]* |
---|
709 | (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc)); |
---|
710 | } |
---|
711 | G4int icall = 2*npoints; |
---|
712 | G4int LH=1; |
---|
713 | G4double S[nst],x[nst],sn[nst],xrn[nst]; |
---|
714 | S[0]=d*a; |
---|
715 | x[0]=LowPoint; |
---|
716 | |
---|
717 | G4bool loopAgain = true; |
---|
718 | |
---|
719 | //Adaptive bipartition scheme |
---|
720 | do{ |
---|
721 | G4double h0=h; |
---|
722 | h=0.5*h; //bipartition |
---|
723 | G4double sumr=0; |
---|
724 | G4int LHN=0; |
---|
725 | G4double si,xa,xb,xc; |
---|
726 | for (G4int i=1;i<=LH;i++){ |
---|
727 | si=S[i-1]; |
---|
728 | xa=x[i-1]; |
---|
729 | xb=xa+h; |
---|
730 | xc=xa+h0; |
---|
731 | a=0.5*(xb-xa); |
---|
732 | b=0.5*(xb+xa); |
---|
733 | c=a*Abscissas[0]; |
---|
734 | G4double d = Weights[0]* |
---|
735 | (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc)); |
---|
736 | |
---|
737 | for (G4int j=1;j<npoints;j++) |
---|
738 | { |
---|
739 | c=a*Abscissas[j]; |
---|
740 | d += Weights[j]* |
---|
741 | (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc)); |
---|
742 | } |
---|
743 | G4double s1=d*a; |
---|
744 | a=0.5*(xc-xb); |
---|
745 | b=0.5*(xc+xb); |
---|
746 | c=a*Abscissas[0]; |
---|
747 | d=Weights[0]* |
---|
748 | (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc)); |
---|
749 | |
---|
750 | for (G4int j=1;j<npoints;j++) |
---|
751 | { |
---|
752 | c=a*Abscissas[j]; |
---|
753 | d += Weights[j]* |
---|
754 | (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc)); |
---|
755 | } |
---|
756 | G4double s2=d*a; |
---|
757 | icall=icall+4*npoints; |
---|
758 | G4double s12=s1+s2; |
---|
759 | if (std::abs(s12-si)<std::max(Ptol*std::abs(s12),1e-35)) |
---|
760 | sumga += s12; |
---|
761 | else |
---|
762 | { |
---|
763 | sumr += s12; |
---|
764 | LHN += 2; |
---|
765 | sn[LHN-1]=s2; |
---|
766 | xrn[LHN-1]=xb; |
---|
767 | sn[LHN-2]=s1; |
---|
768 | xrn[LHN-2]=xa; |
---|
769 | } |
---|
770 | |
---|
771 | if (icall>ncallsmax || LHN>nst) |
---|
772 | { |
---|
773 | G4cout << "G4Penelope08ComptonModel: " << G4endl; |
---|
774 | G4cout << "LowPoint: " << LowPoint << ", High Point: " << HighPoint << G4endl; |
---|
775 | G4cout << "Tolerance: " << MaxError << G4endl; |
---|
776 | G4cout << "Calls: " << icall << ", Integral: " << sumga << ", Error: " << Err << G4endl; |
---|
777 | G4cout << "Number of open subintervals: " << LHN << G4endl; |
---|
778 | G4cout << "WARNING: the required accuracy has not been attained" << G4endl; |
---|
779 | loopAgain = false; |
---|
780 | } |
---|
781 | } |
---|
782 | Err=std::abs(sumr)/std::max(std::abs(sumr+sumga),1e-35); |
---|
783 | if (Err < Ctol || LHN == 0) |
---|
784 | loopAgain = false; //end of cycle |
---|
785 | LH=LHN; |
---|
786 | for (G4int i=0;i<LH;i++) |
---|
787 | { |
---|
788 | S[i]=sn[i]; |
---|
789 | x[i]=xrn[i]; |
---|
790 | } |
---|
791 | }while(Ctol < 1.0 && loopAgain); |
---|
792 | |
---|
793 | |
---|
794 | G4double xs = stre*sumga; |
---|
795 | |
---|
796 | return xs; |
---|
797 | } |
---|
798 | |
---|
799 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
800 | |
---|
801 | G4double G4Penelope08ComptonModel::KleinNishinaCrossSection(G4double energy, |
---|
802 | const G4Material* material) |
---|
803 | { |
---|
804 | // use Klein-Nishina formula |
---|
805 | // total cross section in units of pi*classic_electr_radius^2 |
---|
806 | |
---|
807 | G4double cs = 0; |
---|
808 | |
---|
809 | G4double ek =energy/electron_mass_c2; |
---|
810 | G4double eks = ek*ek; |
---|
811 | G4double ek2 = 1.0+ek+ek; |
---|
812 | G4double ek1 = eks-ek2-1.0; |
---|
813 | |
---|
814 | G4double t0 = 1.0/ek2; |
---|
815 | G4double csl = 0.5*eks*t0*t0+ek2*t0+ek1*std::log(t0)-(1.0/t0); |
---|
816 | |
---|
817 | G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material); |
---|
818 | |
---|
819 | for (size_t i=0;i<theTable->size();i++) |
---|
820 | { |
---|
821 | G4PenelopeOscillator* theOsc = (*theTable)[i]; |
---|
822 | G4double ionEnergy = theOsc->GetIonisationEnergy(); |
---|
823 | G4double tau=(energy-ionEnergy)/energy; |
---|
824 | if (tau > t0) |
---|
825 | { |
---|
826 | G4double csu = 0.5*eks*tau*tau+ek2*tau+ek1*std::log(tau)-(1.0/tau); |
---|
827 | G4double stre = theOsc->GetOscillatorStrength(); |
---|
828 | |
---|
829 | cs += stre*(csu-csl); |
---|
830 | } |
---|
831 | } |
---|
832 | |
---|
833 | cs /= (ek*eks); |
---|
834 | |
---|
835 | return cs; |
---|
836 | |
---|
837 | } |
---|
838 | |
---|