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

Last change on this file since 1228 was 1228, checked in by garnier, 14 years ago

update geant4.9.3 tag
File size: 23.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: G4LivermoreIonisationModel.cc,v 1.7 2009/10/23 09:30:08 pandola Exp $
27// GEANT4 tag $Name: geant4-09-03 $
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// 23 Oct 2009   L Pandola
44//                  - atomic deexcitation managed via G4VEmModel::DeexcitationFlag() is
45//                    set as "true" (default would be false)
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  //
88  //By default: use deexcitation, not auger
89  SetDeexcitationFlag(true);
90  ActivateAuger(false);
91
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{
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    }
673  deexcitationManager.ActivateAugerElectronProduction(val);
674  if (verboseLevel > 1)
675    G4cout << "Auger production set to " << val << G4endl;
676}
677
Note: See TracBrowser for help on using the repository browser.