source: trunk/source/processes/electromagnetic/lowenergy/src/G4Penelope08IonisationModel.cc

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

update to last version 4.9.4
File size: 53.9 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: G4Penelope08IonisationModel.cc,v 1.2 2010/12/15 07:39:12 gunter Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 27 Jul 2010   L Pandola    First complete implementation
34//
35
36#include "G4Penelope08IonisationModel.hh"
37#include "G4ParticleDefinition.hh"
38#include "G4MaterialCutsCouple.hh"
39#include "G4ProductionCutsTable.hh"
40#include "G4DynamicParticle.hh"
41#include "G4AtomicTransitionManager.hh"
42#include "G4AtomicDeexcitation.hh"
43#include "G4AtomicShell.hh"
44#include "G4Gamma.hh"
45#include "G4Electron.hh"
46#include "G4Positron.hh"
47#include "G4AtomicDeexcitation.hh"
48#include "G4PenelopeOscillatorManager.hh"
49#include "G4PenelopeOscillator.hh"
50#include "G4PenelopeCrossSection.hh"
51#include "G4PhysicsFreeVector.hh"
52#include "G4PhysicsLogVector.hh"
53
54//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55 
56 
57G4Penelope08IonisationModel::G4Penelope08IonisationModel(const G4ParticleDefinition*,
58                                             const G4String& nam)
59  :G4VEmModel(nam),isInitialised(false),
60   kineticEnergy1(0.*eV),
61   cosThetaPrimary(1.0),
62   energySecondary(0.*eV),
63   cosThetaSecondary(0.0),
64   targetOscillator(-1),XSTableElectron(0),XSTablePositron(0),
65   theDeltaTable(0),energyGrid(0)
66{
67  fIntrinsicLowEnergyLimit = 100.0*eV;
68  fIntrinsicHighEnergyLimit = 100.0*GeV;
69  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
70  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
71  nBins = 200;
72  //
73  oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
74  //
75  verboseLevel= 0;
76   
77  // Verbosity scale:
78  // 0 = nothing
79  // 1 = warning for energy non-conservation
80  // 2 = details of energy budget
81  // 3 = calculation of cross sections, file openings, sampling of atoms
82  // 4 = entering in methods
83
84  // Atomic deexcitation model activated by default
85  SetDeexcitationFlag(true);
86  ActivateAuger(false);
87}
88
89//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 
91G4Penelope08IonisationModel::~G4Penelope08IonisationModel()
92{
93  ClearTables();
94}
95
96//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97void G4Penelope08IonisationModel::Initialise(const G4ParticleDefinition*,
98                                             const G4DataVector&)
99{
100  if (verboseLevel > 3)
101    G4cout << "Calling G4Penelope08IonisationModel::Initialise()" << G4endl;
102 
103  //Clear and re-build the tables
104  ClearTables();
105  XSTableElectron = new 
106    std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
107  XSTablePositron = new 
108    std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
109
110  theDeltaTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
111
112  //Set the number of bins for the tables. 20 points per decade
113  nBins = (size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit()));
114  nBins = std::max(nBins,(size_t)100);
115
116  energyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
117                                      HighEnergyLimit(), 
118                                      nBins-1); //one hidden bin is added
119
120  if (verboseLevel > 2) {
121    G4cout << "Penelope Ionisation model is initialized " << G4endl
122           << "Energy range: "
123           << LowEnergyLimit() / keV << " keV - "
124           << HighEnergyLimit() / GeV << " GeV. Using " 
125           << nBins << " bins." 
126           << G4endl;
127  }
128 
129  if(isInitialised) return;
130  fParticleChange = GetParticleChangeForLoss();
131  isInitialised = true;
132}
133
134//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135 
136G4double G4Penelope08IonisationModel::CrossSectionPerVolume(const G4Material* material,
137                                           const G4ParticleDefinition* theParticle,
138                                           G4double energy,
139                                           G4double cutEnergy,
140                                           G4double)
141{
142  // Penelope model to calculate the cross section for inelastic collisions above the
143  // threshold. It makes use of the Generalised Oscillator Strength (GOS) model from
144  //  D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
145  //
146  // The total cross section is calculated analytically by taking
147  // into account the atomic oscillators coming into the play for a given threshold.
148  //
149  // For incident e- the maximum energy allowed for the delta-rays is energy/2.
150  // because particles are undistinghishable.
151  //
152  // The contribution is splitted in three parts: distant longitudinal collisions,
153  // distant transverse collisions and close collisions. Each term is described by
154  // its own analytical function.
155  // Fermi density correction is calculated analytically according to
156  //  U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
157  //
158  if (verboseLevel > 3)
159    G4cout << "Calling CrossSectionPerVolume() of G4Penelope08IonisationModel" << G4endl;
160 
161  SetupForMaterial(theParticle, material, energy);
162   
163  G4double totalCross = 0.0;
164  G4double crossPerMolecule = 0.;
165
166  G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
167                                                                cutEnergy);
168
169  if (theXS)
170    crossPerMolecule = theXS->GetHardCrossSection(energy);
171
172  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
173  G4double atPerMol =  oscManager->GetAtomsPerMolecule(material);
174 
175  if (verboseLevel > 3)
176    G4cout << "Material " << material->GetName() << " has " << atPerMol <<
177      "atoms per molecule" << G4endl;
178
179  G4double moleculeDensity = 0.; 
180  if (atPerMol)
181    moleculeDensity = atomDensity/atPerMol;
182 
183  G4double crossPerVolume = crossPerMolecule*moleculeDensity;
184
185  if (verboseLevel > 2)
186  {
187    G4cout << "G4Penelope08IonisationModel " << G4endl;
188    G4cout << "Mean free path for delta emission > " << cutEnergy/keV << " keV at " <<
189      energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl;
190    if (theXS)
191      totalCross = (theXS->GetTotalCrossSection(energy))*moleculeDensity;
192    G4cout << "Total free path for ionisation (no threshold) at " <<
193      energy/keV << " keV = " << (1./totalCross)/mm << " mm" << G4endl;
194  }
195  return crossPerVolume;
196}
197
198//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
199                                                                                                     
200//This is a dummy method. Never inkoved by the tracking, it just issues
201//a warning if one tries to get Cross Sections per Atom via the
202//G4EmCalculator.
203G4double G4Penelope08IonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
204                                                             G4double,
205                                                             G4double,
206                                                             G4double,
207                                                             G4double,
208                                                             G4double)
209{
210  G4cout << "*** G4Penelope08IonisationModel -- WARNING ***" << G4endl;
211  G4cout << "Penelope08 Ionisation model does not calculate cross section _per atom_ " << G4endl;
212  G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
213  G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
214  return 0;
215}
216 
217//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
218
219G4double G4Penelope08IonisationModel::ComputeDEDXPerVolume(const G4Material* material,
220                                           const G4ParticleDefinition* theParticle,
221                                           G4double kineticEnergy,
222                                           G4double cutEnergy)
223{
224  // Penelope model to calculate the stopping power for soft inelastic collisions
225  // below the threshold. It makes use of the Generalised Oscillator Strength (GOS)
226  // model from
227  //  D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
228  //
229  // The stopping power is calculated analytically using the dsigma/dW cross
230  // section from the GOS models, which includes separate contributions from
231  // distant longitudinal collisions, distant transverse collisions and
232  // close collisions. Only the atomic oscillators that come in the play
233  // (according to the threshold) are considered for the calculation.
234  // Differential cross sections have a different form for e+ and e-.
235  //
236  // Fermi density correction is calculated analytically according to
237  //  U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
238 
239  if (verboseLevel > 3)
240    G4cout << "Calling ComputeDEDX() of G4Penelope08IonisationModel" << G4endl;
241 
242
243  G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
244                                                                cutEnergy);
245  G4double sPowerPerMolecule = 0.0;
246  if (theXS)
247    sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy);
248
249 
250  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
251  G4double atPerMol =  oscManager->GetAtomsPerMolecule(material);
252                                                                                       
253  G4double moleculeDensity = 0.; 
254  if (atPerMol)
255    moleculeDensity = atomDensity/atPerMol;
256 
257  G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
258
259  if (verboseLevel > 2)
260    {
261      G4cout << "G4Penelope08IonisationModel " << G4endl;
262      G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
263        kineticEnergy/keV << " keV = " << 
264        sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl;
265    }
266  return sPowerPerVolume;
267}
268 
269//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
270
271G4double G4Penelope08IonisationModel::MinEnergyCut(const G4ParticleDefinition*,
272                                                 const G4MaterialCutsCouple*)
273{
274  return fIntrinsicLowEnergyLimit;
275}
276
277//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
278
279void G4Penelope08IonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
280                                              const G4MaterialCutsCouple* couple,
281                                              const G4DynamicParticle* aDynamicParticle,
282                                              G4double cutE, G4double)
283{
284  // Penelope model to sample the final state following an hard inelastic interaction.
285  // It makes use of the Generalised Oscillator Strength (GOS) model from
286  //  D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
287  //
288  // The GOS model is used to calculate the individual cross sections for all
289  // the atomic oscillators coming in the play, taking into account the three
290  // contributions (distant longitudinal collisions, distant transverse collisions and
291  // close collisions). Then the target shell and the interaction channel are
292  // sampled. Final state of the delta-ray (energy, angle) are generated according
293  // to the analytical distributions (dSigma/dW) for the selected interaction
294  // channels.
295  // For e-, the maximum energy for the delta-ray is initialEnergy/2. (because
296  // particles are indistinghusbable), while it is the full initialEnergy for
297  // e+.
298  // The efficiency of the random sampling algorithm (e.g. for close collisions)
299  // decreases when initial and cutoff energy increase (e.g. 87% for 10-keV primary
300  // and 1 keV threshold, 99% for 10-MeV primary and 10-keV threshold).
301  // Differential cross sections have a different form for e+ and e-.
302  //
303  // WARNING: The model provides an _average_ description of inelastic collisions.
304  // Anyway, the energy spectrum associated to distant excitations of a given
305  // atomic shell is approximated as a single resonance. The simulated energy spectra
306  // show _unphysical_ narrow peaks at energies that are multiple of the shell
307  // resonance energies. The spurious speaks are automatically smoothed out after
308  // multiple inelastic collisions.
309  //
310  // The model determines also the original shell from which the delta-ray is expelled,
311  // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
312  //
313  // Fermi density correction is calculated analytically according to
314  //  U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
315
316  if (verboseLevel > 3)
317    G4cout << "Calling SamplingSecondaries() of G4Penelope08IonisationModel" << G4endl;
318 
319  G4double kineticEnergy0 = aDynamicParticle->GetKineticEnergy();
320  const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition();
321 
322  if (kineticEnergy0 <= fIntrinsicLowEnergyLimit)
323  {
324    fParticleChange->SetProposedKineticEnergy(0.);
325    fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy0);
326    return ;
327  }
328  const G4Material* material = couple->GetMaterial();
329  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(material); 
330
331  G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
332 
333  //Initialise final-state variables. The proper values will be set by the methods
334  // SampleFinalStateElectron() and SampleFinalStatePositron()
335  kineticEnergy1=kineticEnergy0;
336  cosThetaPrimary=1.0;
337  energySecondary=0.0;
338  cosThetaSecondary=1.0;
339  targetOscillator = -1;
340     
341  if (theParticle == G4Electron::Electron())
342    SampleFinalStateElectron(material,cutE,kineticEnergy0);
343  else if (theParticle == G4Positron::Positron())
344    SampleFinalStatePositron(material,cutE,kineticEnergy0);
345  else
346    {
347      G4cout << "G4Penelope08IonisationModel::SamplingSecondaries() - " ;
348      G4cout << "Invalid particle " << theParticle->GetParticleName() << G4endl;
349      G4Exception();   
350    }
351 
352  if (verboseLevel > 3)
353  {
354     G4cout << "G4Penelope08IonisationModel::SamplingSecondaries() for " << 
355        theParticle->GetParticleName() << G4endl;
356      G4cout << "Final eKin = " << kineticEnergy1 << " keV" << G4endl;
357      G4cout << "Final cosTheta = " << cosThetaPrimary << G4endl;
358      G4cout << "Delta-ray eKin = " << energySecondary << " keV" << G4endl;
359      G4cout << "Delta-ray cosTheta = " << cosThetaSecondary << G4endl;
360      G4cout << "Oscillator: " << targetOscillator << G4endl;
361   }
362   
363  //Update the primary particle
364  G4double sint = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary);
365  G4double phiPrimary  = twopi * G4UniformRand();
366  G4double dirx = sint * std::cos(phiPrimary);
367  G4double diry = sint * std::sin(phiPrimary);
368  G4double dirz = cosThetaPrimary;
369 
370  G4ThreeVector electronDirection1(dirx,diry,dirz);
371  electronDirection1.rotateUz(particleDirection0);
372   
373  if (kineticEnergy1 > 0)
374    {
375      fParticleChange->ProposeMomentumDirection(electronDirection1);
376      fParticleChange->SetProposedKineticEnergy(kineticEnergy1);
377    }
378  else
379    {
380      fParticleChange->SetProposedKineticEnergy(0.);
381    }
382 
383  //Generate the delta ray
384  G4double ionEnergyInPenelopeDatabase = 
385    (*theTable)[targetOscillator]->GetIonisationEnergy();
386 
387  //Now, try to handle fluorescence
388  //Notice: merged levels are indicated with Z=0 and flag=30
389  G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag(); 
390  G4int Z = (G4int) (*theTable)[targetOscillator]->GetParentZ();
391
392  //initialize here, then check photons created by Atomic-Deexcitation, and the final state e-
393  std::vector<G4DynamicParticle*>* photonVector=0; 
394  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
395  G4double bindingEnergy = 0.*eV;
396  G4int shellId = 0;
397
398  //Real level
399  if (Z > 0 && shFlag<30)
400    {
401      const G4AtomicShell* shell = transitionManager->Shell(Z,shFlag-1);
402      bindingEnergy = shell->BindingEnergy();
403      shellId = shell->ShellId();
404    }
405
406  //correct the energySecondary to account for the fact that the Penelope
407  //database of ionisation energies is in general (slightly) different
408  //from the fluorescence database used in Geant4.
409  energySecondary += ionEnergyInPenelopeDatabase-bindingEnergy;
410 
411  G4double localEnergyDeposit = bindingEnergy;
412  //testing purposes only
413  G4double energyInFluorescence = 0;
414
415  if (energySecondary < 0)
416    {
417      //It means that there was some problem/mismatch between the two databases.
418      //In this case, the available energy is ok to excite the level according
419      //to the Penelope database, but not according to the Geant4 database
420      //Full residual energy is deposited locally
421      localEnergyDeposit += energySecondary;
422      energySecondary = 0.0;
423    }
424 
425  //the local energy deposit is what remains: part of this may be spent for fluorescence.
426  if(DeexcitationFlag() && Z > 5) 
427    {
428      const G4ProductionCutsTable* theCoupleTable=
429        G4ProductionCutsTable::GetProductionCutsTable();
430
431      size_t index = couple->GetIndex();
432      G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
433      G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
434
435      // Generation of fluorescence
436      // Data in EADL are available only for Z > 5
437      // Protection to avoid generating photons in the unphysical case of
438      // shell binding energy > photon energy
439      if (localEnergyDeposit > cutg || localEnergyDeposit > cute)
440        { 
441          G4DynamicParticle* aPhoton;
442          deexcitationManager.SetCutForSecondaryPhotons(cutg);
443          deexcitationManager.SetCutForAugerElectrons(cute);
444         
445          photonVector = deexcitationManager.GenerateParticles(Z,shellId);
446          if(photonVector) 
447            {
448              size_t nPhotons = photonVector->size();
449              for (size_t k=0; k<nPhotons; k++)
450                {
451                  aPhoton = (*photonVector)[k];
452                  if (aPhoton)
453                    {
454                      G4double itsEnergy = aPhoton->GetKineticEnergy();
455                      if (itsEnergy <= localEnergyDeposit)
456                        {
457                        localEnergyDeposit -= itsEnergy;
458                        if (aPhoton->GetDefinition() == G4Gamma::Gamma()) 
459                          energyInFluorescence += itsEnergy;;
460                        fvect->push_back(aPhoton);                 
461                        }
462                      else
463                        {
464                          delete aPhoton;
465                          (*photonVector)[k]=0;
466                        }
467                    }
468                }
469              delete photonVector;
470            }
471        }
472    }
473
474  //Always produce explicitely the delta electron
475  G4DynamicParticle* electron = 0;
476  G4double sinThetaE = std::sqrt(1.-cosThetaSecondary*cosThetaSecondary);
477  G4double phiEl = phiPrimary+pi; //pi with respect to the primary electron/positron
478  G4double xEl = sinThetaE * std::cos(phiEl);
479  G4double yEl = sinThetaE * std::sin(phiEl);
480  G4double zEl = cosThetaSecondary;
481  G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
482  eDirection.rotateUz(particleDirection0);
483  electron = new G4DynamicParticle (G4Electron::Electron(),
484                                    eDirection,energySecondary) ;
485  fvect->push_back(electron);
486
487
488  if (localEnergyDeposit < 0)
489    {
490      G4cout << "WARNING-" 
491             << "G4Penelope08IonisationModel::SampleSecondaries - Negative energy deposit"
492             << G4endl;
493      localEnergyDeposit=0.;
494    }
495  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
496
497  if (verboseLevel > 1)
498    {
499      G4cout << "-----------------------------------------------------------" << G4endl;
500      G4cout << "Energy balance from G4Penelope08Ionisation" << G4endl;
501      G4cout << "Incoming primary energy: " << kineticEnergy0/keV << " keV" << G4endl;
502      G4cout << "-----------------------------------------------------------" << G4endl;
503      G4cout << "Outgoing primary energy: " << kineticEnergy1/keV << " keV" << G4endl;
504      G4cout << "Delta ray " << energySecondary/keV << " keV" << G4endl;
505      G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
506      G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
507      G4cout << "Total final state: " << (energySecondary+energyInFluorescence+kineticEnergy1+
508                                              localEnergyDeposit)/keV <<
509        " keV" << G4endl;
510      G4cout << "-----------------------------------------------------------" << G4endl;
511    }
512     
513  if (verboseLevel > 0)
514    {
515      G4double energyDiff = std::fabs(energySecondary+energyInFluorescence+kineticEnergy1+
516                                      localEnergyDeposit-kineticEnergy0);
517      if (energyDiff > 0.05*keV)
518        G4cout << "Warning from G4PenelopeIonisation: problem with energy conservation: " << 
519          (energySecondary+energyInFluorescence+kineticEnergy1+localEnergyDeposit)/keV <<
520          " keV (final) vs. " <<
521          kineticEnergy0/keV << " keV (initial)" << G4endl;
522    }
523   
524}
525 
526
527//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
528 
529void G4Penelope08IonisationModel::ActivateAuger(G4bool augerbool)
530{
531  if (!DeexcitationFlag() && augerbool)
532    {
533      G4cout << "WARNING - G4Penelope08IonisationModel" << G4endl;
534      G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
535      G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
536    }
537  deexcitationManager.ActivateAugerElectronProduction(augerbool);
538  if (verboseLevel > 1)
539    G4cout << "Auger production set to " << augerbool << G4endl;
540}
541 
542//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
543
544void G4Penelope08IonisationModel::ClearTables()
545{ 
546  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>::iterator i;
547  if (XSTableElectron)
548    {
549      for (i=XSTableElectron->begin(); i != XSTableElectron->end(); i++)
550        {
551          G4PenelopeCrossSection* tab = i->second;
552          delete tab;
553        }
554      delete XSTableElectron;
555      XSTableElectron = 0;
556    }
557
558  if (XSTablePositron)
559    {
560      for (i=XSTablePositron->begin(); i != XSTablePositron->end(); i++)
561        {
562          G4PenelopeCrossSection* tab = i->second;
563          delete tab;
564        }
565      delete XSTablePositron;
566      XSTablePositron = 0;
567    }
568
569  std::map<const G4Material*,G4PhysicsFreeVector*>::iterator k;
570  if (theDeltaTable)
571    {
572      for (k=theDeltaTable->begin();k!=theDeltaTable->end();k++)       
573        delete k->second;
574      delete theDeltaTable;
575      theDeltaTable = 0;
576    }
577   
578  if (energyGrid)
579    delete energyGrid;
580
581  if (verboseLevel > 2)
582    G4cout << "G4Penelope08IonisationModel: cleared tables" << G4endl;
583  return;
584}
585
586//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
587
588G4PenelopeCrossSection* 
589G4Penelope08IonisationModel::GetCrossSectionTableForCouple(const G4ParticleDefinition* part,
590                                                           const G4Material* mat,
591                                                           G4double cut)
592{
593  if (part != G4Electron::Electron() && part != G4Positron::Positron())
594    {
595      G4cout << "G4Penelope08IonisationModel::GetCrossSectionTableForCouple" << G4endl;
596      G4cout << "Invalid particle: " << part->GetParticleName() << G4endl;
597      G4Exception();
598      return NULL;
599    }
600
601  if (part == G4Electron::Electron())
602    {
603      if (!XSTableElectron)
604        {
605          G4cout << "G4Penelope08IonisationModel::GetCrossSectionTableForCouple()" << G4endl;
606          G4cout << "The Cross Section Table for e- was not initialized correctly!" << G4endl;
607          G4Exception();         
608        }
609      std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
610      if (XSTableElectron->count(theKey)) //table already built
611        return XSTableElectron->find(theKey)->second;
612      else
613        {
614          BuildXSTable(mat,cut,part);
615          if (XSTableElectron->count(theKey)) //now it should be ok!
616            return XSTableElectron->find(theKey)->second;
617          else
618            {
619              G4cout << "G4Penelope08IonisationModel::GetCrossSectionTableForCouple()" << G4endl;
620              G4cout << "Unable to build e- table for " << mat->GetName() << G4endl;
621              G4Exception();       
622            }
623        }
624    }
625
626  if (part == G4Positron::Positron())
627    {
628      if (!XSTablePositron)
629        {
630          G4cout << "G4Penelope08IonisationModel::GetCrossSectionTableForCouple()" << G4endl;
631          G4cout << "The Cross Section Table for e+ was not initialized correctly!" << G4endl;
632          G4Exception();         
633        }
634      std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
635      if (XSTablePositron->count(theKey)) //table already built
636        return XSTablePositron->find(theKey)->second;
637      else
638        {
639          BuildXSTable(mat,cut,part);
640          if (XSTablePositron->count(theKey)) //now it should be ok!
641            return XSTablePositron->find(theKey)->second;
642          else
643            {
644              G4cout << "G4Penelope08IonisationModel::GetCrossSectionTableForCouple()" << G4endl;
645              G4cout << "Unable to build e+ table for " << mat->GetName() << G4endl;
646              G4Exception();       
647            }
648        }
649    }
650  return NULL;
651}
652
653//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
654
655void G4Penelope08IonisationModel::BuildXSTable(const G4Material* mat,G4double cut,
656                                               const G4ParticleDefinition* part)
657{
658  //
659  //This method fills the G4PenelopeCrossSection containers for electrons or positrons
660  //and for the given material/cut couple. The calculation is done as sum over the
661  //individual shells.
662  //Equivalent of subroutines EINaT and PINaT of Penelope
663  //
664  if (verboseLevel > 2)
665    {
666      G4cout << "G4Penelope08IonisationModel: going to build cross section table " << G4endl;
667      G4cout << "for " << part->GetParticleName() << " in " << mat->GetName() << G4endl;
668    }
669
670  //Tables have been already created (checked by GetCrossSectionTableForCouple)
671  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
672  size_t numberOfOscillators = theTable->size();
673
674  if (energyGrid->GetVectorLength() != nBins) 
675    {
676      G4cout << "G4Penelope08IonisationModel::BuildXSTable" << G4endl;
677      G4cout << "Energy Grid looks not initialized" << G4endl;
678      G4cout << nBins << " " << energyGrid->GetVectorLength() << G4endl;
679      G4Exception();
680    }
681
682  G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(nBins,numberOfOscillators);
683 
684  //loop on the energy grid
685  for (size_t bin=0;bin<nBins;bin++)
686    {
687       G4double energy = energyGrid->GetLowEdgeEnergy(bin);
688       G4double XH0=0, XH1=0, XH2=0;
689       G4double XS0=0, XS1=0, XS2=0;
690   
691       //oscillator loop
692       for (size_t iosc=0;iosc<numberOfOscillators;iosc++)
693         {
694           G4DataVector* tempStorage = 0;
695
696           G4PenelopeOscillator* theOsc = (*theTable)[iosc];
697           G4double delta = GetDensityCorrection(mat,energy);
698           if (part == G4Electron::Electron())       
699             tempStorage = ComputeShellCrossSectionsElectron(theOsc,energy,cut,delta);
700           else if (part == G4Positron::Positron())
701             tempStorage = ComputeShellCrossSectionsPositron(theOsc,energy,cut,delta);
702           //check results are all right
703           if (!tempStorage)
704             {
705               G4cout << "G4Penelope08IonisationModel::BuildXSTable" << G4endl;
706               G4cout << "Problem in calculating the shell XS " << G4endl;
707               G4Exception();
708             }
709           if (tempStorage->size() != 6)
710             {
711               G4cout << "G4Penelope08IonisationModel::BuildXSTable" << G4endl;
712               G4cout << "Problem in calculating the shell XS " << G4endl;
713               G4cout << "Result has dimension " << tempStorage->size() << " instead of 6" << G4endl;
714               G4Exception();
715             }
716           G4double stre = theOsc->GetOscillatorStrength();
717
718           XH0 += stre*(*tempStorage)[0];
719           XH1 += stre*(*tempStorage)[1];
720           XH2 += stre*(*tempStorage)[2];
721           XS0 += stre*(*tempStorage)[3];
722           XS1 += stre*(*tempStorage)[4];
723           XS2 += stre*(*tempStorage)[5];
724           XSEntry->AddShellCrossSectionPoint(bin,iosc,energy,stre*(*tempStorage)[0]);
725           if (tempStorage)
726             {
727               delete tempStorage;
728               tempStorage = 0;
729             }
730         }       
731       XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2);
732    }
733
734
735  //Normalize shell cross sections before insertion in the table
736  XSEntry->NormalizeShellCrossSections();
737
738
739  //Insert in the appropriate table
740  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
741  if (part == G4Electron::Electron()) 
742    XSTableElectron->insert(std::make_pair(theKey,XSEntry));
743  else if (part == G4Positron::Positron())
744    XSTablePositron->insert(std::make_pair(theKey,XSEntry));
745 
746  return;
747}
748
749
750//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751
752G4double G4Penelope08IonisationModel::GetDensityCorrection(const G4Material* mat,
753                                                           G4double energy)
754{
755  G4double result = 0;
756  if (!theDeltaTable)
757    {
758      G4cout << "G4Penelope08IonisationModel::GetDensityCorrection()" << G4endl;
759      G4cout << "Delta Table not initialized. Was Initialise() run?" << G4endl;
760      G4Exception();
761    }
762  if (energy <= 0*eV)
763    {
764      G4cout << "G4Penelope08IonisationModel::GetDensityCorrection()" << G4endl;
765      G4cout << "Invalid energy " << energy/eV << " eV " << G4endl;
766      return 0;
767    }
768  G4double logene = std::log(energy);
769
770  //check if the material has been built
771  if (!(theDeltaTable->count(mat)))
772    BuildDeltaTable(mat);
773
774  if (theDeltaTable->count(mat))
775    {
776      G4PhysicsFreeVector* vec = theDeltaTable->find(mat)->second;
777      result = vec->Value(logene); //the table has delta vs. ln(E)     
778    }
779  else
780    {
781      G4cout << "G4Penelope08IonisationModel::GetDensityCorrection()" << G4endl;
782      G4cout << "Unable to build table for " << mat->GetName() << G4endl;
783      G4Exception();
784    }
785
786  return result;
787}
788
789//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
790
791void G4Penelope08IonisationModel::BuildDeltaTable(const G4Material* mat)
792{
793  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
794  G4double plasmaSq = oscManager->GetPlasmaEnergySquared(mat);
795  G4double totalZ = oscManager->GetTotalZ(mat);
796  size_t numberOfOscillators = theTable->size();
797
798  if (energyGrid->GetVectorLength() != nBins) 
799    {
800      G4cout << "G4Penelope08IonisationModel::BuildDeltaTable" << G4endl;
801      G4cout << "Energy Grid for Delta looks not initialized" << G4endl;
802      G4cout << nBins << " " << energyGrid->GetVectorLength() << G4endl;
803      G4Exception();
804    }
805
806  G4PhysicsFreeVector* theVector = new G4PhysicsFreeVector(nBins);
807
808  //loop on the energy grid
809  for (size_t bin=0;bin<nBins;bin++)
810    {
811      G4double delta = 0.;
812      G4double energy = energyGrid->GetLowEdgeEnergy(bin);
813
814      //Here calculate delta
815      G4double gam = 1.0+(energy/electron_mass_c2);
816      G4double gamSq = gam*gam;
817
818      G4double TST = totalZ/(gamSq*plasmaSq);
819      G4double wl2 = 0;
820      G4double fdel = 0;
821
822      //loop on oscillators
823      for (size_t i=0;i<numberOfOscillators;i++)
824        {
825          G4PenelopeOscillator* theOsc = (*theTable)[i];
826          G4double wri = theOsc->GetResonanceEnergy();
827          fdel += theOsc->GetOscillatorStrength()/(wri*wri+wl2);
828        }     
829      if (fdel >= TST) //if fdel < TST, delta = 0
830        {
831          //get last oscillator
832          G4PenelopeOscillator* theOsc = (*theTable)[numberOfOscillators-1]; 
833          wl2 = theOsc->GetResonanceEnergy()*theOsc->GetResonanceEnergy();
834         
835          //First iteration
836          G4bool loopAgain = false;
837          do
838            {
839              loopAgain = false;
840              wl2 += wl2;
841              fdel = 0.;
842              for (size_t i=0;i<numberOfOscillators;i++)
843                {
844                  G4PenelopeOscillator* theOsc = (*theTable)[i];
845                  G4double wri = theOsc->GetResonanceEnergy();
846                  fdel += theOsc->GetOscillatorStrength()/(wri*wri+wl2);
847                }
848              if (fdel > TST)
849                loopAgain = true;
850            }while(loopAgain);
851
852          G4double wl2l = 0;
853          G4double wl2u = wl2;
854          //second iteration
855          do
856            {       
857              loopAgain = false;
858              wl2 = 0.5*(wl2l+wl2u);
859              fdel = 0;
860              for (size_t i=0;i<numberOfOscillators;i++)
861                {
862                  G4PenelopeOscillator* theOsc = (*theTable)[i];
863                  G4double wri = theOsc->GetResonanceEnergy();
864                  fdel += theOsc->GetOscillatorStrength()/(wri*wri+wl2);
865                }
866              if (fdel > TST)
867                wl2l = wl2;
868              else
869                wl2u = wl2;
870              if ((wl2u-wl2l)>1e-12*wl2)
871                loopAgain = true;
872            }while(loopAgain);
873         
874          //Eventually get density correction
875          delta = 0.;
876          for (size_t i=0;i<numberOfOscillators;i++)
877            {
878              G4PenelopeOscillator* theOsc = (*theTable)[i];
879              G4double wri = theOsc->GetResonanceEnergy();
880              delta += theOsc->GetOscillatorStrength()*
881                std::log(1.0+(wl2/(wri*wri)));               
882            }
883          delta = (delta/totalZ)-wl2/(gamSq*plasmaSq);
884        }
885      energy = std::max(1e-9*eV,energy); //prevents log(0)
886      theVector->PutValue(bin,std::log(energy),delta);
887    }
888  theDeltaTable->insert(std::make_pair(mat,theVector));
889  return;
890}
891                                                       
892//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
893G4DataVector* G4Penelope08IonisationModel::ComputeShellCrossSectionsElectron(G4PenelopeOscillator* theOsc,
894                                                                             G4double energy,
895                                                                             G4double cut,
896                                                                             G4double delta)
897{
898  //
899  //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for
900  //the given oscillator/cut and at the given energy.
901  //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2)
902  //Equivalent of subroutines EINaT1 of Penelope
903  //
904  // Results are _per target electron_
905  //
906  G4DataVector* result = new G4DataVector();
907  for (size_t i=0;i<6;i++)
908    result->push_back(0.);
909  G4double ionEnergy = theOsc->GetIonisationEnergy();
910 
911  //return a set of zero's if the energy it too low to excite the current oscillator
912  if (energy < ionEnergy)
913    return result;
914
915  G4double H0=0.,H1=0.,H2=0.;
916  G4double S0=0.,S1=0.,S2=0.;
917
918  //Define useful constants to be used in the calculation
919  G4double gamma = 1.0+energy/electron_mass_c2;
920  G4double gammaSq = gamma*gamma;
921  G4double beta = (gammaSq-1.0)/gammaSq;
922  G4double pielr2 = pi*classic_electr_radius*classic_electr_radius; //pi*re^2
923  G4double constant = pielr2*2.0*electron_mass_c2/beta;
924  G4double XHDT0 = std::log(gammaSq)-beta;
925
926  G4double cpSq = energy*(energy+2.0*electron_mass_c2);
927  G4double cp = std::sqrt(cpSq);
928  G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2));
929
930  //
931  // Distant interactions
932  //
933  G4double resEne = theOsc->GetResonanceEnergy();
934  G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy();
935  if (energy > resEne)
936    {
937      G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2);
938      G4double cp1 = std::sqrt(cp1Sq);
939     
940      //Distant longitudinal interactions
941      G4double QM = 0;
942      if (resEne > 1e-6*energy)
943        QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
944      else
945        {
946          QM = resEne*resEne/(beta*2.0*electron_mass_c2);
947          QM = QM*(1.0-0.5*QM/electron_mass_c2);
948        }
949      G4double SDL1 = 0;
950      if (QM < cutoffEne)
951        SDL1 = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)));
952     
953      //Distant transverse interactions
954      if (SDL1)
955        {
956          G4double SDT1 = std::max(XHDT0-delta,0.0);
957          G4double SD1 = SDL1+SDT1;
958          if (cut > resEne)
959            {
960              S1 = SD1; //XS1
961              S0 = SD1/resEne; //XS0
962              S2 = SD1*resEne; //XS2
963            }
964          else
965            {
966              H1 = SD1; //XH1
967              H0 = SD1/resEne; //XH0
968              H2 = SD1*resEne; //XH2
969            }
970        }
971    }
972  //
973  // Close collisions (Moller's cross section)
974  //
975  G4double wl = std::max(cut,cutoffEne);
976  G4double ee = energy + ionEnergy;
977  G4double wu = 0.5*ee;
978  if (wl < wu-(1e-5*eV))
979    {
980      H0 += (1.0/(ee-wu)) - (1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) + 
981        (1.0-amol)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee + 
982        amol*(wu-wl)/(ee*ee);
983      H1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) + 
984        (2.0-amol)*std::log((ee-wu)/(ee-wl)) + 
985        amol*(wu*wu-wl*wl)/(2.0*ee*ee);
986      H2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) - 
987        (wl*(2.0*ee-wl)/(ee-wl)) + 
988        (3.0-amol)*ee*std::log((ee-wu)/(ee-wl)) +       
989        amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee);
990      wu = wl;
991    }
992  wl = cutoffEne;
993 
994  if (wl > wu-(1e-5*eV))
995    {
996      (*result)[0] = constant*H0;
997      (*result)[1] = constant*H1;
998      (*result)[2] = constant*H2;
999      (*result)[3] = constant*S0;
1000      (*result)[4] = constant*S1;
1001      (*result)[5] = constant*S2;
1002      return result;
1003    }
1004
1005  S0 += (1.0/(ee-wu))-(1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) + 
1006    (1.0-amol)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee +
1007    amol*(wu-wl)/(ee*ee);
1008  S1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) + 
1009    (2.0-amol)*std::log((ee-wu)/(ee-wl)) + 
1010    amol*(wu*wu-wl*wl)/(2.0*ee*ee);
1011  S2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) - 
1012    (wl*(2.0*ee-wl)/(ee-wl)) + 
1013    (3.0-amol)*ee*std::log((ee-wu)/(ee-wl)) + 
1014    amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee);
1015
1016  (*result)[0] = constant*H0;
1017  (*result)[1] = constant*H1;
1018  (*result)[2] = constant*H2;
1019  (*result)[3] = constant*S0;
1020  (*result)[4] = constant*S1;
1021  (*result)[5] = constant*S2;
1022  return result;
1023}
1024
1025//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1026G4DataVector* G4Penelope08IonisationModel::ComputeShellCrossSectionsPositron(G4PenelopeOscillator* theOsc,
1027                                                                             G4double energy,
1028                                                                             G4double cut,
1029                                                                             G4double delta)
1030{
1031  //
1032  //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for
1033  //the given oscillator/cut and at the given energy.
1034  //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2)
1035  //Equivalent of subroutines PINaT1 of Penelope
1036  //
1037  // Results are _per target electron_
1038  //
1039  G4DataVector* result = new G4DataVector();
1040  for (size_t i=0;i<6;i++)
1041    result->push_back(0.);
1042  G4double ionEnergy = theOsc->GetIonisationEnergy();
1043 
1044  //return a set of zero's if the energy it too low to excite the current oscillator
1045  if (energy < ionEnergy)
1046    return result;
1047
1048  G4double H0=0.,H1=0.,H2=0.;
1049  G4double S0=0.,S1=0.,S2=0.;
1050
1051  //Define useful constants to be used in the calculation
1052  G4double gamma = 1.0+energy/electron_mass_c2;
1053  G4double gammaSq = gamma*gamma;
1054  G4double beta = (gammaSq-1.0)/gammaSq;
1055  G4double pielr2 = pi*classic_electr_radius*classic_electr_radius; //pi*re^2
1056  G4double constant = pielr2*2.0*electron_mass_c2/beta;
1057  G4double XHDT0 = std::log(gammaSq)-beta;
1058
1059  G4double cpSq = energy*(energy+2.0*electron_mass_c2);
1060  G4double cp = std::sqrt(cpSq);
1061  G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2));
1062  G4double g12 = (gamma+1.0)*(gamma+1.0);
1063  //Bhabha coefficients
1064  G4double bha1 = amol*(2.0*g12-1.0)/(gammaSq-1.0);
1065  G4double bha2 = amol*(3.0+1.0/g12);
1066  G4double bha3 = amol*2.0*gamma*(gamma-1.0)/g12;
1067  G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/g12;
1068
1069  //
1070  // Distant interactions
1071  //
1072  G4double resEne = theOsc->GetResonanceEnergy();
1073  G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy();
1074  if (energy > resEne)
1075    {
1076      G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2);
1077      G4double cp1 = std::sqrt(cp1Sq);
1078     
1079      //Distant longitudinal interactions
1080      G4double QM = 0;
1081      if (resEne > 1e-6*energy)
1082        QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
1083      else
1084        {
1085          QM = resEne*resEne/(beta*2.0*electron_mass_c2);
1086          QM = QM*(1.0-0.5*QM/electron_mass_c2);
1087        }
1088      G4double SDL1 = 0;
1089      if (QM < cutoffEne)
1090        SDL1 = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)));
1091     
1092      //Distant transverse interactions
1093      if (SDL1)
1094        {
1095          G4double SDT1 = std::max(XHDT0-delta,0.0);
1096          G4double SD1 = SDL1+SDT1;
1097          if (cut > resEne)
1098            {
1099              S1 = SD1; //XS1
1100              S0 = SD1/resEne; //XS0
1101              S2 = SD1*resEne; //XS2
1102            }
1103          else
1104            {
1105              H1 = SD1; //XH1
1106              H0 = SD1/resEne; //XH0
1107              H2 = SD1*resEne; //XH2
1108            }
1109        }
1110    }
1111
1112  //
1113  // Close collisions (Bhabha's cross section)
1114  //
1115  G4double wl = std::max(cut,cutoffEne);
1116  G4double wu = energy; 
1117  G4double energySq = energy*energy;
1118  if (wl < wu-(1e-5*eV))
1119    {
1120      G4double wlSq = wl*wl;
1121      G4double wuSq = wu*wu;
1122      H0 += (1.0/wl) - (1.0/wu)- bha1*std::log(wu/wl)/energy 
1123        + bha2*(wu-wl)/energySq 
1124        - bha3*(wuSq-wlSq)/(2.0*energySq*energy)
1125        + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq);
1126      H1 += std::log(wu/wl) - bha1*(wu-wl)/energy
1127        + bha2*(wuSq-wlSq)/(2.0*energySq)
1128        - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energy)
1129        + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq);
1130      H2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy)
1131        + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
1132        - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy)
1133        + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq);
1134      wu = wl;
1135    }
1136  wl = cutoffEne;
1137 
1138  if (wl > wu-(1e-5*eV))
1139    {
1140      (*result)[0] = constant*H0;
1141      (*result)[1] = constant*H1;
1142      (*result)[2] = constant*H2;
1143      (*result)[3] = constant*S0;
1144      (*result)[4] = constant*S1;
1145      (*result)[5] = constant*S2;
1146      return result;
1147    }
1148
1149  G4double wlSq = wl*wl;
1150  G4double wuSq = wu*wu;
1151
1152  S0 += (1.0/wl) - (1.0/wu) - bha1*std::log(wu/wl)/energy
1153    + bha2*(wu-wl)/energySq 
1154    - bha3*(wuSq-wlSq)/(2.0*energySq*energy)
1155    + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq);
1156
1157  S1 += std::log(wu/wl) - bha1*(wu-wl)/energy
1158    + bha2*(wuSq-wlSq)/(2.0*energySq)
1159    - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energy)
1160    + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq);
1161
1162  S2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy)
1163    + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
1164    - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy)
1165    + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq);
1166
1167 (*result)[0] = constant*H0;
1168 (*result)[1] = constant*H1;
1169 (*result)[2] = constant*H2;
1170 (*result)[3] = constant*S0;
1171 (*result)[4] = constant*S1;
1172 (*result)[5] = constant*S2;
1173
1174 return result;
1175}
1176
1177//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1178void G4Penelope08IonisationModel::SampleFinalStateElectron(const G4Material* mat,
1179                                                           G4double cutEnergy,
1180                                                           G4double kineticEnergy)
1181{
1182  // This method sets the final ionisation parameters
1183  // kineticEnergy1, cosThetaPrimary (= updates of the primary e-)
1184  // energySecondary, cosThetaSecondary (= info of delta-ray)
1185  // targetOscillator (= ionised oscillator)
1186  //
1187  // The method implements SUBROUTINE EINa of Penelope
1188  //
1189
1190  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
1191  size_t numberOfOscillators = theTable->size();
1192  G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(G4Electron::Electron(),mat,
1193                                                                cutEnergy);
1194  G4double delta = GetDensityCorrection(mat,kineticEnergy);
1195 
1196  // Selection of the active oscillator
1197  G4double TST = G4UniformRand();
1198  targetOscillator = numberOfOscillators-1; //initialization, last oscillator
1199  G4double XSsum = 0.;
1200
1201  for (size_t i=0;i<numberOfOscillators-1;i++)
1202    {
1203      XSsum += theXS->GetShellCrossSection(i,kineticEnergy); 
1204       
1205      if (XSsum > TST)
1206        {
1207          targetOscillator = (G4int) i;
1208          break;
1209        }
1210    }
1211
1212 
1213  if (verboseLevel > 3)
1214    {
1215      G4cout << "SampleFinalStateElectron: sampled oscillator #" << targetOscillator << "." << G4endl;
1216      G4cout << "Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/eV << 
1217        " eV " << G4endl;
1218      G4cout << "Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/eV << " eV "
1219             << G4endl;
1220    }
1221  //Constants
1222  G4double rb = kineticEnergy + 2.0*electron_mass_c2;
1223  G4double gam = 1.0+kineticEnergy/electron_mass_c2;
1224  G4double gam2 = gam*gam;
1225  G4double beta2 = (gam2-1.0)/gam2;
1226  G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam);
1227
1228  //Partial cross section of the active oscillator
1229  G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
1230  G4double invResEne = 1.0/resEne;
1231  G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
1232  G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
1233  G4double XHDL = 0.;
1234  G4double XHDT = 0.;
1235  G4double QM = 0.;
1236  G4double cps = 0.;
1237  G4double cp = 0.;
1238
1239  //Distant excitations
1240  if (resEne > cutEnergy && resEne < kineticEnergy)
1241    {
1242      cps = kineticEnergy*rb;
1243      cp = std::sqrt(cps);
1244      G4double XHDT0 = std::max(std::log(gam2)-beta2-delta,0.);
1245      if (resEne > 1.0e-6*kineticEnergy)
1246        {
1247          G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
1248          QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
1249        }
1250      else
1251        {
1252          QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
1253          QM *= (1.0-QM*0.5/electron_mass_c2);
1254        }
1255      if (QM < cutoffEne)
1256        {
1257          XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
1258            *invResEne;
1259          XHDT = XHDT0*invResEne;         
1260        }
1261      else
1262        {
1263          QM = cutoffEne;
1264          XHDL = 0.;
1265          XHDT = 0.;
1266        }
1267    }
1268  else
1269    {
1270      QM = cutoffEne;
1271      cps = 0.;
1272      cp = 0.;
1273      XHDL = 0.;
1274      XHDT = 0.;
1275    }
1276
1277  //Close collisions
1278  G4double EE = kineticEnergy + ionEne;
1279  G4double wmaxc = 0.5*EE;
1280  G4double wcl = std::max(cutEnergy,cutoffEne);
1281  G4double rcl = wcl/EE;
1282  G4double XHC = 0.;
1283  if (wcl < wmaxc)
1284    {
1285      G4double rl1 = 1.0-rcl;
1286      G4double rrl1 = 1.0/rl1;
1287      XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+
1288             (1.0-amol)*std::log(rcl*rrl1))/EE;
1289    }
1290
1291  //Total cross section per molecule for the active shell, in cm2
1292  G4double XHTOT = XHC + XHDL + XHDT;
1293
1294  //very small cross section, do nothing
1295  if (XHTOT < 1.e-14*barn)
1296    {
1297      kineticEnergy1=kineticEnergy;
1298      cosThetaPrimary=1.0;
1299      energySecondary=0.0;
1300      cosThetaSecondary=1.0;
1301      targetOscillator = numberOfOscillators-1;
1302      return;
1303    }
1304 
1305  //decide which kind of interaction we'll have
1306  TST = XHTOT*G4UniformRand();
1307 
1308
1309  // Hard close collision
1310  G4double TS1 = XHC;
1311 
1312  if (TST < TS1)
1313    {
1314      G4double A = 5.0*amol;
1315      G4double ARCL = A*0.5*rcl;
1316      G4double rk=0.;
1317      G4bool loopAgain = false;
1318      do
1319        {
1320          loopAgain = false;
1321          G4double fb = (1.0+ARCL)*G4UniformRand();
1322          if (fb < 1)
1323            rk = rcl/(1.0-fb*(1.0-(rcl+rcl)));       
1324          else
1325            rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL;
1326          G4double rk2 = rk*rk;
1327          G4double rkf = rk/(1.0-rk);
1328          G4double phi = 1.0+rkf*rkf-rkf+amol*(rk2+rkf);
1329          if (G4UniformRand()*(1.0+A*rk2) > phi)
1330            loopAgain = true;
1331        }while(loopAgain);
1332      //energy and scattering angle (primary electron)
1333      G4double deltaE = rk*EE;
1334     
1335      kineticEnergy1 = kineticEnergy - deltaE;
1336      cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
1337      //energy and scattering angle of the delta ray
1338      energySecondary = deltaE - ionEne; //subtract ionisation energy
1339      cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
1340      if (verboseLevel > 3)
1341        G4cout << "SampleFinalStateElectron: sampled close collision " << G4endl;
1342      return;                               
1343    }
1344
1345  //Hard distant longitudinal collisions
1346  TS1 += XHDL;
1347  G4double deltaE = resEne;
1348  kineticEnergy1 = kineticEnergy - deltaE;
1349
1350  if (TST < TS1)
1351    {
1352      G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
1353      G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand())
1354                       - (QS*0.5/electron_mass_c2));
1355      G4double QTREV = Q*(Q+2.0*electron_mass_c2);
1356      G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
1357      cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
1358      if (cosThetaPrimary > 1.) 
1359        cosThetaPrimary = 1.0;
1360      //energy and emission angle of the delta ray
1361      energySecondary = deltaE - ionEne;
1362      cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
1363      if (cosThetaSecondary > 1.0)
1364        cosThetaSecondary = 1.0;
1365      if (verboseLevel > 3)
1366        G4cout << "SampleFinalStateElectron: sampled distant longitudinal collision " << G4endl;
1367      return;     
1368    }
1369
1370  //Hard distant transverse collisions
1371  cosThetaPrimary = 1.0;
1372  //energy and emission angle of the delta ray
1373  energySecondary = deltaE - ionEne;
1374  cosThetaSecondary = 0.5;
1375  if (verboseLevel > 3)
1376    G4cout << "SampleFinalStateElectron: sampled distant transverse collision " << G4endl;
1377
1378  return;
1379}
1380
1381
1382
1383//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1384void G4Penelope08IonisationModel::SampleFinalStatePositron(const G4Material* mat,
1385                                                           G4double cutEnergy,
1386                                                           G4double kineticEnergy)
1387{
1388  // This method sets the final ionisation parameters
1389  // kineticEnergy1, cosThetaPrimary (= updates of the primary e-)
1390  // energySecondary, cosThetaSecondary (= info of delta-ray)
1391  // targetOscillator (= ionised oscillator)
1392  //
1393  // The method implements SUBROUTINE PINa of Penelope
1394  //
1395 
1396  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
1397  size_t numberOfOscillators = theTable->size();
1398  G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(G4Positron::Positron(),mat,
1399                                                                cutEnergy);
1400  G4double delta = GetDensityCorrection(mat,kineticEnergy);
1401
1402  // Selection of the active oscillator
1403  G4double TST = G4UniformRand();
1404  targetOscillator = numberOfOscillators-1; //initialization, last oscillator
1405  G4double XSsum = 0.;
1406  for (size_t i=0;i<numberOfOscillators-1;i++)
1407    {
1408      XSsum += theXS->GetShellCrossSection(i,kineticEnergy);   
1409      if (XSsum > TST)
1410        {
1411          targetOscillator = (G4int) i;
1412          break;
1413        }
1414    }
1415
1416  if (verboseLevel > 3)
1417    {
1418      G4cout << "SampleFinalStatePositron: sampled oscillator #" << targetOscillator << "." << G4endl;
1419      G4cout << "Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/eV
1420             << " eV " << G4endl; 
1421      G4cout << "Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/eV
1422             << " eV " << G4endl;
1423    }
1424
1425  //Constants
1426  G4double rb = kineticEnergy + 2.0*electron_mass_c2;
1427  G4double gam = 1.0+kineticEnergy/electron_mass_c2;
1428  G4double gam2 = gam*gam;
1429  G4double beta2 = (gam2-1.0)/gam2;
1430  G4double g12 = (gam+1.0)*(gam+1.0);
1431  G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam);
1432  //Bhabha coefficients
1433  G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0);
1434  G4double bha2 = amol*(3.0+1.0/g12);
1435  G4double bha3 = amol*2.0*gam*(gam-1.0)/g12;
1436  G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12;
1437
1438  //
1439  //Partial cross section of the active oscillator
1440  //
1441  G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
1442  G4double invResEne = 1.0/resEne;
1443  G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
1444  G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
1445
1446  G4double XHDL = 0.;
1447  G4double XHDT = 0.;
1448  G4double QM = 0.;
1449  G4double cps = 0.;
1450  G4double cp = 0.;
1451
1452  //Distant excitations XS (same as for electrons)
1453  if (resEne > cutEnergy && resEne < kineticEnergy)
1454    {
1455      cps = kineticEnergy*rb;
1456      cp = std::sqrt(cps);
1457      G4double XHDT0 = std::max(std::log(gam2)-beta2-delta,0.);
1458      if (resEne > 1.0e-6*kineticEnergy)
1459        {
1460          G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
1461          QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
1462        }
1463      else
1464        {
1465          QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
1466          QM *= (1.0-QM*0.5/electron_mass_c2);
1467        }
1468      if (QM < cutoffEne)
1469        {
1470          XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
1471            *invResEne;
1472          XHDT = XHDT0*invResEne;         
1473        }
1474      else
1475        {
1476          QM = cutoffEne;
1477          XHDL = 0.;
1478          XHDT = 0.;
1479        }
1480    }
1481  else
1482    {
1483      QM = cutoffEne;
1484      cps = 0.;
1485      cp = 0.;
1486      XHDL = 0.;
1487      XHDT = 0.;
1488    }
1489  //Close collisions (Bhabha)
1490  G4double wmaxc = kineticEnergy;
1491  G4double wcl = std::max(cutEnergy,cutoffEne);
1492  G4double rcl = wcl/kineticEnergy;
1493  G4double XHC = 0.;
1494  if (wcl < wmaxc)
1495    {
1496      G4double rl1 = 1.0-rcl;
1497      XHC = ((1.0/rcl-1.0)+bha1*std::log(rcl)+bha2*rl1
1498             + (bha3/2.0)*(rcl*rcl-1.0) 
1499             + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineticEnergy;
1500    }
1501
1502  //Total cross section per molecule for the active shell, in cm2
1503  G4double XHTOT = XHC + XHDL + XHDT;
1504
1505  //very small cross section, do nothing
1506  if (XHTOT < 1.e-14*barn)
1507    {
1508      kineticEnergy1=kineticEnergy;
1509      cosThetaPrimary=1.0;
1510      energySecondary=0.0;
1511      cosThetaSecondary=1.0;
1512      targetOscillator = numberOfOscillators-1;
1513      return;
1514    }
1515
1516  //decide which kind of interaction we'll have
1517  TST = XHTOT*G4UniformRand();
1518 
1519  // Hard close collision
1520  G4double TS1 = XHC;
1521  if (TST < TS1)
1522    {
1523      G4double rl1 = 1.0-rcl;
1524      G4double rk=0.;
1525      G4bool loopAgain = false;
1526      do
1527        {
1528          loopAgain = false;
1529          rk = rcl/(1.0-G4UniformRand()*rl1);
1530          G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
1531          if (G4UniformRand() > phi)
1532            loopAgain = true;
1533        }while(loopAgain);
1534      //energy and scattering angle (primary electron)
1535      G4double deltaE = rk*kineticEnergy;     
1536      kineticEnergy1 = kineticEnergy - deltaE;
1537      cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
1538      //energy and scattering angle of the delta ray
1539      energySecondary = deltaE - ionEne; //subtract ionisation energy
1540      cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
1541      if (verboseLevel > 3)
1542        G4cout << "SampleFinalStatePositron: sampled close collision " << G4endl;
1543      return;                               
1544    }
1545
1546  //Hard distant longitudinal collisions
1547  TS1 += XHDL;
1548  G4double deltaE = resEne;
1549  kineticEnergy1 = kineticEnergy - deltaE;
1550  if (TST < TS1)
1551    {
1552      G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
1553      G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand())
1554                       - (QS*0.5/electron_mass_c2));
1555      G4double QTREV = Q*(Q+2.0*electron_mass_c2);
1556      G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
1557      cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
1558      if (cosThetaPrimary > 1.) 
1559        cosThetaPrimary = 1.0;
1560      //energy and emission angle of the delta ray
1561      energySecondary = deltaE - ionEne;
1562      cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
1563      if (cosThetaSecondary > 1.0)
1564        cosThetaSecondary = 1.0;
1565      if (verboseLevel > 3)
1566        G4cout << "SampleFinalStatePositron: sampled distant longitudinal collision " << G4endl;
1567      return;     
1568    }
1569
1570  //Hard distant transverse collisions
1571  cosThetaPrimary = 1.0;
1572  //energy and emission angle of the delta ray
1573  energySecondary = deltaE - ionEne;
1574  cosThetaSecondary = 0.5;
1575
1576  if (verboseLevel > 3)   
1577    G4cout << "SampleFinalStatePositron: sampled distant transverse collision " << G4endl;
1578   
1579  return;
1580}
1581
Note: See TracBrowser for help on using the repository browser.