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

Last change on this file since 1058 was 1058, checked in by garnier, 15 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.