source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeCompton.cc@ 1201

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 26.7 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: G4PenelopeCompton.cc,v 1.36 2009/06/11 15:47:08 mantero Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 12 Feb 2003 MG Pia const argument in SelectRandomAtomForCompton
34// Migration to "cuts per region"
35// 14 Feb 2003 MG Pia Corrected compilation errors and warnings
36// from SUN
37// Modified some variables to lowercase initial
38// 10 Mar 2003 V.Ivanchenko Remove CutPerMaterial warning
39// 13 Mar 2003 L.Pandola Code "cleaned"
40// 20 Mar 2003 L.Pandola ReadData() changed (performance improved)
41// 26 Mar 2003 L.Pandola Added fluorescence
42// 24 May 2003 MGP Removed memory leak
43// 09 Mar 2004 L.Pandola Bug fixed in the generation of final state
44// (bug report # 585)
45// 17 Mar 2004 L.Pandola Removed unnecessary calls to std::pow(a,b)
46// 18 Mar 2004 L.Pandola Use of std::map (code review)
47// 26 Mar 2008 L.Pandola Add boolean flag to control atomic de-excitation
48// 27 Mar 2008 L.Pandola Re-named some variables to improve readability,
49// and check for strict energy conservation
50// 03 Jun 2008 L.Pandola Added further protection against non-conservation
51// of energy: it may happen because ionization energy
52// from de-excitation manager and from Penelope internal
53// database do not match (difference is <10 eV, but may
54// give a e- with negative kinetic energy).
55//
56// -------------------------------------------------------------------
57
58#include "G4PenelopeCompton.hh"
59#include "Randomize.hh"
60#include "G4ParticleDefinition.hh"
61#include "G4Track.hh"
62#include "G4Step.hh"
63#include "G4ForceCondition.hh"
64#include "G4Gamma.hh"
65#include "G4Electron.hh"
66#include "G4DynamicParticle.hh"
67#include "G4VParticleChange.hh"
68#include "G4ThreeVector.hh"
69#include "G4EnergyLossTables.hh"
70#include "G4VCrossSectionHandler.hh"
71#include "G4CrossSectionHandler.hh"
72#include "G4VEMDataSet.hh"
73#include "G4EMDataSet.hh"
74#include "G4CompositeEMDataSet.hh"
75#include "G4VDataSetAlgorithm.hh"
76#include "G4LogLogInterpolation.hh"
77#include "G4VRangeTest.hh"
78#include "G4RangeTest.hh"
79#include "G4ProductionCutsTable.hh"
80#include "G4AtomicTransitionManager.hh"
81#include "G4AtomicShell.hh"
82#include "G4AtomicDeexcitation.hh"
83#include "G4PenelopeIntegrator.hh"
84#include "G4MaterialCutsCouple.hh"
85
86
87G4PenelopeCompton::G4PenelopeCompton(const G4String& processName)
88 : G4VDiscreteProcess(processName),
89 lowEnergyLimit(250*eV),
90 highEnergyLimit(100*GeV),
91 intrinsicLowEnergyLimit(10*eV),
92 intrinsicHighEnergyLimit(100*GeV),
93 energyForIntegration(0.0),
94 ZForIntegration(1),
95 nBins(200),
96 cutForLowEnergySecondaryPhotons(250.0*eV),
97 fUseAtomicDeexcitation(true)
98{
99 if (lowEnergyLimit < intrinsicLowEnergyLimit ||
100 highEnergyLimit > intrinsicHighEnergyLimit)
101 {
102 G4Exception("G4PenelopeCompton::G4PenelopeCompton - energy outside intrinsic process validity range");
103 }
104
105 meanFreePathTable = 0;
106 ionizationEnergy = new std::map<G4int,G4DataVector*>;
107 hartreeFunction = new std::map<G4int,G4DataVector*>;
108 occupationNumber = new std::map<G4int,G4DataVector*>;
109
110 rangeTest = new G4RangeTest;
111
112 ReadData(); //Read data from file
113
114 if (verboseLevel > 0)
115 {
116 G4cout << GetProcessName() << " is created " << G4endl
117 << "Energy range: "
118 << lowEnergyLimit / keV << " keV - "
119 << highEnergyLimit / GeV << " GeV"
120 << G4endl;
121 }
122
123
124 G4cout << G4endl;
125 G4cout << "*******************************************************************************" << G4endl;
126 G4cout << "*******************************************************************************" << G4endl;
127 G4cout << " The class G4PenelopeCompton is NOT SUPPORTED ANYMORE. " << G4endl;
128 G4cout << " It will be REMOVED with the next major release of Geant4. " << G4endl;
129 G4cout << " Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
130 G4cout << "*******************************************************************************" << G4endl;
131 G4cout << "*******************************************************************************" << G4endl;
132 G4cout << G4endl;
133
134}
135
136G4PenelopeCompton::~G4PenelopeCompton()
137{
138 delete meanFreePathTable;
139 delete rangeTest;
140
141 for (size_t i=0;i<matCrossSections->size();i++)
142 {
143 delete (*matCrossSections)[i];
144 }
145
146 delete matCrossSections;
147
148 for (G4int Z=1;Z<100;Z++)
149 {
150 if (ionizationEnergy->count(Z)) delete (ionizationEnergy->find(Z)->second);
151 if (hartreeFunction->count(Z)) delete (hartreeFunction->find(Z)->second);
152 if (occupationNumber->count(Z)) delete (occupationNumber->find(Z)->second);
153 }
154 delete ionizationEnergy;
155 delete hartreeFunction;
156 delete occupationNumber;
157}
158
159void G4PenelopeCompton::BuildPhysicsTable(const G4ParticleDefinition& )
160{
161 G4DataVector energyVector;
162 G4double dBin = std::log10(highEnergyLimit/lowEnergyLimit)/nBins;
163 for (G4int i=0;i<nBins;i++)
164 {
165 energyVector.push_back(std::pow(10.,std::log10(lowEnergyLimit)+i*dBin));
166 }
167
168 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
169 G4int nMaterials = G4Material::GetNumberOfMaterials();
170 G4VDataSetAlgorithm* algo = new G4LogLogInterpolation();
171
172 //size_t nOfBins = energyVector.size();
173 //size_t bin=0;
174
175 G4DataVector* energies;
176 G4DataVector* data;
177
178 matCrossSections = new std::vector<G4VEMDataSet*>;
179
180 for (G4int m=0; m<nMaterials; m++)
181 {
182 const G4Material* material= (*materialTable)[m];
183 G4int nElements = material->GetNumberOfElements();
184 const G4ElementVector* elementVector = material->GetElementVector();
185 const G4double* nAtomsPerVolume = material->GetAtomicNumDensityVector();
186
187 G4VEMDataSet* setForMat = new G4CompositeEMDataSet(algo,1.,1.);
188
189 for (G4int i=0; i<nElements; i++) {
190
191 G4int Z = (G4int) (*elementVector)[i]->GetZ();
192 G4double density = nAtomsPerVolume[i];
193 G4double cross=0.0;
194 energies = new G4DataVector;
195 data = new G4DataVector;
196
197
198 for (size_t bin=0; bin<energyVector.size(); bin++)
199 {
200 G4double e = energyVector[bin];
201 energies->push_back(e);
202 cross = density * CrossSection(e,Z);
203 data->push_back(cross);
204 }
205
206 G4VEMDataSet* elSet = new G4EMDataSet(i,energies,data,algo,1.,1.);
207 setForMat->AddComponent(elSet);
208 }
209
210 matCrossSections->push_back(setForMat);
211 }
212
213
214 //Build the mean free path table!
215 G4double matCS = 0.0;
216 G4VEMDataSet* matCrossSet = new G4CompositeEMDataSet(algo,1.,1.);
217 G4VEMDataSet* materialSet = new G4CompositeEMDataSet(algo,1.,1.);
218
219
220 for (G4int m=0; m<nMaterials; m++)
221 {
222 energies = new G4DataVector;
223 data = new G4DataVector;
224 const G4Material* material= (*materialTable)[m];
225 material= (*materialTable)[m];
226 for (size_t bin=0; bin<energyVector.size(); bin++)
227 {
228 G4double energy = energyVector[bin];
229 energies->push_back(energy);
230 matCrossSet = (*matCrossSections)[m];
231 matCS = 0.0;
232 G4int nElm = matCrossSet->NumberOfComponents();
233 for(G4int j=0; j<nElm; j++) {
234 matCS += matCrossSet->GetComponent(j)->FindValue(energy);
235 }
236 if (matCS > 0.)
237 {
238 data->push_back(1./matCS);
239 }
240 else
241 {
242 data->push_back(DBL_MAX);
243 }
244 }
245 G4VEMDataSet* dataSet = new G4EMDataSet(m,energies,data,algo,1.,1.);
246 materialSet->AddComponent(dataSet);
247 }
248 meanFreePathTable = materialSet;
249}
250
251G4VParticleChange* G4PenelopeCompton::PostStepDoIt(const G4Track& aTrack,
252 const G4Step& aStep)
253{
254 //Penelope model
255
256 aParticleChange.Initialize(aTrack);
257
258 // Dynamic particle quantities
259 const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
260 G4double photonEnergy0 = incidentPhoton->GetKineticEnergy();
261
262 if (photonEnergy0 <= lowEnergyLimit)
263 {
264 aParticleChange.ProposeTrackStatus(fStopAndKill);
265 aParticleChange.ProposeEnergy(0.);
266 aParticleChange.ProposeLocalEnergyDeposit(photonEnergy0);
267 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
268 }
269
270 G4ParticleMomentum photonDirection0 = incidentPhoton->GetMomentumDirection();
271
272 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
273 const G4Material* material = couple->GetMaterial();
274
275 G4int Z = SelectRandomAtomForCompton(material,photonEnergy0);
276 const G4int nmax = 64;
277 G4double rn[nmax],pac[nmax];
278
279 G4double ki,ki1,ki2,ki3,taumin,a1,a2;
280 G4double tau,TST;
281 G4double S=0.0;
282 G4double epsilon,cosTheta;
283 G4double harFunc = 0.0;
284 G4int occupNb= 0;
285 G4double ionEnergy=0.0;
286 G4int nosc = occupationNumber->find(Z)->second->size();
287 G4int iosc = nosc;
288 ki = photonEnergy0/electron_mass_c2;
289 ki2 = 2*ki+1.0;
290 ki3 = ki*ki;
291 ki1 = ki3-ki2-1.0;
292 taumin = 1.0/ki2;
293 a1 = std::log(ki2);
294 a2 = a1+2.0*ki*(1.0+ki)/(ki2*ki2);
295 //If the incoming photon is above 5 MeV, the quicker approach based on the
296 //pure Klein-Nishina formula is used
297 if (photonEnergy0 > 5*MeV)
298 {
299 do{
300 do{
301 if ((a2*G4UniformRand()) < a1)
302 {
303 tau = std::pow(taumin,G4UniformRand());
304 }
305 else
306 {
307 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
308 }
309 //rejection function
310 TST = (1+tau*(ki1+tau*(ki2+tau*ki3)))/(ki3*tau*(1.0+tau*tau));
311 }while (G4UniformRand()> TST);
312 epsilon=tau;
313 cosTheta = 1.0 - (1.0-tau)/(ki*tau);
314 //Target shell electrons
315 TST = Z*G4UniformRand();
316 iosc = nosc;
317 S=0.0;
318 for (G4int j=0;j<nosc;j++)
319 {
320 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[j];
321 S = S + occupNb;
322 if (S > TST) iosc = j;
323 if (S > TST) break;
324 }
325 ionEnergy = (*(ionizationEnergy->find(Z)->second))[iosc];
326 }while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0);
327 }
328 else //photonEnergy0<5 MeV
329 {
330 //Incoherent scattering function for theta=PI
331 G4double s0=0.0;
332 G4double pzomc=0.0,rni=0.0;
333 G4double aux=0.0;
334 for (G4int i=0;i<nosc;i++){
335 ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
336 if (photonEnergy0 > ionEnergy)
337 {
338 G4double aux = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
339 harFunc = (*(hartreeFunction->find(Z)->second))[i]/fine_structure_const;
340 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
341 pzomc = harFunc*(aux-electron_mass_c2*ionEnergy)/
342 (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
343 if (pzomc > 0)
344 {
345 rni = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*(std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
346 }
347 else
348 {
349 rni = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*(std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
350 }
351 s0 = s0 + occupNb*rni;
352 }
353 }
354
355 //Sampling tau
356 G4double cdt1;
357 do
358 {
359 if ((G4UniformRand()*a2) < a1)
360 {
361 tau = std::pow(taumin,G4UniformRand());
362 }
363 else
364 {
365 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
366 }
367 cdt1 = (1.0-tau)/(ki*tau);
368 S=0.0;
369 //Incoherent scattering function
370 for (G4int i=0;i<nosc;i++){
371 ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
372 if (photonEnergy0 > ionEnergy) //sum only on excitable levels
373 {
374 aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1;
375 harFunc = (*(hartreeFunction->find(Z)->second))[i]/fine_structure_const;
376 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
377 pzomc = harFunc*(aux-electron_mass_c2*ionEnergy)/
378 (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
379 if (pzomc > 0)
380 {
381 rn[i] = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
382 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
383 }
384 else
385 {
386 rn[i] = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
387 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
388 }
389 S = S + occupNb*rn[i];
390 pac[i] = S;
391 }
392 else
393 {
394 pac[i] = S-(1e-06);
395 }
396 }
397 //Rejection function
398 TST = S*(1.0+tau*(ki1+tau*(ki2+tau*ki3)))/(ki3*tau*(1.0+tau*tau));
399 }while ((G4UniformRand()*s0) > TST);
400
401 //Target electron shell
402 cosTheta = 1.0 - cdt1;
403 G4double fpzmax=0.0,fpz=0.0;
404 G4double A=0.0;
405 do
406 {
407 do
408 {
409 TST =S*G4UniformRand();
410 iosc=nosc;
411 for (G4int i=0;i<nosc;i++){
412 if (pac[i]>TST) iosc = i;
413 if (pac[i]>TST) break;
414 }
415 A = G4UniformRand()*rn[iosc];
416 harFunc = (*(hartreeFunction->find(Z)->second))[iosc]/fine_structure_const;
417 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[iosc];
418 if (A < 0.5) {
419 pzomc = (std::sqrt(0.5)-std::sqrt(0.5-std::log(2.0*A)))/
420 (std::sqrt(2.0)*harFunc);
421 }
422 else
423 {
424 pzomc = (std::sqrt(0.5-std::log(2.0-2.0*A))-std::sqrt(0.5))/
425 (std::sqrt(2.0)*harFunc);
426 }
427 } while (pzomc < -1);
428 // F(EP) rejection
429 G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
430 G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC);
431 if (AF > 0) {
432 fpzmax = 1.0+AF*0.2;
433 }
434 else
435 {
436 fpzmax = 1.0-AF*0.2;
437 }
438 fpz = 1.0+AF*std::max(std::min(pzomc,0.2),-0.2);
439 }while ((fpzmax*G4UniformRand())>fpz);
440
441 //Energy of the scattered photon
442 G4double T = pzomc*pzomc;
443 G4double b1 = 1.0-T*tau*tau;
444 G4double b2 = 1.0-T*tau*cosTheta;
445 if (pzomc > 0.0)
446 {
447 epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
448 }
449 else
450 {
451 epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
452 }
453 }
454
455 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
456 G4double phi = twopi * G4UniformRand() ;
457 G4double dirx = sinTheta * std::cos(phi);
458 G4double diry = sinTheta * std::sin(phi);
459 G4double dirz = cosTheta ;
460
461 // Update G4VParticleChange for the scattered photon
462
463 G4ThreeVector photonDirection1(dirx,diry,dirz);
464 photonDirection1.rotateUz(photonDirection0);
465 aParticleChange.ProposeMomentumDirection(photonDirection1) ;
466 G4double photonEnergy1 = epsilon * photonEnergy0;
467
468 if (photonEnergy1 > 0.)
469 {
470 aParticleChange.ProposeEnergy(photonEnergy1) ;
471 }
472 else
473 {
474 aParticleChange.ProposeEnergy(0.) ;
475 aParticleChange.ProposeTrackStatus(fStopAndKill);
476 }
477
478
479 // Kinematics of the scattered electron
480 G4double diffEnergy = photonEnergy0*(1-epsilon);
481 ionEnergy = (*(ionizationEnergy->find(Z)->second))[iosc];
482 G4double Q2 = photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
483 G4double cosThetaE; //scattering angle for the electron
484 if (Q2 > 1.0e-12)
485 {
486 cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2);
487 }
488 else
489 {
490 cosThetaE = 1.0;
491 }
492 G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE);
493 //initialize here, then check photons created by Atomic-Deexcitation, and the final state e-
494 G4int nbOfSecondaries = 0;
495
496 std::vector<G4DynamicParticle*>* photonVector=0;
497
498 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
499 const G4AtomicShell* shell = transitionManager->Shell(Z,iosc);
500 G4double bindingEnergy = shell->BindingEnergy();
501 G4int shellId = shell->ShellId();
502 G4double ionEnergyInPenelopeDatabase = ionEnergy;
503 ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase); //protection against energy non-conservation
504
505 G4double eKineticEnergy = diffEnergy - ionEnergy; //subtract the excitation energy. If not emitted by fluorescence,
506 //the ionization energy is deposited as local energy deposition
507 G4double localEnergyDeposit = ionEnergy;
508 G4double energyInFluorescence = 0.; //testing purposes only
509
510 if (eKineticEnergy < 0)
511 {
512 //It means that there was some problem/mismatch between the two databases. Try to make it work
513 //In this case available Energy (diffEnergy) < ionEnergy
514 //Full residual energy is deposited locally
515 localEnergyDeposit = diffEnergy;
516 eKineticEnergy = 0.0;
517 }
518
519 //the local energy deposit is what remains: part of this may be spent for fluorescence.
520
521 if (fUseAtomicDeexcitation)
522 {
523 G4int nPhotons=0;
524
525 const G4ProductionCutsTable* theCoupleTable=
526 G4ProductionCutsTable::GetProductionCutsTable();
527 size_t indx = couple->GetIndex();
528
529 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
530 cutg = std::max(cutForLowEnergySecondaryPhotons,cutg);
531
532 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
533 cute = std::max(cutForLowEnergySecondaryPhotons,cute);
534
535 G4DynamicParticle* aPhoton;
536 G4AtomicDeexcitation deexcitationManager;
537
538 if (Z>5 && (localEnergyDeposit > cutg || localEnergyDeposit > cute))
539 {
540 photonVector = deexcitationManager.GenerateParticles(Z,shellId);
541 for (size_t k=0;k<photonVector->size();k++){
542 aPhoton = (*photonVector)[k];
543 if (aPhoton)
544 {
545 G4double itsCut = cutg;
546 if (aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
547 G4double itsEnergy = aPhoton->GetKineticEnergy();
548 if (itsEnergy > itsCut && itsEnergy <= localEnergyDeposit)
549 {
550 nPhotons++;
551 localEnergyDeposit -= itsEnergy;
552 energyInFluorescence += itsEnergy;
553 }
554 else
555 {
556 delete aPhoton;
557 (*photonVector)[k]=0;
558 }
559 }
560 }
561 }
562 nbOfSecondaries=nPhotons;
563 }
564
565
566 // Generate the electron only if with large enough range w.r.t. cuts and safety
567 G4double safety = aStep.GetPostStepPoint()->GetSafety();
568 G4DynamicParticle* electron = 0;
569 if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety) &&
570 eKineticEnergy>cutForLowEnergySecondaryPhotons)
571 {
572 G4double xEl = sinThetaE * std::cos(phi+pi);
573 G4double yEl = sinThetaE * std::sin(phi+pi);
574 G4double zEl = cosThetaE;
575 G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
576 eDirection.rotateUz(photonDirection0);
577 electron = new G4DynamicParticle (G4Electron::Electron(),
578 eDirection,eKineticEnergy) ;
579 nbOfSecondaries++;
580 }
581 else
582 {
583 localEnergyDeposit += eKineticEnergy;
584 }
585
586 aParticleChange.SetNumberOfSecondaries(nbOfSecondaries);
587 if (electron) aParticleChange.AddSecondary(electron);
588
589 //This block below is executed only if there is at least one secondary photon produced by
590 //AtomicDeexcitation
591 if (photonVector)
592 {
593 for (size_t ll=0;ll<photonVector->size();ll++)
594 {
595 if ((*photonVector)[ll]) aParticleChange.AddSecondary((*photonVector)[ll]);
596 }
597 }
598 delete photonVector;
599 if (localEnergyDeposit < 0)
600 {
601 G4cout << "WARNING-"
602 << "G4PenelopeCompton::PostStepDoIt - Negative energy deposit"
603 << G4endl;
604 localEnergyDeposit=0.;
605 }
606 aParticleChange.ProposeLocalEnergyDeposit(localEnergyDeposit);
607
608
609 if (verboseLevel > 1)
610 {
611 G4cout << "-----------------------------------------------------------" << G4endl;
612 G4cout << "Energy balance from G4PenelopeCompton" << G4endl;
613 G4cout << "Incoming photon energy: " << photonEnergy0/keV << " keV" << G4endl;
614 G4cout << "-----------------------------------------------------------" << G4endl;
615 G4cout << "Scattered photon: " << photonEnergy1/keV << " keV" << G4endl;
616 G4double electronEnergy = 0.;
617 if (electron)
618 electronEnergy = eKineticEnergy;
619 G4cout << "Scattered electron " << electronEnergy/keV << " keV" << G4endl;
620 G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
621 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
622 G4cout << "Total final state: " << (photonEnergy1+electronEnergy+energyInFluorescence+localEnergyDeposit)/keV <<
623 " keV" << G4endl;
624 G4cout << "-----------------------------------------------------------" << G4endl;
625 }
626
627 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
628}
629
630G4bool G4PenelopeCompton::IsApplicable(const G4ParticleDefinition& particle)
631{
632 return ( &particle == G4Gamma::Gamma() );
633}
634
635G4double G4PenelopeCompton::GetMeanFreePath(const G4Track& track,
636 G4double, // previousStepSize
637 G4ForceCondition*)
638{
639 const G4DynamicParticle* photon = track.GetDynamicParticle();
640 G4double energy = photon->GetKineticEnergy();
641 G4Material* material = track.GetMaterial();
642 size_t materialIndex = material->GetIndex();
643
644 G4double meanFreePath;
645 if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
646 else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
647 else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
648 return meanFreePath;
649}
650
651
652void G4PenelopeCompton::ReadData()
653{
654 char* path = getenv("G4LEDATA");
655 if (!path)
656 {
657 G4String excep = "G4PenelopeCompton - G4LEDATA environment variable not set!";
658 G4Exception(excep);
659 }
660 G4String pathString(path);
661 G4String pathFile = pathString + "/penelope/compton-pen.dat";
662 std::ifstream file(pathFile);
663 std::filebuf* lsdp = file.rdbuf();
664
665 if (!(lsdp->is_open()))
666 {
667 G4String excep = "G4PenelopeCompton - data file " + pathFile + " not found!";
668 G4Exception(excep);
669 }
670
671 G4int k1,test,test1;
672 G4double a1,a2;
673 G4int Z=1,nLevels=0;
674 G4DataVector* f;
675 G4DataVector* u;
676 G4DataVector* j;
677
678 do{
679 f = new G4DataVector;
680 u = new G4DataVector;
681 j = new G4DataVector;
682 file >> Z >> nLevels;
683 for (G4int h=0;h<nLevels;h++){
684 file >> k1 >> a1 >> a2;
685 f->push_back((G4double) k1);
686 u->push_back(a1);
687 j->push_back(a2);
688 }
689 ionizationEnergy->insert(std::make_pair(Z,u));
690 hartreeFunction->insert(std::make_pair(Z,j));
691 occupationNumber->insert(std::make_pair(Z,f));
692 file >> test >> test1; //-1 -1 close the data for each Z
693 if (test > 0) {
694 G4String excep = "G4PenelopeCompton - data file corrupted!";
695 G4Exception(excep);
696 }
697 }while (test != -2); //the very last Z is closed with -2 instead of -1
698}
699
700G4double G4PenelopeCompton::CrossSection(G4double energy,G4int Z)
701{
702 G4double cs=0.0;
703 energyForIntegration=energy;
704 ZForIntegration = Z;
705 if (energy< 5*MeV)
706 {
707 G4PenelopeIntegrator<G4PenelopeCompton,G4double (G4PenelopeCompton::*)(G4double)> theIntegrator;
708 cs = theIntegrator.Calculate(this,&G4PenelopeCompton::DifferentialCrossSection,-1.0,1.0,1e-05);
709 }
710 else
711 {
712 G4double ki=energy/electron_mass_c2;
713 G4double ki3=ki*ki;
714 G4double ki2=1.0+2*ki;
715 G4double ki1=ki3-ki2-1.0;
716 G4double t0=1.0/(ki2);
717 G4double csl = 0.5*ki3*t0*t0+ki2*t0+ki1*std::log(t0)-(1.0/t0);
718 G4int nosc = occupationNumber->find(Z)->second->size();
719 for (G4int i=0;i<nosc;i++)
720 {
721 G4double ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
722 G4double tau=(energy-ionEnergy)/energy;
723 if (tau > t0)
724 {
725 G4double csu = 0.5*ki3*tau*tau+ki2*tau+ki1*std::log(tau)-(1.0/tau);
726 G4int f = (G4int) (*(occupationNumber->find(Z)->second))[i];
727 cs = cs + f*(csu-csl);
728 }
729 }
730 cs=pi*classic_electr_radius*classic_electr_radius*cs/(ki*ki3);
731 }
732 return cs;
733}
734
735
736G4double G4PenelopeCompton::DifferentialCrossSection(G4double cosTheta)
737{
738 const G4double k2 = std::sqrt(2.0);
739 const G4double k1 = std::sqrt(0.5);
740 const G4double k12 = 0.5;
741 G4double cdt1 = 1.0-cosTheta;
742 G4double energy = energyForIntegration;
743 G4int Z = ZForIntegration;
744 G4double ionEnergy=0.0,Pzimax=0.0,XKN=0.0;
745 G4double diffCS=0.0;
746 G4double x=0.0,siap=0.0;
747 G4double harFunc=0.0;
748 G4int occupNb;
749 //energy of Compton line;
750 G4double EOEC = 1.0+(energy/electron_mass_c2)*cdt1;
751 G4double ECOE = 1.0/EOEC;
752 //Incoherent scattering function (analytical profile)
753 G4double sia = 0.0;
754 G4int nosc = occupationNumber->find(Z)->second->size();
755 for (G4int i=0;i<nosc;i++){
756 ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
757 //Sum only of those shells for which E>Eion
758 if (energy > ionEnergy)
759 {
760 G4double aux = energy * (energy-ionEnergy)*cdt1;
761 Pzimax = (aux - electron_mass_c2*ionEnergy)/(electron_mass_c2*std::sqrt(2*aux+ionEnergy*ionEnergy));
762 harFunc = (*(hartreeFunction->find(Z)->second))[i]/fine_structure_const;
763 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
764 x = harFunc*Pzimax;
765 if (x > 0)
766 {
767 siap = 1.0-0.5*std::exp(k12-(k1+k2*x)*(k1+k2*x));
768 }
769 else
770 {
771 siap = 0.5*std::exp(k12-(k1-k2*x)*(k1-k2*x));
772 }
773 sia = sia + occupNb*siap; //sum of all contributions;
774 }
775 }
776 XKN = EOEC+ECOE-1+cosTheta*cosTheta;
777 diffCS = pi*classic_electr_radius*classic_electr_radius*ECOE*ECOE*XKN*sia;
778 return diffCS;
779}
780
781G4int G4PenelopeCompton::SelectRandomAtomForCompton(const G4Material* material,G4double energy) const
782{
783 G4int nElements = material->GetNumberOfElements();
784 //Special case: the material consists of one element
785 if (nElements == 1)
786 {
787 G4int Z = (G4int) material->GetZ();
788 return Z;
789 }
790
791 //Composite material
792 const G4ElementVector* elementVector = material->GetElementVector();
793 size_t materialIndex = material->GetIndex();
794
795 G4VEMDataSet* materialSet = (*matCrossSections)[materialIndex];
796 G4double materialCrossSection0 = 0.0;
797 G4DataVector cross;
798 cross.clear();
799 G4int i;
800 for (i=0;i<nElements;i++)
801 {
802 G4double cr = (materialSet->GetComponent(i))->FindValue(energy);
803 materialCrossSection0 += cr;
804 cross.push_back(materialCrossSection0); //cumulative cross section
805 }
806
807 G4double random = G4UniformRand()*materialCrossSection0;
808 for (i=0;i<nElements;i++)
809 {
810 if (random <= cross[i]) return (G4int) (*elementVector)[i]->GetZ();
811 }
812 //It should never get here
813 return 0;
814}
815
Note: See TracBrowser for help on using the repository browser.