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

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

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