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

Last change on this file since 991 was 968, checked in by garnier, 15 years ago

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