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

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

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