source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeComptonModel.cc

Last change on this file was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 27.3 KB
Line 
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: G4PenelopeComptonModel.cc,v 1.11 2010/12/01 15:20:26 pandola Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 02 Oct 2008 L Pandola Migration from process to model
34// 28 Oct 2008 L Pandola Treat the database data from Penelope according to the
35// original model, namely merging levels below 15 eV in
36// a single one. Still, it is not fully compliant with the
37// original Penelope model, because plasma excitation is not
38// considered.
39// 22 Nov 2008 L Pandola Make unit of measurements explicit for binding energies
40// that are read from the external files.
41// 24 Nov 2008 L Pandola Find a cleaner way to delete vectors.
42// 16 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
43// - apply internal high-energy limit only in constructor
44// - do not apply low-energy limit (default is 0)
45// - do not apply production threshold on level of the model
46// 21 Oct 2009 L Pandola Remove un-necessary fUseAtomicDeexcitation flag - now managed by
47// G4VEmModel::DeexcitationFlag()
48// Add ActivateAuger() method
49//
50
51#include "G4PenelopeComptonModel.hh"
52#include "G4ParticleDefinition.hh"
53#include "G4MaterialCutsCouple.hh"
54#include "G4ProductionCutsTable.hh"
55#include "G4DynamicParticle.hh"
56#include "G4VEMDataSet.hh"
57#include "G4PhysicsTable.hh"
58#include "G4ElementTable.hh"
59#include "G4Element.hh"
60#include "G4PhysicsLogVector.hh"
61#include "G4PenelopeIntegrator.hh"
62#include "G4AtomicTransitionManager.hh"
63#include "G4AtomicShell.hh"
64#include "G4Gamma.hh"
65#include "G4Electron.hh"
66
67//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68
69
70G4PenelopeComptonModel::G4PenelopeComptonModel(const G4ParticleDefinition*,
71 const G4String& nam)
72 :G4VEmModel(nam),ionizationEnergy(0),hartreeFunction(0),
73 occupationNumber(0),isInitialised(false)
74{
75 fIntrinsicLowEnergyLimit = 100.0*eV;
76 fIntrinsicHighEnergyLimit = 100.0*GeV;
77 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
78 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
79 //
80 energyForIntegration = 0.0;
81 ZForIntegration = 1;
82
83 verboseLevel= 0;
84 // Verbosity scale:
85 // 0 = nothing
86 // 1 = warning for energy non-conservation
87 // 2 = details of energy budget
88 // 3 = calculation of cross sections, file openings, sampling of atoms
89 // 4 = entering in methods
90
91 //by default, the model will use atomic deexcitation
92 SetDeexcitationFlag(true);
93 ActivateAuger(false);
94
95 //These vectors do not change when materials or cut change.
96 //Therefore I can read it at the constructor
97 ionizationEnergy = new std::map<G4int,G4DataVector*>;
98 hartreeFunction = new std::map<G4int,G4DataVector*>;
99 occupationNumber = new std::map<G4int,G4DataVector*>;
100
101 ReadData(); //Read data from file
102
103}
104
105//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
106
107G4PenelopeComptonModel::~G4PenelopeComptonModel()
108{
109 std::map <G4int,G4DataVector*>::iterator i;
110 if (ionizationEnergy)
111 {
112 for (i=ionizationEnergy->begin();i != ionizationEnergy->end();i++)
113 if (i->second) delete i->second;
114 delete ionizationEnergy;
115 }
116 if (hartreeFunction)
117 {
118 for (i=hartreeFunction->begin();i != hartreeFunction->end();i++)
119 if (i->second) delete i->second;
120 delete hartreeFunction;
121 }
122 if (occupationNumber)
123 {
124 for (i=occupationNumber->begin();i != occupationNumber->end();i++)
125 if (i->second) delete i->second;
126 delete occupationNumber;
127 }
128}
129
130//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131
132void G4PenelopeComptonModel::Initialise(const G4ParticleDefinition* particle,
133 const G4DataVector& cuts)
134{
135 if (verboseLevel > 3)
136 G4cout << "Calling G4PenelopeComptonModel::Initialise()" << G4endl;
137
138 InitialiseElementSelectors(particle,cuts);
139
140 if (verboseLevel > 0) {
141 G4cout << "Penelope Compton model is initialized " << G4endl
142 << "Energy range: "
143 << LowEnergyLimit() / keV << " keV - "
144 << HighEnergyLimit() / GeV << " GeV"
145 << G4endl;
146 }
147
148 if(isInitialised) return;
149 fParticleChange = GetParticleChangeForGamma();
150 isInitialised = true;
151}
152
153//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
154
155G4double G4PenelopeComptonModel::ComputeCrossSectionPerAtom(
156 const G4ParticleDefinition*,
157 G4double energy,
158 G4double Z, G4double,
159 G4double, G4double)
160{
161 // Penelope model to calculate the Compton scattering cross section:
162 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
163 //
164 // The cross section for Compton scattering is calculated according to the Klein-Nishina
165 // formula for energy > 5 MeV.
166 // For E < 5 MeV it is used a parametrization for the differential cross-section dSigma/dOmega,
167 // which is integrated numerically in cos(theta), G4PenelopeComptonModel::DifferentialCrossSection().
168 // The parametrization includes the J(p)
169 // distribution profiles for the atomic shells, that are tabulated from Hartree-Fock calculations
170 // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
171 //
172 if (verboseLevel > 3)
173 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeComptonModel" << G4endl;
174
175 G4int iZ = (G4int) Z;
176 G4double cs=0.0;
177 energyForIntegration=energy;
178 ZForIntegration = iZ;
179 if (energy< 5*MeV)
180 {
181 // numerically integrate differential cross section dSigma/dOmega
182 G4PenelopeIntegrator<G4PenelopeComptonModel,G4double (G4PenelopeComptonModel::*)(G4double)>
183 theIntegrator;
184 cs = theIntegrator.Calculate(this,&G4PenelopeComptonModel::DifferentialCrossSection,-1.0,1.0,1e-05);
185 }
186 else
187 {
188 // use Klein-Nishina formula
189 G4double ki=energy/electron_mass_c2;
190 G4double ki3=ki*ki;
191 G4double ki2=1.0+2*ki;
192 G4double ki1=ki3-ki2-1.0;
193 G4double t0=1.0/(ki2);
194 G4double csl = 0.5*ki3*t0*t0+ki2*t0+ki1*std::log(t0)-(1.0/t0);
195 G4int nosc = occupationNumber->find(iZ)->second->size();
196 for (G4int i=0;i<nosc;i++)
197 {
198 G4double ionEnergy = (*(ionizationEnergy->find(iZ)->second))[i];
199 G4double tau=(energy-ionEnergy)/energy;
200 if (tau > t0)
201 {
202 G4double csu = 0.5*ki3*tau*tau+ki2*tau+ki1*std::log(tau)-(1.0/tau);
203 G4int f = (G4int) (*(occupationNumber->find(iZ)->second))[i];
204 cs = cs + f*(csu-csl);
205 }
206 }
207 cs=pi*classic_electr_radius*classic_electr_radius*cs/(ki*ki3);
208 }
209
210 if (verboseLevel > 2)
211 G4cout << "Compton cross Section at " << energy/keV << " keV for Z=" << Z <<
212 " = " << cs/barn << " barn" << G4endl;
213 return cs;
214}
215
216//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
217
218void G4PenelopeComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
219 const G4MaterialCutsCouple* couple,
220 const G4DynamicParticle* aDynamicGamma,
221 G4double,
222 G4double)
223{
224
225 // Penelope model to sample the Compton scattering final state.
226 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
227 // The model determines also the original shell from which the electron is expelled,
228 // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
229 //
230 // The final state for Compton scattering is calculated according to the Klein-Nishina
231 // formula for energy > 5 MeV. In this case, the Doppler broadening is negligible and
232 // one can assume that the target electron is at rest.
233 // For E < 5 MeV it is used the parametrization for the differential cross-section dSigma/dOmega,
234 // to sample the scattering angle and the energy of the emerging electron, which is
235 // G4PenelopeComptonModel::DifferentialCrossSection(). The rejection method is
236 // used to sample cos(theta). The efficiency increases monotonically with photon energy and is
237 // nearly independent on the Z; typical values are 35%, 80% and 95% for 1 keV, 1 MeV and 10 MeV,
238 // respectively.
239 // The parametrization includes the J(p) distribution profiles for the atomic shells, that are
240 // tabulated
241 // from Hartree-Fock calculations from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201.
242 // Doppler broadening is included.
243 //
244
245 if (verboseLevel > 3)
246 G4cout << "Calling SampleSecondaries() of G4PenelopeComptonModel" << G4endl;
247
248 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
249
250 if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
251 {
252 fParticleChange->ProposeTrackStatus(fStopAndKill);
253 fParticleChange->SetProposedKineticEnergy(0.);
254 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
255 return ;
256 }
257
258 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
259
260 // Select randomly one element in the current material
261 if (verboseLevel > 2)
262 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
263 // atom can be selected effitiantly if element selectors are initialised
264 const G4Element* anElement =
265 SelectRandomAtom(couple,G4Gamma::GammaDefinition(),photonEnergy0);
266 G4int Z = (G4int) anElement->GetZ();
267 if (verboseLevel > 2)
268 G4cout << "Selected " << anElement->GetName() << G4endl;
269
270 const G4int nmax = 64;
271 G4double rn[nmax]={0.0};
272 G4double pac[nmax]={0.0};
273
274 G4double ki,ki1,ki2,ki3,taumin,a1,a2;
275 G4double tau,TST;
276 G4double S=0.0;
277 G4double epsilon,cosTheta;
278 G4double harFunc = 0.0;
279 G4int occupNb= 0;
280 G4double ionEnergy=0.0;
281 G4int nosc = occupationNumber->find(Z)->second->size();
282 G4int iosc = nosc;
283 ki = photonEnergy0/electron_mass_c2;
284 ki2 = 2*ki+1.0;
285 ki3 = ki*ki;
286 ki1 = ki3-ki2-1.0;
287 taumin = 1.0/ki2;
288 a1 = std::log(ki2);
289 a2 = a1+2.0*ki*(1.0+ki)/(ki2*ki2);
290 //If the incoming photon is above 5 MeV, the quicker approach based on the
291 //pure Klein-Nishina formula is used
292 if (photonEnergy0 > 5*MeV)
293 {
294 do{
295 do{
296 if ((a2*G4UniformRand()) < a1)
297 {
298 tau = std::pow(taumin,G4UniformRand());
299 }
300 else
301 {
302 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
303 }
304 //rejection function
305 TST = (1+tau*(ki1+tau*(ki2+tau*ki3)))/(ki3*tau*(1.0+tau*tau));
306 }while (G4UniformRand()> TST);
307 epsilon=tau;
308 cosTheta = 1.0 - (1.0-tau)/(ki*tau);
309 //Target shell electrons
310 TST = Z*G4UniformRand();
311 iosc = nosc;
312 S=0.0;
313 for (G4int j=0;j<nosc;j++)
314 {
315 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[j];
316 S = S + occupNb;
317 if (S > TST) iosc = j;
318 if (S > TST) break;
319 }
320 ionEnergy = (*(ionizationEnergy->find(Z)->second))[iosc];
321 }while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0);
322 }
323 else //photonEnergy0 < 5 MeV
324 {
325 //Incoherent scattering function for theta=PI
326 G4double s0=0.0;
327 G4double pzomc=0.0,rni=0.0;
328 G4double aux=0.0;
329 for (G4int i=0;i<nosc;i++){
330 ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
331 if (photonEnergy0 > ionEnergy)
332 {
333 G4double aux = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
334 harFunc = (*(hartreeFunction->find(Z)->second))[i]/fine_structure_const;
335 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
336 pzomc = harFunc*(aux-electron_mass_c2*ionEnergy)/
337 (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
338 if (pzomc > 0)
339 {
340 rni = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
341 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
342 }
343 else
344 {
345 rni = 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 }
348 s0 = s0 + occupNb*rni;
349 }
350 }
351
352 //Sampling tau
353 G4double cdt1;
354 do
355 {
356 if ((G4UniformRand()*a2) < a1)
357 {
358 tau = std::pow(taumin,G4UniformRand());
359 }
360 else
361 {
362 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
363 }
364 cdt1 = (1.0-tau)/(ki*tau);
365 S=0.0;
366 //Incoherent scattering function
367 for (G4int i=0;i<nosc;i++){
368 ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
369 if (photonEnergy0 > ionEnergy) //sum only on excitable levels
370 {
371 aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1;
372 harFunc = (*(hartreeFunction->find(Z)->second))[i]/fine_structure_const;
373 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
374 pzomc = harFunc*(aux-electron_mass_c2*ionEnergy)/
375 (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
376 if (pzomc > 0)
377 {
378 rn[i] = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
379 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
380 }
381 else
382 {
383 rn[i] = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
384 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
385 }
386 S = S + occupNb*rn[i];
387 pac[i] = S;
388 }
389 else
390 {
391 pac[i] = S-(1e-06);
392 }
393 }
394 //Rejection function
395 TST = S*(1.0+tau*(ki1+tau*(ki2+tau*ki3)))/(ki3*tau*(1.0+tau*tau));
396 }while ((G4UniformRand()*s0) > TST);
397 //Target electron shell
398 cosTheta = 1.0 - cdt1;
399 G4double fpzmax=0.0,fpz=0.0;
400 G4double A=0.0;
401 do
402 {
403 do
404 {
405 TST =S*G4UniformRand();
406 iosc=nosc;
407 for (G4int i=0;i<nosc;i++){
408 if (pac[i]>TST) iosc = i;
409 if (pac[i]>TST) break;
410 }
411 A = G4UniformRand()*rn[iosc];
412 harFunc = (*(hartreeFunction->find(Z)->second))[iosc]/fine_structure_const;
413 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[iosc];
414 if (A < 0.5) {
415 pzomc = (std::sqrt(0.5)-std::sqrt(0.5-std::log(2.0*A)))/
416 (std::sqrt(2.0)*harFunc);
417 }
418 else
419 {
420 pzomc = (std::sqrt(0.5-std::log(2.0-2.0*A))-std::sqrt(0.5))/
421 (std::sqrt(2.0)*harFunc);
422 }
423 } while (pzomc < -1);
424 // F(EP) rejection
425 G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
426 G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC);
427 if (AF > 0) {
428 fpzmax = 1.0+AF*0.2;
429 }
430 else
431 {
432 fpzmax = 1.0-AF*0.2;
433 }
434 fpz = 1.0+AF*std::max(std::min(pzomc,0.2),-0.2);
435 }while ((fpzmax*G4UniformRand())>fpz);
436
437 //Energy of the scattered photon
438 G4double T = pzomc*pzomc;
439 G4double b1 = 1.0-T*tau*tau;
440 G4double b2 = 1.0-T*tau*cosTheta;
441 if (pzomc > 0.0)
442 {
443 epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
444 }
445 else
446 {
447 epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
448 }
449 }
450
451
452 //Ok, the kinematics has been calculated.
453 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
454 G4double phi = twopi * G4UniformRand() ;
455 G4double dirx = sinTheta * std::cos(phi);
456 G4double diry = sinTheta * std::sin(phi);
457 G4double dirz = cosTheta ;
458
459 // Update G4VParticleChange for the scattered photon
460 G4ThreeVector photonDirection1(dirx,diry,dirz);
461 photonDirection1.rotateUz(photonDirection0);
462 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
463 G4double photonEnergy1 = epsilon * photonEnergy0;
464
465 if (photonEnergy1 > 0.)
466 {
467 fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
468 }
469 else
470 {
471 fParticleChange->SetProposedKineticEnergy(0.) ;
472 fParticleChange->ProposeTrackStatus(fStopAndKill);
473 }
474
475 //Create scattered electron
476 G4double diffEnergy = photonEnergy0*(1-epsilon);
477 ionEnergy = (*(ionizationEnergy->find(Z)->second))[iosc];
478 G4double Q2 =
479 photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
480 G4double cosThetaE; //scattering angle for the electron
481 if (Q2 > 1.0e-12)
482 {
483 cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2);
484 }
485 else
486 {
487 cosThetaE = 1.0;
488 }
489 G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE);
490
491 //initialize here, then check photons created by Atomic-Deexcitation, and the final state e-
492 std::vector<G4DynamicParticle*>* photonVector=0;
493
494 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
495 const G4AtomicShell* shell = transitionManager->Shell(Z,iosc);
496 G4double bindingEnergy = shell->BindingEnergy();
497 G4int shellId = shell->ShellId();
498 G4double ionEnergyInPenelopeDatabase = ionEnergy;
499 //protection against energy non-conservation
500 ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase);
501
502 //subtract the excitation energy. If not emitted by fluorescence
503 //the ionization energy is deposited as local energy deposition
504 G4double eKineticEnergy = diffEnergy - ionEnergy;
505 G4double localEnergyDeposit = ionEnergy;
506 G4double energyInFluorescence = 0.; //testing purposes only
507
508 if (eKineticEnergy < 0)
509 {
510 //It means that there was some problem/mismatch between the two databases.
511 //Try to make it work
512 //In this case available Energy (diffEnergy) < ionEnergy
513 //Full residual energy is deposited locally
514 localEnergyDeposit = diffEnergy;
515 eKineticEnergy = 0.0;
516 }
517
518 //the local energy deposit is what remains: part of this may be spent for fluorescence.
519 if(DeexcitationFlag() && Z > 5) {
520
521 const G4ProductionCutsTable* theCoupleTable=
522 G4ProductionCutsTable::GetProductionCutsTable();
523
524 size_t index = couple->GetIndex();
525 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
526 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
527
528 // Generation of fluorescence
529 // Data in EADL are available only for Z > 5
530 // Protection to avoid generating photons in the unphysical case of
531 // shell binding energy > photon energy
532 if (localEnergyDeposit > cutg || localEnergyDeposit > cute)
533 {
534 G4DynamicParticle* aPhoton;
535 deexcitationManager.SetCutForSecondaryPhotons(cutg);
536 deexcitationManager.SetCutForAugerElectrons(cute);
537
538 photonVector = deexcitationManager.GenerateParticles(Z,shellId);
539 if(photonVector)
540 {
541 size_t nPhotons = photonVector->size();
542 for (size_t k=0; k<nPhotons; k++)
543 {
544 aPhoton = (*photonVector)[k];
545 if (aPhoton)
546 {
547 G4double itsEnergy = aPhoton->GetKineticEnergy();
548 if (itsEnergy <= localEnergyDeposit)
549 {
550 localEnergyDeposit -= itsEnergy;
551 if (aPhoton->GetDefinition() == G4Gamma::Gamma())
552 energyInFluorescence += itsEnergy;;
553 fvect->push_back(aPhoton);
554 }
555 else
556 {
557 delete aPhoton;
558 (*photonVector)[k]=0;
559 }
560 }
561 }
562 delete photonVector;
563 }
564 }
565 }
566
567 //Produce explicitely the electron only if its proposed kinetic energy is
568 //above the cut, otherwise add local energy deposition
569 G4DynamicParticle* electron = 0;
570 // if (eKineticEnergy > cutE) // VI: may be consistency problem if cut is applied here
571 if (eKineticEnergy > 0.0)
572 {
573 G4double xEl = sinThetaE * std::cos(phi+pi);
574 G4double yEl = sinThetaE * std::sin(phi+pi);
575 G4double zEl = cosThetaE;
576 G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
577 eDirection.rotateUz(photonDirection0);
578 electron = new G4DynamicParticle (G4Electron::Electron(),
579 eDirection,eKineticEnergy) ;
580 fvect->push_back(electron);
581 }
582 else
583 {
584 localEnergyDeposit += eKineticEnergy;
585 }
586
587 if (localEnergyDeposit < 0)
588 {
589 G4cout << "WARNING-"
590 << "G4PenelopeComptonModel::SampleSecondaries - Negative energy deposit"
591 << G4endl;
592 localEnergyDeposit=0.;
593 }
594 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
595
596 G4double electronEnergy = 0.;
597 if (verboseLevel > 1)
598 {
599 G4cout << "-----------------------------------------------------------" << G4endl;
600 G4cout << "Energy balance from G4PenelopeCompton" << G4endl;
601 G4cout << "Incoming photon energy: " << photonEnergy0/keV << " keV" << G4endl;
602 G4cout << "-----------------------------------------------------------" << G4endl;
603 G4cout << "Scattered photon: " << photonEnergy1/keV << " keV" << G4endl;
604 if (electron)
605 electronEnergy = eKineticEnergy;
606 G4cout << "Scattered electron " << electronEnergy/keV << " keV" << G4endl;
607 G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
608 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
609 G4cout << "Total final state: " << (photonEnergy1+electronEnergy+energyInFluorescence+
610 localEnergyDeposit)/keV <<
611 " keV" << G4endl;
612 G4cout << "-----------------------------------------------------------" << G4endl;
613 }
614 if (verboseLevel > 0)
615 {
616 G4double energyDiff = std::fabs(photonEnergy1+
617 electronEnergy+energyInFluorescence+
618 localEnergyDeposit-photonEnergy0);
619 if (energyDiff > 0.05*keV)
620 G4cout << "Warning from G4PenelopeCompton: problem with energy conservation: " <<
621 (photonEnergy1+electronEnergy+energyInFluorescence+localEnergyDeposit)/keV <<
622 " keV (final) vs. " <<
623 photonEnergy0/keV << " keV (initial)" << G4endl;
624 }
625}
626
627//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
628
629void G4PenelopeComptonModel::ReadData()
630{
631 char* path = getenv("G4LEDATA");
632 if (!path)
633 {
634 G4String excep = "G4PenelopeComptonModel - G4LEDATA environment variable not set!";
635 G4Exception(excep);
636 return;
637 }
638 G4String pathString(path);
639 G4String pathFile = pathString + "/penelope/compton-pen.dat";
640 std::ifstream file(pathFile);
641
642 if (!file.is_open())
643 {
644 G4String excep = "G4PenelopeComptonModel - data file " + pathFile + " not found!";
645 G4Exception(excep);
646 }
647
648 G4int k1,test,test1;
649 G4double a1,a2;
650 G4int Z=1,nLevels=0;
651
652 if (!ionizationEnergy || !hartreeFunction || !occupationNumber)
653 {
654 G4String excep = "G4PenelopeComptonModel: problem with reading data from file";
655 G4Exception(excep);
656 return;
657 }
658
659 do{
660 G4double harOfElectronsBelowThreshold = 0;
661 G4int nbOfElectronsBelowThreshold = 0;
662 file >> Z >> nLevels;
663 //Check for nLevels validity, before using it in a loop
664 if (nLevels<0 || nLevels>64)
665 {
666 G4String excep = "G4PenelopeComptonModel: corrupted data file?";
667 G4Exception(excep);
668 return;
669 }
670 G4DataVector* occVector = new G4DataVector;
671 G4DataVector* harVector = new G4DataVector;
672 G4DataVector* bindingEVector = new G4DataVector;
673 for (G4int h=0;h<nLevels;h++)
674 {
675 file >> k1 >> a1 >> a2;
676 //Make explicit unit of measurements for ionisation energy, which is MeV
677 a1 *= MeV;
678 if (a1 > 15*eV)
679 {
680 occVector->push_back((G4double) k1);
681 bindingEVector->push_back(a1);
682 harVector->push_back(a2);
683 }
684 else
685 {
686 nbOfElectronsBelowThreshold += k1;
687 harOfElectronsBelowThreshold += k1*a2;
688 }
689 }
690 //Add the "final" level
691 if (nbOfElectronsBelowThreshold)
692 {
693 occVector->push_back(nbOfElectronsBelowThreshold);
694 bindingEVector->push_back(0*eV);
695 G4double averageHartree =
696 harOfElectronsBelowThreshold/((G4double) nbOfElectronsBelowThreshold);
697 harVector->push_back(averageHartree);
698 }
699 //Ok, done for element Z
700 occupationNumber->insert(std::make_pair(Z,occVector));
701 ionizationEnergy->insert(std::make_pair(Z,bindingEVector));
702 hartreeFunction->insert(std::make_pair(Z,harVector));
703 file >> test >> test1; //-1 -1 close the data for each Z
704 if (test > 0) {
705 G4String excep = "G4PenelopeComptonModel - data file corrupted!";
706 G4Exception(excep);
707 }
708 }while (test != -2); //the very last Z is closed with -2 instead of -1
709 file.close();
710 if (verboseLevel > 2)
711 {
712 G4cout << "Data from G4PenelopeComptonModel read " << G4endl;
713 }
714}
715
716//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
717
718G4double G4PenelopeComptonModel::DifferentialCrossSection(G4double cosTheta)
719{
720 //
721 // Penelope model for the Compton scattering differential cross section
722 // dSigma/dOmega.
723 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
724 // The parametrization includes the J(p) distribution profiles for the atomic shells,
725 // that are tabulated from Hartree-Fock calculations
726 // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
727 //
728 const G4double k2 = std::sqrt(2.0);
729 const G4double k1 = std::sqrt(0.5);
730 const G4double k12 = 0.5;
731 G4double cdt1 = 1.0-cosTheta;
732 G4double energy = energyForIntegration;
733 G4int Z = ZForIntegration;
734 //energy of Compton line;
735 G4double EOEC = 1.0+(energy/electron_mass_c2)*cdt1;
736 G4double ECOE = 1.0/EOEC;
737 //Incoherent scattering function (analytical profile)
738 G4double sia = 0.0;
739 G4int nosc = occupationNumber->find(Z)->second->size();
740 for (G4int i=0;i<nosc;i++){
741 G4double ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
742 //Sum only of those shells for which E>Eion
743 if (energy > ionEnergy)
744 {
745 G4double aux = energy * (energy-ionEnergy)*cdt1;
746 G4double Pzimax =
747 (aux - electron_mass_c2*ionEnergy)/(electron_mass_c2*std::sqrt(2*aux+ionEnergy*ionEnergy));
748 G4double harFunc = (*(hartreeFunction->find(Z)->second))[i]/fine_structure_const;
749 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
750 G4double x = harFunc*Pzimax;
751 G4double siap = 0;
752 if (x > 0)
753 {
754 siap = 1.0-0.5*std::exp(k12-(k1+k2*x)*(k1+k2*x));
755 }
756 else
757 {
758 siap = 0.5*std::exp(k12-(k1-k2*x)*(k1-k2*x));
759 }
760 sia = sia + occupNb*siap; //sum of all contributions;
761 }
762 }
763 G4double XKN = EOEC+ECOE-1+cosTheta*cosTheta;
764 G4double diffCS = pi*classic_electr_radius*classic_electr_radius*ECOE*ECOE*XKN*sia;
765 return diffCS;
766}
767
768//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
769
770void G4PenelopeComptonModel::ActivateAuger(G4bool augerbool)
771{
772 if (!DeexcitationFlag() && augerbool)
773 {
774 G4cout << "WARNING - G4PenelopeComptonModel" << G4endl;
775 G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
776 G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
777 }
778 deexcitationManager.ActivateAugerElectronProduction(augerbool);
779 if (verboseLevel > 1)
780 G4cout << "Auger production set to " << augerbool << G4endl;
781}
782
783//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.