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

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

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