source: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreIonisationModel.cc@ 1190

Last change on this file since 1190 was 1058, checked in by garnier, 17 years ago

file release beta
File size: 22.8 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: G4LivermoreIonisationModel.cc,v 1.4 2009/05/19 14:57:01 pandola Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 12 Jan 2009 L Pandola Migration from process to model
34// 03 Mar 2009 L Pandola Bug fix (release memory in the destructor)
35// 15 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
36// - apply internal high-energy limit only in constructor
37// - do not apply low-energy limit (default is 0)
38// - simplify sampling of deexcitation by using cut in energy
39// - set activation of Auger "false"
40// - remove initialisation of element selectors
41// 19 May 2009 L Pandola Explicitely set to zero pointers deleted in
42// Initialise(), since they might be checked later on
43//
44
45#include "G4LivermoreIonisationModel.hh"
46#include "G4ParticleDefinition.hh"
47#include "G4MaterialCutsCouple.hh"
48#include "G4ProductionCutsTable.hh"
49#include "G4DynamicParticle.hh"
50#include "G4Element.hh"
51#include "G4AtomicTransitionManager.hh"
52#include "G4AtomicDeexcitation.hh"
53#include "G4AtomicShell.hh"
54#include "G4Gamma.hh"
55#include "G4Electron.hh"
56#include "G4CrossSectionHandler.hh"
57#include "G4AtomicDeexcitation.hh"
58#include "G4ProcessManager.hh"
59#include "G4VEMDataSet.hh"
60#include "G4EMDataSet.hh"
61#include "G4eIonisationCrossSectionHandler.hh"
62#include "G4eIonisationSpectrum.hh"
63#include "G4VDataSetAlgorithm.hh"
64#include "G4SemiLogInterpolation.hh"
65#include "G4ShellVacancy.hh"
66#include "G4VDataSetAlgorithm.hh"
67#include "G4LogLogInterpolation.hh"
68#include "G4CompositeEMDataSet.hh"
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72
73G4LivermoreIonisationModel::G4LivermoreIonisationModel(const G4ParticleDefinition*,
74 const G4String& nam)
75 :G4VEmModel(nam),isInitialised(false),crossSectionHandler(0),
76 energySpectrum(0),shellVacancy(0)
77{
78 fIntrinsicLowEnergyLimit = 10.0*eV;
79 fIntrinsicHighEnergyLimit = 100.0*GeV;
80 fNBinEnergyLoss = 360;
81 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
82 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
83 //
84 verboseLevel = 0;
85 //
86 fUseAtomicDeexcitation = true;
87 ActivateAuger(false);
88 //
89 // Notice: the fluorescence along step is generated only if it is
90 // set by the PROCESS (e.g. G4eIonisation) via the command
91 // process->ActivateDeexcitation(true);
92 //
93}
94
95//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96
97G4LivermoreIonisationModel::~G4LivermoreIonisationModel()
98{
99 if (energySpectrum) delete energySpectrum;
100 if (crossSectionHandler) delete crossSectionHandler;
101 if (shellVacancy) delete shellVacancy;
102}
103
104//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105
106void G4LivermoreIonisationModel::Initialise(const G4ParticleDefinition* particle,
107 const G4DataVector& cuts)
108{
109 //Check that the Livermore Ionisation is NOT attached to e+
110 if (particle != G4Electron::Electron())
111 {
112 G4cout << "ERROR: Livermore Ionisation Model is applicable only to electrons" << G4endl;
113 G4cout << "It cannot be registered to " << particle->GetParticleName() << G4endl;
114 G4Exception();
115 }
116
117 //Read energy spectrum
118 if (energySpectrum)
119 {
120 delete energySpectrum;
121 energySpectrum = 0;
122 }
123 energySpectrum = new G4eIonisationSpectrum();
124 if (verboseLevel > 0)
125 G4cout << "G4VEnergySpectrum is initialized" << G4endl;
126
127 //Initialize cross section handler
128 if (crossSectionHandler)
129 {
130 delete crossSectionHandler;
131 crossSectionHandler = 0;
132 }
133
134 G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation();
135 crossSectionHandler = new G4eIonisationCrossSectionHandler(energySpectrum,interpolation,
136 LowEnergyLimit(),HighEnergyLimit(),
137 fNBinEnergyLoss);
138 crossSectionHandler->Clear();
139 crossSectionHandler->LoadShellData("ioni/ion-ss-cs-");
140 //This is used to retrieve cross section values later on
141 crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
142
143 //Fluorescence data
144 transitionManager = G4AtomicTransitionManager::Instance();
145 if (shellVacancy) delete shellVacancy;
146 shellVacancy = new G4ShellVacancy();
147 InitialiseFluorescence();
148
149 if (verboseLevel > 0)
150 {
151 G4cout << "Livermore Ionisation model is initialized " << G4endl
152 << "Energy range: "
153 << LowEnergyLimit() / keV << " keV - "
154 << HighEnergyLimit() / GeV << " GeV"
155 << G4endl;
156 }
157
158 if (verboseLevel > 1)
159 {
160 G4cout << "Cross section data: " << G4endl;
161 crossSectionHandler->PrintData();
162 G4cout << "Parameters: " << G4endl;
163 energySpectrum->PrintData();
164 }
165
166 if(isInitialised) return;
167 fParticleChange = GetParticleChangeForLoss();
168 isInitialised = true;
169}
170
171//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172
173G4double G4LivermoreIonisationModel::MinEnergyCut(const G4ParticleDefinition*,
174 const G4MaterialCutsCouple*)
175{
176 return 250.*eV;
177}
178
179//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180
181G4double G4LivermoreIonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
182 G4double energy,
183 G4double Z, G4double,
184 G4double cutEnergy,
185 G4double)
186{
187 G4int iZ = (G4int) Z;
188 if (!crossSectionHandler)
189 {
190 G4cout << "G4LivermoreIonisationModel::ComputeCrossSectionPerAtom" << G4endl;
191 G4cout << "The cross section handler is not correctly initialized" << G4endl;
192 G4Exception();
193 }
194
195 //The cut is already included in the crossSectionHandler
196 G4double cs =
197 crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,cutEnergy,iZ);
198
199 if (verboseLevel > 1)
200 {
201 G4cout << "G4LivermoreIonisationModel " << G4endl;
202 G4cout << "Cross section for delta emission > " << cutEnergy/keV << " keV at " <<
203 energy/keV << " keV and Z = " << iZ << " --> " << cs/barn << " barn" << G4endl;
204 }
205 return cs;
206}
207
208
209//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
210
211G4double G4LivermoreIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
212 const G4ParticleDefinition* ,
213 G4double kineticEnergy,
214 G4double cutEnergy)
215{
216 G4double sPower = 0.0;
217
218 const G4ElementVector* theElementVector = material->GetElementVector();
219 size_t NumberOfElements = material->GetNumberOfElements() ;
220 const G4double* theAtomicNumDensityVector =
221 material->GetAtomicNumDensityVector();
222
223 // loop for elements in the material
224 for (size_t iel=0; iel<NumberOfElements; iel++ )
225 {
226 G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
227 G4int nShells = transitionManager->NumberOfShells(iZ);
228 for (G4int n=0; n<nShells; n++)
229 {
230 G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
231 kineticEnergy, n);
232 G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n);
233 sPower += e * cs * theAtomicNumDensityVector[iel];
234 }
235 G4double esp = energySpectrum->Excitation(iZ,kineticEnergy);
236 sPower += esp * theAtomicNumDensityVector[iel];
237 }
238
239 if (verboseLevel > 2)
240 {
241 G4cout << "G4LivermoreIonisationModel " << G4endl;
242 G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
243 kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl;
244 }
245
246 return sPower;
247}
248
249//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
250
251void G4LivermoreIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
252 const G4MaterialCutsCouple* couple,
253 const G4DynamicParticle* aDynamicParticle,
254 G4double cutE,
255 G4double maxE)
256{
257
258 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
259
260 if (kineticEnergy <= fIntrinsicLowEnergyLimit)
261 {
262 fParticleChange->SetProposedKineticEnergy(0.);
263 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
264 return ;
265 }
266
267 // Select atom and shell
268 G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
269 G4int shell = crossSectionHandler->SelectRandomShell(Z, kineticEnergy);
270 const G4AtomicShell* atomicShell =
271 (G4AtomicTransitionManager::Instance())->Shell(Z, shell);
272 G4double bindingEnergy = atomicShell->BindingEnergy();
273 G4int shellId = atomicShell->ShellId();
274
275 // Sample delta energy using energy interval for delta-electrons
276 G4double energyMax =
277 std::min(maxE,energySpectrum->MaxEnergyOfSecondaries(kineticEnergy));
278 G4double energyDelta = energySpectrum->SampleEnergy(Z, cutE, energyMax,
279 kineticEnergy, shell);
280
281 if (energyDelta == 0.) //nothing happens
282 return;
283
284 // Transform to shell potential
285 G4double deltaKinE = energyDelta + 2.0*bindingEnergy;
286 G4double primaryKinE = kineticEnergy + 2.0*bindingEnergy;
287
288 // sampling of scattering angle neglecting atomic motion
289 G4double deltaMom = std::sqrt(deltaKinE*(deltaKinE + 2.0*electron_mass_c2));
290 G4double primaryMom = std::sqrt(primaryKinE*(primaryKinE + 2.0*electron_mass_c2));
291
292 G4double cost = deltaKinE * (primaryKinE + 2.0*electron_mass_c2)
293 / (deltaMom * primaryMom);
294 if (cost > 1.) cost = 1.;
295 G4double sint = std::sqrt(1. - cost*cost);
296 G4double phi = twopi * G4UniformRand();
297 G4double dirx = sint * std::cos(phi);
298 G4double diry = sint * std::sin(phi);
299 G4double dirz = cost;
300
301 // Rotate to incident electron direction
302 G4ThreeVector primaryDirection = aDynamicParticle->GetMomentumDirection();
303 G4ThreeVector deltaDir(dirx,diry,dirz);
304 deltaDir.rotateUz(primaryDirection);
305 //Updated components
306 dirx = deltaDir.x();
307 diry = deltaDir.y();
308 dirz = deltaDir.z();
309
310 // Take into account atomic motion del is relative momentum of the motion
311 // kinetic energy of the motion == bindingEnergy in V.Ivanchenko model
312 cost = 2.0*G4UniformRand() - 1.0;
313 sint = std::sqrt(1. - cost*cost);
314 phi = twopi * G4UniformRand();
315 G4double del = std::sqrt(bindingEnergy *(bindingEnergy + 2.0*electron_mass_c2))
316 / deltaMom;
317 dirx += del* sint * std::cos(phi);
318 diry += del* sint * std::sin(phi);
319 dirz += del* cost;
320
321 // Find out new primary electron direction
322 G4double finalPx = primaryMom*primaryDirection.x() - deltaMom*dirx;
323 G4double finalPy = primaryMom*primaryDirection.y() - deltaMom*diry;
324 G4double finalPz = primaryMom*primaryDirection.z() - deltaMom*dirz;
325
326 //Ok, ready to create the delta ray
327 G4DynamicParticle* theDeltaRay = new G4DynamicParticle();
328 theDeltaRay->SetKineticEnergy(energyDelta);
329 G4double norm = 1.0/std::sqrt(dirx*dirx + diry*diry + dirz*dirz);
330 dirx *= norm;
331 diry *= norm;
332 dirz *= norm;
333 theDeltaRay->SetMomentumDirection(dirx, diry, dirz);
334 theDeltaRay->SetDefinition(G4Electron::Electron());
335 fvect->push_back(theDeltaRay);
336
337 //This is the amount of energy available for fluorescence
338 G4double theEnergyDeposit = bindingEnergy;
339
340 // fill ParticleChange
341 // changed energy and momentum of the actual particle
342 G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit;
343 if(finalKinEnergy < 0.0)
344 {
345 theEnergyDeposit += finalKinEnergy;
346 finalKinEnergy = 0.0;
347 }
348 else
349 {
350 G4double norm = 1.0/std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
351 finalPx *= norm;
352 finalPy *= norm;
353 finalPz *= norm;
354 fParticleChange->ProposeMomentumDirection(finalPx, finalPy, finalPz);
355 }
356 fParticleChange->SetProposedKineticEnergy(finalKinEnergy);
357
358 // deexcitation may be active per G4Region
359 if(DeexcitationFlag() && Z > 5)
360 {
361 G4ProductionCutsTable* theCoupleTable =
362 G4ProductionCutsTable::GetProductionCutsTable();
363 // Retrieve cuts for gammas
364 G4double cutG = (*theCoupleTable->GetEnergyCutsVector(0))[couple->GetIndex()];
365 if(theEnergyDeposit > cutG || theEnergyDeposit > cutE)
366 {
367 deexcitationManager.SetCutForSecondaryPhotons(cutG);
368 deexcitationManager.SetCutForAugerElectrons(cutE);
369 std::vector<G4DynamicParticle*>* secondaryVector =
370 deexcitationManager.GenerateParticles(Z, shellId);
371 G4DynamicParticle* aSecondary;
372
373 if (secondaryVector)
374 {
375 for (size_t i = 0; i<secondaryVector->size(); i++)
376 {
377 aSecondary = (*secondaryVector)[i];
378 //Check if it is a valid secondary
379 if (aSecondary)
380 {
381 G4double e = aSecondary->GetKineticEnergy();
382 if (e < theEnergyDeposit)
383 {
384 theEnergyDeposit -= e;
385 fvect->push_back(aSecondary);
386 }
387 else
388 {
389 delete aSecondary;
390 (*secondaryVector)[i] = 0;
391 }
392 }
393 }
394 }
395 }
396 }
397
398 if (theEnergyDeposit < 0)
399 {
400 G4cout << "G4LivermoreIonisationModel: Negative energy deposit: "
401 << theEnergyDeposit/eV << " eV" << G4endl;
402 theEnergyDeposit = 0.0;
403 }
404
405 //Assign local energy deposit
406 fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit);
407
408 if (verboseLevel > 1)
409 {
410 G4cout << "-----------------------------------------------------------" << G4endl;
411 G4cout << "Energy balance from G4LivermoreIonisation" << G4endl;
412 G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
413 G4cout << "-----------------------------------------------------------" << G4endl;
414 G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl;
415 G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl;
416 G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl;
417 G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl;
418 G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy+
419 theEnergyDeposit)/keV << " keV" << G4endl;
420 G4cout << "-----------------------------------------------------------" << G4endl;
421 }
422 return;
423}
424
425//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
426
427
428void G4LivermoreIonisationModel::SampleDeexcitationAlongStep(const G4Material* theMaterial,
429 const G4Track& theTrack,
430 G4double& eloss)
431{
432 //No call if there is no deexcitation along step
433 if (!DeexcitationFlag()) return;
434
435 //This method gets the energy loss calculated "Along the Step" and
436 //(including fluctuations) and produces explicit fluorescence/Auger
437 //secondaries. The eloss value is updated.
438 G4double energyLossBefore = eloss;
439 if (verboseLevel > 2)
440 G4cout << "Energy loss along step before deexcitation : " << energyLossBefore/keV <<
441 " keV" << G4endl;
442
443 G4double incidentEnergy = theTrack.GetDynamicParticle()->GetKineticEnergy();
444
445 G4ProductionCutsTable* theCoupleTable =
446 G4ProductionCutsTable::GetProductionCutsTable();
447 const G4MaterialCutsCouple* couple = theTrack.GetMaterialCutsCouple();
448 size_t index = couple->GetIndex();
449 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
450 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
451
452 //Notice: in LowEnergyIonisation, fluorescence is always generated above 250 eV
453 //not above the tracking cut.
454 //G4double cutForLowEnergySecondaryParticles = 250.0*eV;
455
456 std::vector<G4DynamicParticle*>* deexcitationProducts =
457 new std::vector<G4DynamicParticle*>;
458
459 if(eloss > cute || eloss > cutg)
460 {
461 const G4AtomicTransitionManager* transitionManager =
462 G4AtomicTransitionManager::Instance();
463 deexcitationManager.SetCutForSecondaryPhotons(cutg);
464 deexcitationManager.SetCutForAugerElectrons(cute);
465
466 size_t nElements = theMaterial->GetNumberOfElements();
467 const G4ElementVector* theElementVector = theMaterial->GetElementVector();
468
469 std::vector<G4DynamicParticle*>* secVector = 0;
470 G4DynamicParticle* aSecondary = 0;
471 //G4ParticleDefinition* type = 0;
472 G4double e;
473 G4ThreeVector position;
474 G4int shell, shellId;
475
476 // sample secondaries
477 G4double eTot = 0.0;
478 std::vector<G4int> n =
479 shellVacancy->GenerateNumberOfIonisations(couple,
480 incidentEnergy,eloss);
481 for (size_t i=0; i<nElements; i++)
482 {
483 G4int Z = (G4int)((*theElementVector)[i]->GetZ());
484 size_t nVacancies = n[i];
485 G4double maxE = transitionManager->Shell(Z, 0)->BindingEnergy();
486 if (nVacancies > 0 && Z > 5 && (maxE > cute || maxE > cutg))
487 {
488 for (size_t j=0; j<nVacancies; j++)
489 {
490 shell = crossSectionHandler->SelectRandomShell(Z, incidentEnergy);
491 shellId = transitionManager->Shell(Z, shell)->ShellId();
492 G4double maxEShell =
493 transitionManager->Shell(Z, shell)->BindingEnergy();
494 if (maxEShell > cute || maxEShell > cutg )
495 {
496 secVector = deexcitationManager.GenerateParticles(Z, shellId);
497 if (secVector)
498 {
499 for (size_t l = 0; l<secVector->size(); l++) {
500 aSecondary = (*secVector)[l];
501 if (aSecondary)
502 {
503 e = aSecondary->GetKineticEnergy();
504 if ( eTot + e <= eloss )
505 {
506 eTot += e;
507 deexcitationProducts->push_back(aSecondary);
508 }
509 else
510 {
511 delete aSecondary;
512 }
513 }
514 }
515 delete secVector;
516 }
517 }
518 }
519 }
520 }
521 }
522
523 size_t nSecondaries = deexcitationProducts->size();
524 if (nSecondaries > 0)
525 {
526 fParticleChange->SetNumberOfSecondaries(nSecondaries);
527 const G4StepPoint* preStep = theTrack.GetStep()->GetPreStepPoint();
528 const G4StepPoint* postStep = theTrack.GetStep()->GetPostStepPoint();
529 G4ThreeVector r = preStep->GetPosition();
530 G4ThreeVector deltaR = postStep->GetPosition();
531 deltaR -= r;
532 G4double t = preStep->GetGlobalTime();
533 G4double deltaT = postStep->GetGlobalTime();
534 deltaT -= t;
535 G4double time, q;
536 G4ThreeVector position;
537
538 for (size_t i=0; i<nSecondaries; i++)
539 {
540 G4DynamicParticle* part = (*deexcitationProducts)[i];
541 if (part)
542 {
543 G4double eSecondary = part->GetKineticEnergy();
544 eloss -= eSecondary;
545 if (eloss > 0.)
546 {
547 q = G4UniformRand();
548 time = deltaT*q + t;
549 position = deltaR*q;
550 position += r;
551 G4Track* newTrack = new G4Track(part, time, position);
552 pParticleChange->AddSecondary(newTrack);
553 }
554 else
555 {
556 eloss += eSecondary;
557 delete part;
558 part = 0;
559 }
560 }
561 }
562 }
563 delete deexcitationProducts;
564
565 if (verboseLevel > 2)
566 G4cout << "Energy loss along step after deexcitation : " << eloss/keV <<
567 " keV" << G4endl;
568}
569
570//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
571
572void G4LivermoreIonisationModel::InitialiseFluorescence()
573{
574 G4DataVector* ksi = 0;
575 G4DataVector* energyVector = 0;
576 size_t binForFluo = fNBinEnergyLoss/10;
577
578 G4PhysicsLogVector* eVector = new G4PhysicsLogVector(LowEnergyLimit(),HighEnergyLimit(),
579 binForFluo);
580 const G4ProductionCutsTable* theCoupleTable=
581 G4ProductionCutsTable::GetProductionCutsTable();
582 size_t numOfCouples = theCoupleTable->GetTableSize();
583
584 // Loop on couples
585 for (size_t m=0; m<numOfCouples; m++)
586 {
587 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(m);
588 const G4Material* material= couple->GetMaterial();
589
590 const G4ElementVector* theElementVector = material->GetElementVector();
591 size_t NumberOfElements = material->GetNumberOfElements() ;
592 const G4double* theAtomicNumDensityVector =
593 material->GetAtomicNumDensityVector();
594
595 G4VDataSetAlgorithm* interp = new G4LogLogInterpolation();
596 G4VEMDataSet* xsis = new G4CompositeEMDataSet(interp, 1., 1.);
597 //loop on elements
598 G4double energyCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[m];
599 for (size_t iel=0; iel<NumberOfElements; iel++ )
600 {
601 G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
602 energyVector = new G4DataVector();
603 ksi = new G4DataVector();
604 //Loop on energy
605 for (size_t j = 0; j<binForFluo; j++)
606 {
607 G4double energy = eVector->GetLowEdgeEnergy(j);
608 G4double cross = 0.;
609 G4double eAverage= 0.;
610 G4int nShells = transitionManager->NumberOfShells(iZ);
611
612 for (G4int n=0; n<nShells; n++)
613 {
614 G4double e = energySpectrum->AverageEnergy(iZ, 0.0,energyCut,
615 energy, n);
616 G4double pro = energySpectrum->Probability(iZ, 0.0,energyCut,
617 energy, n);
618 G4double cs= crossSectionHandler->FindValue(iZ, energy, n);
619 eAverage += e * cs * theAtomicNumDensityVector[iel];
620 cross += cs * pro * theAtomicNumDensityVector[iel];
621 if(verboseLevel > 1)
622 {
623 G4cout << "Z= " << iZ
624 << " shell= " << n
625 << " E(keV)= " << energy/keV
626 << " Eav(keV)= " << e/keV
627 << " pro= " << pro
628 << " cs= " << cs
629 << G4endl;
630 }
631 }
632
633 G4double coeff = 0.0;
634 if(eAverage > 0.)
635 {
636 coeff = cross/eAverage;
637 eAverage /= cross;
638 }
639
640 if(verboseLevel > 1)
641 {
642 G4cout << "Ksi Coefficient for Z= " << iZ
643 << " E(keV)= " << energy/keV
644 << " Eav(keV)= " << eAverage/keV
645 << " coeff= " << coeff
646 << G4endl;
647 }
648 energyVector->push_back(energy);
649 ksi->push_back(coeff);
650 }
651 G4VEMDataSet* p = new G4EMDataSet(iZ,energyVector,ksi,interp,1.,1.);
652 xsis->AddComponent(p);
653 }
654 if(verboseLevel>0) xsis->PrintData();
655 shellVacancy->AddXsiTable(xsis);
656 }
657}
658
659//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
660
661void G4LivermoreIonisationModel::ActivateAuger(G4bool val)
662{
663 deexcitationManager.ActivateAugerElectronProduction(val);
664}
665
Note: See TracBrowser for help on using the repository browser.