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

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

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