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

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

geant4 tag 9.4
File size: 62.3 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: G4PenelopeIonisationModel.cc,v 1.18 2010/12/01 15:20:35 pandola Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 26 Nov 2008   L Pandola    Migration from process to model
34// 17 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
35//                  - apply internal high-energy limit only in constructor
36//                  - do not apply low-energy limit (default is 0)
37//                  - added MinEnergyCut method
38//                  - do not change track status
39// 19 May 2009   L Pandola    Explicitely set to zero pointers deleted in
40//                            Initialise(), since they might be checked later on
41// 21 Oct 2009   L Pandola    Remove un-necessary fUseAtomicDeexcitation flag - now managed by
42//                            G4VEmModel::DeexcitationFlag()
43//                            Add ActivateAuger() method
44// 15 Mar 2010   L Pandola    Explicitely initialize Auger to false
45// 29 Mar 2010   L Pandola    Added a dummy ComputeCrossSectionPerAtom() method issueing a
46//                            warning if users try to access atomic cross sections via
47//                            G4EmCalculator
48// 15 Apr 2010   L. Pandola   Implemented model's own version of MinEnergyCut()
49// 23 Apr 2010   L. Pandola   Removed InitialiseElementSelectors() call. Useless here and
50//                            triggers fake warning messages
51//
52
53#include "G4PenelopeIonisationModel.hh"
54#include "G4ParticleDefinition.hh"
55#include "G4MaterialCutsCouple.hh"
56#include "G4ProductionCutsTable.hh"
57#include "G4DynamicParticle.hh"
58#include "G4Element.hh"
59#include "G4AtomicTransitionManager.hh"
60#include "G4AtomicDeexcitation.hh"
61#include "G4AtomicShell.hh"
62#include "G4Gamma.hh"
63#include "G4Electron.hh"
64#include "G4Positron.hh"
65#include "G4CrossSectionHandler.hh"
66#include "G4AtomicDeexcitation.hh"
67#include "G4VEMDataSet.hh"
68
69//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70
71
72G4PenelopeIonisationModel::G4PenelopeIonisationModel(const G4ParticleDefinition*,
73                                             const G4String& nam)
74  :G4VEmModel(nam),isInitialised(false),
75   kineticEnergy1(0.*eV),
76   cosThetaPrimary(1.0),
77   energySecondary(0.*eV),
78   cosThetaSecondary(0.0),
79   iOsc(-1),
80   crossSectionHandler(0),ionizationEnergy(0),
81   resonanceEnergy(0),occupationNumber(0),shellFlag(0),
82   theXSTable(0)
83{
84  fIntrinsicLowEnergyLimit = 100.0*eV;
85  fIntrinsicHighEnergyLimit = 100.0*GeV;
86  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
87  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
88  //
89  //
90  verboseLevel= 0;
91 
92  // Verbosity scale:
93  // 0 = nothing
94  // 1 = warning for energy non-conservation
95  // 2 = details of energy budget
96  // 3 = calculation of cross sections, file openings, sampling of atoms
97  // 4 = entering in methods
98
99  // Atomic deexcitation model activated by default
100  SetDeexcitationFlag(true);
101  ActivateAuger(false);
102 
103  //These vectors do not change when materials or cut change.
104  //Therefore I can read it at the constructor
105  ionizationEnergy = new std::map<G4int,G4DataVector*>;
106  resonanceEnergy  = new std::map<G4int,G4DataVector*>;
107  occupationNumber = new std::map<G4int,G4DataVector*>;
108  shellFlag = new std::map<G4int,G4DataVector*>;
109
110  ReadData(); //Read data from file
111
112}
113
114//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115
116G4PenelopeIonisationModel::~G4PenelopeIonisationModel()
117{
118 
119  if (crossSectionHandler) delete crossSectionHandler;
120
121  if (theXSTable)
122    {
123      for (size_t i=0; i<theXSTable->size(); i++)
124        delete (*theXSTable)[i];
125      delete theXSTable;
126     }
127
128 
129  std::map <G4int,G4DataVector*>::iterator i;
130  if (ionizationEnergy)
131    {
132      for (i=ionizationEnergy->begin();i != ionizationEnergy->end();i++)
133        if (i->second) delete i->second;
134      delete ionizationEnergy;
135    }
136  if (resonanceEnergy)
137    {
138      for (i=resonanceEnergy->begin();i != resonanceEnergy->end();i++)
139        if (i->second) delete i->second;
140      delete resonanceEnergy;
141    }
142  if (occupationNumber)
143    {
144      for (i=occupationNumber->begin();i != occupationNumber->end();i++)
145        if (i->second) delete i->second;
146      delete occupationNumber;
147    }
148  if (shellFlag)
149    {
150      for (i=shellFlag->begin();i != shellFlag->end();i++)
151        if (i->second) delete i->second;
152      delete shellFlag;
153    }
154}
155
156//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157
158void G4PenelopeIonisationModel::Initialise(const G4ParticleDefinition* particle,
159                                       const G4DataVector& )
160{
161  if (verboseLevel > 3)
162    G4cout << "Calling G4PenelopeIonisationModel::Initialise()" << G4endl;
163
164  //Delete and re-initialize the cross section handler
165  if (crossSectionHandler)
166    {
167      crossSectionHandler->Clear();
168      delete crossSectionHandler;
169      crossSectionHandler = 0;
170    }
171
172   if (theXSTable)
173    {
174      for (size_t i=0; i<theXSTable->size(); i++)
175        {
176          delete (*theXSTable)[i];
177          (*theXSTable)[i] = 0;
178        }
179      delete theXSTable;
180      theXSTable = 0;
181     }
182
183  crossSectionHandler = new G4CrossSectionHandler();
184  crossSectionHandler->Clear();
185  G4String crossSectionFile = "NULL";
186 
187  if (particle == G4Electron::Electron())
188     crossSectionFile = "penelope/ion-cs-el-";
189  else if (particle == G4Positron::Positron())
190     crossSectionFile = "penelope/ion-cs-po-"; 
191  crossSectionHandler->LoadData(crossSectionFile);
192  //This is used to retrieve cross section values later on
193  G4VEMDataSet* emdata = 
194    crossSectionHandler->BuildMeanFreePathForMaterials();
195  //The method BuildMeanFreePathForMaterials() is required here only to force
196  //the building of an internal table: the output pointer can be deleted
197  delete emdata;
198
199  if (verboseLevel > 2) 
200    G4cout << "Loaded cross section files for PenelopeIonisationModel" << G4endl;
201
202  if (verboseLevel > 2) {
203    G4cout << "Penelope Ionisation model is initialized " << G4endl
204           << "Energy range: "
205           << LowEnergyLimit() / keV << " keV - "
206           << HighEnergyLimit() / GeV << " GeV"
207           << G4endl;
208  }
209
210  if(isInitialised) return;
211  fParticleChange = GetParticleChangeForLoss();
212  isInitialised = true; 
213}
214
215//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
216
217G4double G4PenelopeIonisationModel::CrossSectionPerVolume(const G4Material* material,
218                                           const G4ParticleDefinition* theParticle,
219                                           G4double energy,
220                                           G4double cutEnergy,
221                                           G4double)           
222{ 
223  // Penelope model to calculate the cross section for inelastic collisions above the
224  // threshold. It makes use of the Generalised Oscillator Strength (GOS) model from
225  //  D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
226  //
227  // The total cross section (hard+soft) is read from a database file (element per
228  // element), while the ratio hard-to-total is calculated analytically by taking
229  // into account the atomic oscillators coming into the play for a given threshold.
230  // This is done by the method CalculateCrossSectionsRatio().
231  // For incident e- the maximum energy allowed for the delta-rays is energy/2.
232  // because particles are undistinghishable.
233  //
234  // The contribution is splitted in three parts: distant longitudinal collisions,
235  // distant transverse collisions and close collisions. Each term is described by
236  // its own analytical function.
237  // Fermi density correction is calculated analytically according to
238  //  U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
239  //
240  if (verboseLevel > 3)
241    G4cout << "Calling CrossSectionPerVolume() of G4PenelopeIonisationModel" << G4endl;
242
243  SetupForMaterial(theParticle, material, energy);
244
245  // VI - should be check at initialisation not in run time
246  /*
247  if (!crossSectionHandler)
248    {
249      G4cout << "G4PenelopeIonisationModel::CrossSectionPerVolume" << G4endl;
250      G4cout << "The cross section handler is not correctly initialized" << G4endl;
251      G4Exception();
252    }
253  */ 
254  if (!theXSTable)
255    {
256      if (verboseLevel > 2)
257        {
258          G4cout << "G4PenelopeIonisationModel::CrossSectionPerVolume" << G4endl;
259          G4cout << "Going to build Cross Section table " << G4endl;
260        }
261      theXSTable = new std::vector<G4VEMDataSet*>;
262      theXSTable = BuildCrossSectionTable(theParticle); 
263    }
264 
265  G4double totalCross = 0.0; 
266  G4double cross = 0.0;
267  const G4ElementVector* theElementVector = material->GetElementVector();
268  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
269  G4double electronVolumeDensity =
270    material->GetTotNbOfElectPerVolume();  //electron density
271
272  if (verboseLevel > 4)
273    G4cout << "Electron volume density of " << material->GetName() << ": " << 
274      electronVolumeDensity*cm3 << " electrons/cm3" << G4endl;
275
276  G4int nelm = material->GetNumberOfElements();
277  for (G4int i=0; i<nelm; i++) 
278    {
279      G4int iZ = (G4int) (*theElementVector)[i]->GetZ();
280      G4double ratio = CalculateCrossSectionsRatio(energy,cutEnergy,iZ,electronVolumeDensity,
281                                                   theParticle);
282      cross += theAtomNumDensityVector[i]*
283        crossSectionHandler->FindValue(iZ,energy)*ratio;
284      totalCross += theAtomNumDensityVector[i]*
285        crossSectionHandler->FindValue(iZ,energy);
286    }
287
288  if (verboseLevel > 2)
289  {
290    G4cout << "G4PenelopeIonisationModel " << G4endl;
291    G4cout << "Mean free path for delta emission > " << cutEnergy/keV << " keV at " << 
292      energy/keV << " keV = " << (1./cross)/mm << " mm" << G4endl;
293    G4cout << "Total free path for ionisation (no threshold) at " << 
294      energy/keV << " keV = " << (1./totalCross)/mm << " mm" << G4endl;
295  }
296  return cross;
297}
298
299//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
300
301//This is a dummy method. Never inkoved by the tracking, it just issues
302//a warning if one tries to get Cross Sections per Atom via the
303//G4EmCalculator.
304G4double G4PenelopeIonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
305                                                             G4double,
306                                                             G4double,
307                                                             G4double,
308                                                             G4double,
309                                                             G4double)
310{
311  G4cout << "*** G4PenelopeIonisationModel -- WARNING ***" << G4endl;
312  G4cout << "Penelope Ionisation model does not calculate cross section _per atom_ " << G4endl;
313  G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
314  G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
315  return 0;
316}
317
318//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
319
320G4double G4PenelopeIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
321                                           const G4ParticleDefinition* theParticle,
322                                           G4double kineticEnergy,
323                                           G4double cutEnergy)
324{
325  // Penelope model to calculate the stopping power for soft inelastic collisions
326  // below the threshold. It makes use of the Generalised Oscillator Strength (GOS)
327  // model from
328  //  D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
329  //
330  // The stopping power is calculated analytically using the dsigma/dW cross
331  // section from the GOS models, which includes separate contributions from
332  // distant longitudinal collisions, distant transverse collisions and
333  // close collisions. Only the atomic oscillators that come in the play
334  // (according to the threshold) are considered for the calculation.
335  // Differential cross sections have a different form for e+ and e-.
336  //
337  // Fermi density correction is calculated analytically according to
338  //  U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
339
340  if (verboseLevel > 3)
341    G4cout << "Calling ComputeDEDX() of G4PenelopeIonisationModel" << G4endl;
342
343  G4double sPower = 0.0;
344  const G4ElementVector* theElementVector = material->GetElementVector();
345  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
346  G4double electronVolumeDensity =
347    material->GetTotNbOfElectPerVolume();  //electron density
348
349  //Loop on the elements in the material
350  G4int nelm = material->GetNumberOfElements();
351  for (G4int i=0; i<nelm; i++) 
352    {
353      G4int iZ = (G4int) (*theElementVector)[i]->GetZ();
354   
355      //Calculate stopping power contribution from each element
356      //Constants
357      G4double gamma = 1.0+kineticEnergy/electron_mass_c2;
358      G4double gamma2 = gamma*gamma;
359      G4double beta2 = (gamma2-1.0)/gamma2;
360      G4double constant = pi*classic_electr_radius*classic_electr_radius
361        *2.0*electron_mass_c2/beta2;
362 
363      G4double delta = CalculateDeltaFermi(kineticEnergy,iZ,electronVolumeDensity);
364      G4int nbOsc = (G4int) resonanceEnergy->find(iZ)->second->size();
365      G4double stoppingPowerForElement = 0.0;
366      //Loop on oscillators of element Z
367      for (G4int iosc=0;iosc<nbOsc;iosc++)
368        {
369          G4double S1 = 0.0;
370          G4double resEnergy = (*(resonanceEnergy->find(iZ)->second))[iosc];
371          if (theParticle == G4Electron::Electron())
372            S1 = ComputeStoppingPowerForElectrons(kineticEnergy,cutEnergy,delta,resEnergy);
373          else if (theParticle == G4Positron::Positron())
374            S1 = ComputeStoppingPowerForPositrons(kineticEnergy,cutEnergy,delta,resEnergy);
375          G4double occupNb = (*(occupationNumber->find(iZ)->second))[iosc];
376          stoppingPowerForElement += occupNb*constant*S1;
377        }
378
379      sPower += stoppingPowerForElement*theAtomNumDensityVector[i];
380    }
381  if (verboseLevel > 2)
382    {
383      G4cout << "G4PenelopeIonisationModel " << G4endl;
384      G4cout << "Stopping power < " << cutEnergy/keV << " keV at " << 
385        kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl;
386    } 
387  return sPower;
388}
389
390//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
391
392G4double G4PenelopeIonisationModel::MinEnergyCut(const G4ParticleDefinition*,
393                                                 const G4MaterialCutsCouple*)
394{
395  return fIntrinsicLowEnergyLimit;
396}
397
398//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
399
400void G4PenelopeIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
401                                              const G4MaterialCutsCouple* couple,
402                                              const G4DynamicParticle* aDynamicParticle,
403                                              G4double cutE, G4double)
404{
405  // Penelope model to sample the final state following an hard inelastic interaction.
406  // It makes use of the Generalised Oscillator Strength (GOS) model from
407  //  D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
408  //
409  // The GOS model is used to calculate the individual cross sections for all
410  // the atomic oscillators coming in the play, taking into account the three
411  // contributions (distant longitudinal collisions, distant transverse collisions and
412  // close collisions). Then the target shell and the interaction channel are
413  // sampled. Final state of the delta-ray (energy, angle) are generated according
414  // to the analytical distributions (dSigma/dW) for the selected interaction
415  // channels.
416  // For e-, the maximum energy for the delta-ray is initialEnergy/2. (because
417  // particles are indistinghusbable), while it is the full initialEnergy for
418  // e+.
419  // The efficiency of the random sampling algorithm (e.g. for close collisions)
420  // decreases when initial and cutoff energy increase (e.g. 87% for 10-keV primary
421  // and 1 keV threshold, 99% for 10-MeV primary and 10-keV threshold).
422  // Differential cross sections have a different form for e+ and e-.
423  //
424  // WARNING: The model provides an _average_ description of inelastic collisions.
425  // Anyway, the energy spectrum associated to distant excitations of a given
426  // atomic shell is approximated as a single resonance. The simulated energy spectra
427  // show _unphysical_ narrow peaks at energies that are multiple of the shell
428  // resonance energies. The spurious speaks are automatically smoothed out after
429  // multiple inelastic collisions.
430  //
431  // The model determines also the original shell from which the delta-ray is expelled,
432  // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
433  //
434  // Fermi density correction is calculated analytically according to
435  //  U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
436
437  if (verboseLevel > 3)
438    G4cout << "Calling SamplingSecondaries() of G4PenelopeIonisationModel" << G4endl;
439
440  G4double kineticEnergy0 = aDynamicParticle->GetKineticEnergy();
441  const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition();
442
443  if (kineticEnergy0 <= fIntrinsicLowEnergyLimit)
444  {
445    fParticleChange->SetProposedKineticEnergy(0.);
446    fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy0);
447    return ;
448  }
449  const G4double electronVolumeDensity = 
450    couple->GetMaterial()->GetTotNbOfElectPerVolume(); 
451  G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
452
453  //Initialise final-state variables. The proper values will be set by the methods
454  // CalculateDiscreteForElectrons() and CalculateDiscreteForPositrons()
455  kineticEnergy1=kineticEnergy0;
456  cosThetaPrimary=1.0;
457  energySecondary=0.0;
458  cosThetaSecondary=1.0;
459
460  // Select randomly one element in the current material
461  if (verboseLevel > 2)
462    G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
463
464  //Use sampler of G4CrossSectionHandler for now
465  // atom can be selected effitiantly if element selectors are initialised   
466  //const G4Element* anElement = SelectRandomAtom(couple,theParticle,kineticEnergy0);
467
468  G4int iZ = SampleRandomAtom(couple,kineticEnergy0);
469 
470  const G4ProductionCutsTable* theCoupleTable=
471    G4ProductionCutsTable::GetProductionCutsTable();
472  size_t indx = couple->GetIndex();
473  G4double cutG = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
474 
475  if (verboseLevel > 2)
476    G4cout << "Selected Z = " << iZ << G4endl;
477 
478  // The method CalculateDiscreteForXXXX() set the private variables:
479  // kineticEnergy1 = energy of the primary electron/positron after the interaction
480  // cosThetaPrimary = std::cos(theta) of the primary after the interaction
481  // energySecondary = energy of the secondary electron
482  // cosThetaSecondary = std::cos(theta) of the secondary
483  if (theParticle == G4Electron::Electron())
484    CalculateDiscreteForElectrons(kineticEnergy0,cutE,iZ,electronVolumeDensity);
485  else if (theParticle == G4Positron::Positron())
486    CalculateDiscreteForPositrons(kineticEnergy0,cutE,iZ,electronVolumeDensity);
487
488  // if (energySecondary == 0) return;
489
490  //Update the primary particle
491  G4double sint = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary);
492  G4double phi  = twopi * G4UniformRand();
493  G4double dirx = sint * std::cos(phi);
494  G4double diry = sint * std::sin(phi);
495  G4double dirz = cosThetaPrimary;
496
497  G4ThreeVector electronDirection1(dirx,diry,dirz);
498  electronDirection1.rotateUz(particleDirection0);
499 
500  if (kineticEnergy1 > 0)
501    {
502      fParticleChange->ProposeMomentumDirection(electronDirection1);
503      fParticleChange->SetProposedKineticEnergy(kineticEnergy1);
504    }
505  else
506    {
507      fParticleChange->SetProposedKineticEnergy(0.);
508    }
509 
510  //Generate the delta ray
511  G4int iosc2 = 0;
512  G4double ioniEnergy = 0.0;
513  if (iOsc > 0) {
514    ioniEnergy=(*(ionizationEnergy->find(iZ)->second))[iOsc];
515    iosc2 = (ionizationEnergy->find(iZ)->second->size()) - iOsc; //they are in reversed order
516  }
517
518  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
519  G4double bindingEnergy = 0.0;
520  G4int shellId = 0;
521  if (iOsc > 0)
522    {
523      const G4AtomicShell* shell = transitionManager->Shell(iZ,iosc2-1); // Modified by Alf
524      bindingEnergy = shell->BindingEnergy();
525      shellId = shell->ShellId();
526    }
527
528  G4double ionEnergy = bindingEnergy; //energy spent to ionise the atom according to G4dabatase
529  G4double eKineticEnergy = energySecondary;
530
531  //This is an awful thing: Penelope generates the fluorescence only for L and K shells
532  //(i.e. Osc = 1 --> 4). For high-Z, the other shells can be quite relevant. In this case
533  //one MUST ensure ''by hand'' the energy conservation. Then there is the other problem that
534  //the fluorescence database of Penelope doesn not match that of Geant4.
535
536  G4double energyBalance = kineticEnergy0 - kineticEnergy1 - energySecondary; //Penelope Balance
537
538  if (std::abs(energyBalance) < 1*eV)
539    //in this case Penelope didn't subtract the fluorescence energy: do here by hand
540    eKineticEnergy = energySecondary - bindingEnergy;   
541  else 
542    //Penelope subtracted the fluorescence, but one has to match the databases
543    eKineticEnergy = energySecondary+ioniEnergy-bindingEnergy;
544 
545  G4double localEnergyDeposit = ionEnergy; 
546  G4double energyInFluorescence = 0.0*eV;
547
548  if(DeexcitationFlag() && iZ > 5) 
549    {
550      if (ionEnergy > cutG || ionEnergy > cutE)
551        {
552          deexcitationManager.SetCutForSecondaryPhotons(cutG);
553          deexcitationManager.SetCutForAugerElectrons(cutE);
554          std::vector<G4DynamicParticle*> *photonVector =
555            deexcitationManager.GenerateParticles(iZ,shellId);
556          //Check for secondaries
557          if(photonVector) 
558            {
559              for (size_t k=0;k<photonVector->size();k++)
560                {
561                  G4DynamicParticle* aPhoton = (*photonVector)[k];
562                  if (aPhoton)
563                    {
564                      G4double itsEnergy = aPhoton->GetKineticEnergy();
565                      if (itsEnergy <= localEnergyDeposit)
566                        {
567                          if(aPhoton->GetDefinition() == G4Gamma::Gamma())
568                            energyInFluorescence += itsEnergy;
569                          localEnergyDeposit -= itsEnergy;
570                          fvect->push_back(aPhoton);
571                        }
572                      else
573                        {
574                          delete aPhoton;
575                          (*photonVector)[k] = 0;
576                        }
577                    }
578                }
579              delete photonVector;
580            }
581        }
582    }
583
584  // Generate the delta ray
585  G4double sin2 = std::sqrt(1. - cosThetaSecondary*cosThetaSecondary);
586  G4double phi2  = twopi * G4UniformRand();
587 
588  G4double xEl = sin2 * std::cos(phi2); 
589  G4double yEl = sin2 * std::sin(phi2);
590  G4double zEl = cosThetaSecondary;
591  G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
592  eDirection.rotateUz(particleDirection0);
593
594  G4DynamicParticle* deltaElectron = new G4DynamicParticle (G4Electron::Electron(),
595                                                            eDirection,eKineticEnergy) ;
596  fvect->push_back(deltaElectron);
597 
598  if (localEnergyDeposit < 0)
599    {
600      G4cout << "WARNING-" 
601             << "G4PenelopeIonisationModel::SampleSecondaries - Negative energy deposit"
602             << G4endl;
603      localEnergyDeposit=0.;
604    }
605  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
606 
607  if (verboseLevel > 1)
608    {
609      G4cout << "-----------------------------------------------------------" << G4endl;
610      G4cout << "Energy balance from G4PenelopeIonisation" << G4endl;
611      G4cout << "Incoming primary energy: " << kineticEnergy0/keV << " keV" << G4endl;
612      G4cout << "-----------------------------------------------------------" << G4endl;
613      G4cout << "Outgoing primary energy: " << kineticEnergy1/keV << " keV" << G4endl;
614      G4cout << "Delta ray " << eKineticEnergy/keV << " keV" << G4endl;
615      G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
616      G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
617      G4cout << "Total final state: " << (eKineticEnergy+energyInFluorescence+kineticEnergy1+
618                                          localEnergyDeposit)/keV << 
619        " keV" << G4endl;
620      G4cout << "-----------------------------------------------------------" << G4endl;
621    }
622  if (verboseLevel > 0)
623    {
624      G4double energyDiff = std::fabs(eKineticEnergy+energyInFluorescence+kineticEnergy1+
625                                      localEnergyDeposit-kineticEnergy0);
626      if (energyDiff > 0.05*keV)
627        G4cout << "Warning from G4PenelopeIonisation: problem with energy conservation: " << 
628          (eKineticEnergy+energyInFluorescence+kineticEnergy1+localEnergyDeposit)/keV << 
629          " keV (final) vs. " << 
630          kineticEnergy0/keV << " keV (initial)" << G4endl;
631    }
632}
633
634//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
635
636void G4PenelopeIonisationModel::ReadData()
637{
638  if (verboseLevel > 2)
639    {
640      G4cout << "Data from G4PenelopeIonisationModel read " << G4endl;
641    }
642  char* path = getenv("G4LEDATA");
643  if (!path)
644    {
645      G4String excep = "G4PenelopeIonisationModel - G4LEDATA environment variable not set!";
646      G4Exception(excep);
647      return;
648    }
649  G4String pathString(path);
650  G4String pathFile = pathString + "/penelope/ion-pen.dat";
651  std::ifstream file(pathFile);
652 
653  if (!file.is_open())
654    {
655      G4String excep = "G4PenelopeIonisationModel - data file " + pathFile + " not found!";
656      G4Exception(excep);
657    }
658
659 if (!ionizationEnergy || !resonanceEnergy || !occupationNumber || !shellFlag)
660    {
661      G4String excep = "G4PenelopeIonisationModel: problem with reading data from file";
662      G4Exception(excep);
663      return;
664    }
665
666 G4int Z=1,nLevels=0;
667 G4int test,test1;
668
669 do{
670   file >> Z >> nLevels;
671   //Check for nLevels validity, before using it in a loop
672   if (nLevels<0 || nLevels>64)
673     {
674       G4String excep = "G4PenelopeIonisationModel: corrupted data file ?";
675       G4Exception(excep);
676       return;
677     }
678   //Allocate space for storage
679   G4DataVector* occVector = new G4DataVector;
680   G4DataVector* ionEVector = new G4DataVector;
681   G4DataVector* resEVector = new G4DataVector;
682   G4DataVector* shellIndVector = new G4DataVector;
683   //
684   G4double a1,a2,a3,a4;
685   G4int k1,k2,k3;
686   for (G4int h=0;h<nLevels;h++)
687     {
688       //index,occup number,ion energy,res energy,fj0,kz,shell flag
689       file >> k1 >> a1 >> a2 >> a3 >> a4 >> k2 >> k3;
690       //Make explicit unit of measurements for ionisation and resonance energy,
691       // which is MeV
692       a2 *= MeV;
693       a3 *= MeV;
694       //
695       occVector->push_back(a1);
696       ionEVector->push_back(a2);
697       resEVector->push_back(a3);
698       shellIndVector->push_back((G4double) k3);
699     }
700   //Ok, done for element Z
701   occupationNumber->insert(std::make_pair(Z,occVector));
702   ionizationEnergy->insert(std::make_pair(Z,ionEVector));
703   resonanceEnergy->insert(std::make_pair(Z,resEVector));
704   shellFlag->insert(std::make_pair(Z,shellIndVector));
705   file >> test >> test1; //-1 -1 close the data for each Z
706   if (test > 0) {
707     G4String excep = "G4PenelopeIonisationModel - data file corrupted!";
708     G4Exception(excep);
709   }
710 }while (test != -2); //the very last Z is closed with -2 instead of -1
711}
712                                   
713                                           
714//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
715
716G4double G4PenelopeIonisationModel::CalculateDeltaFermi(G4double kinEnergy ,G4int Z,
717                                                        G4double electronVolumeDensity)
718{
719  G4double plasmaEnergyCoefficient = 1.377e-39*(MeV*MeV*m3); //(e*hbar)^2/(epsilon0*electron_mass)
720  G4double plasmaEnergySquared = plasmaEnergyCoefficient*electronVolumeDensity;
721  // std::sqrt(plasmaEnergySquared) is the plasma energy of the solid (MeV)
722  G4double gam = 1.0+kinEnergy/electron_mass_c2;
723  G4double gam2=gam*gam;
724  G4double delta = 0.0;
725
726  //Density effect
727  G4double TST = ((G4double) Z)/(gam2*plasmaEnergySquared); 
728 
729  G4double wl2 = 0.0;
730  G4double fdel=0.0;
731  G4double wr=0;
732  size_t nbOsc = resonanceEnergy->find(Z)->second->size();
733  for(size_t i=0;i<nbOsc;i++)
734    {
735      G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
736      wr = (*(resonanceEnergy->find(Z)->second))[i];
737      fdel += occupNb/(wr*wr+wl2);
738    }
739  if (fdel < TST) return delta;
740  G4double help1 = (*(resonanceEnergy->find(Z)->second))[nbOsc-1];
741  wl2 = help1*help1;
742  do{
743    wl2=wl2*2.0;
744    fdel = 0.0;
745    for (size_t ii=0;ii<nbOsc;ii++){
746      G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[ii];
747      wr = (*(resonanceEnergy->find(Z)->second))[ii];
748      fdel += occupNb/(wr*wr+wl2);
749    }
750  }while (fdel > TST);
751  G4double wl2l=0.0;
752  G4double wl2u = wl2;
753  G4double control = 0.0;
754  do{
755    wl2=0.5*(wl2l+wl2u);
756    fdel = 0.0;
757    for (size_t jj=0;jj<nbOsc;jj++){
758      G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[jj];
759      wr = (*(resonanceEnergy->find(Z)->second))[jj];
760      fdel += occupNb/(wr*wr+wl2);
761    }
762    if (fdel > TST)
763      wl2l = wl2;
764    else
765      wl2u = wl2;
766    control = wl2u-wl2l-wl2*1e-12; 
767  }while(control>0);
768
769  //Density correction effect
770   for (size_t kk=0;kk<nbOsc;kk++){
771      G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[kk];
772      wr = (*(resonanceEnergy->find(Z)->second))[kk];
773      delta += occupNb*std::log(1.0+wl2/(wr*wr));
774    }
775   delta = (delta/((G4double) Z))-wl2/(gam2*plasmaEnergySquared); 
776   return delta;
777}
778
779//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
780
781void 
782G4PenelopeIonisationModel::CalculateDiscreteForElectrons(G4double kinEnergy,
783                                                         G4double cutoffEnergy,
784                                                         G4int Z,
785                                                         G4double electronVolumeDensity)
786{
787  if (verboseLevel > 2)
788    G4cout << "Entering in CalculateDiscreteForElectrons() for energy " << 
789      kinEnergy/keV << " keV " << G4endl;
790
791  //Initialization of variables to be calculated in this routine
792  kineticEnergy1=kinEnergy;
793  cosThetaPrimary=1.0;
794  energySecondary=0.0;
795  cosThetaSecondary=1.0;
796  iOsc=-1;
797
798  //constants
799  G4double rb=kinEnergy+2.0*electron_mass_c2;
800  G4double gamma = 1.0+kinEnergy/electron_mass_c2;
801  G4double gamma2 = gamma*gamma;
802  G4double beta2 = (gamma2-1.0)/gamma2;
803  G4double amol = (gamma-1.0)*(gamma-1.0)/gamma2;
804  G4double cps = kinEnergy*rb;
805  G4double cp = std::sqrt(cps);
806 
807  G4double delta = CalculateDeltaFermi(kinEnergy,Z,electronVolumeDensity);
808  G4double distantTransvCS0 = std::max(std::log(gamma2)-beta2-delta,0.0);
809
810  G4double rl,rl1;
811
812  if (cutoffEnergy > kinEnergy) return; //delta rays are not generated
813
814  G4DataVector* qm = new G4DataVector();
815  G4DataVector* cumulHardCS = new G4DataVector();
816  std::vector<G4int> *typeOfInteraction = new std::vector<G4int>;
817  G4DataVector* nbOfLevel = new G4DataVector();
818 
819  //Hard close collisions with outer shells
820  G4double wmaxc = 0.5*kinEnergy;
821  G4double closeCS0 = 0.0;
822  G4double closeCS = 0.0;
823  if (cutoffEnergy>0.1*eV) 
824    {
825      rl=cutoffEnergy/kinEnergy;
826      rl1=1.0-rl;
827      if (rl < 0.5)
828        closeCS0 = (amol*(0.5-rl)+(1.0/rl)-(1.0/rl1)+(1.0-amol)*std::log(rl/rl1))/kinEnergy;
829    }
830
831  // Cross sections for the different oscillators
832
833  // totalHardCS contains the cumulative hard interaction cross section for the different
834  // excitable levels and the different interaction channels (close, distant, etc.),
835  // i.e.
836  // cumulHardCS[0] = 0.0
837  // cumulHardCS[1] = 1st excitable level (distant longitudinal only)
838  // cumulHardCS[2] = 1st excitable level (distant longitudinal + transverse)
839  // cumulHardCS[3] = 1st excitable level (distant longitudinal + transverse + close)
840  // cumulHardCS[4] = 1st excitable level (all channels) + 2nd excitable level (distant long only)
841  // etc.
842  // This is used for sampling the atomic level which is ionised and the channel of the
843  // interaction.
844  //
845  // For each index iFill of the cumulHardCS vector,
846  // nbOfLevel[iFill] contains the current excitable atomic level and
847  // typeOfInteraction[iFill] contains the current interaction channel, with the legenda:
848  //   1 = distant longitudinal interaction
849  //   2 = distant transverse interaction
850  //   3 = close collision
851  //   4 = close collision with outer shells (in this case nbOfLevel < 0 --> no binding energy)
852
853
854  G4int nOscil = ionizationEnergy->find(Z)->second->size();
855  G4double totalHardCS = 0.0;
856  G4double involvedElectrons = 0.0;
857  for (G4int i=0;i<nOscil;i++){
858    G4double wi = (*(resonanceEnergy->find(Z)->second))[i];
859    G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
860 
861    //Distant excitations
862    if (wi>cutoffEnergy && wi<kinEnergy)
863      {
864        if (wi>(1e-6*kinEnergy))
865          {
866            G4double cpp=std::sqrt((kinEnergy-wi)*(kinEnergy-wi+2.0*electron_mass_c2));
867            qm->push_back(std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2);
868          }
869        else
870          {
871            qm->push_back((wi*wi)/(beta2*2.0*electron_mass_c2));
872          }
873        if ((*qm)[i] < wi)
874          {
875           
876            G4double distantLongitCS =  occupNb*std::log(wi*((*qm)[i]+2.0*electron_mass_c2)/
877                                         ((*qm)[i]*(wi+2.0*electron_mass_c2)))/wi;
878            cumulHardCS->push_back(totalHardCS);
879            typeOfInteraction->push_back(1); //distant longitudinal
880            nbOfLevel->push_back((G4double) i); //only excitable level are counted
881            totalHardCS += distantLongitCS;
882           
883            G4double distantTransvCS = occupNb*distantTransvCS0/wi;
884           
885            cumulHardCS->push_back(totalHardCS);
886            typeOfInteraction->push_back(2); //distant tranverse
887            nbOfLevel->push_back((G4double) i);
888            totalHardCS += distantTransvCS;
889          }
890      }
891    else 
892      qm->push_back(wi);
893     
894    //close collisions
895    if(wi < wmaxc)
896      {
897        if (wi < cutoffEnergy) 
898          involvedElectrons += occupNb;
899        else
900          {
901            rl=wi/kinEnergy;
902            rl1=1.0-rl;
903            closeCS = occupNb*(amol*(0.5-rl)+(1.0/rl)-(1.0/rl1)+(1.0-amol)*std::log(rl/rl1))/kinEnergy;
904            cumulHardCS->push_back(totalHardCS);
905            typeOfInteraction->push_back(3); //close
906            nbOfLevel->push_back((G4double) i);
907            totalHardCS += closeCS;
908          }
909      }
910  } // loop on the levels
911 
912  cumulHardCS->push_back(totalHardCS);
913  typeOfInteraction->push_back(4); //close interaction with outer shells
914  nbOfLevel->push_back(-1.0);
915  totalHardCS += involvedElectrons*closeCS0;
916  cumulHardCS->push_back(totalHardCS); //this is the final value of the totalHardCS
917
918  if (totalHardCS < 1e-30) {
919    kineticEnergy1=kinEnergy;
920    cosThetaPrimary=1.0;
921    energySecondary=0.0;
922    cosThetaSecondary=0.0;
923    iOsc=-1;
924    delete qm;
925    delete cumulHardCS;
926    delete typeOfInteraction;
927    delete nbOfLevel;
928    return;
929  }
930
931  //Testing purposes
932  if (verboseLevel > 6)
933    {
934      for (size_t t=0;t<cumulHardCS->size();t++)
935        G4cout << (*cumulHardCS)[t] << " " << (*typeOfInteraction)[t] << 
936          " " << (*nbOfLevel)[t] << G4endl;
937    }
938
939  //Selection of the active oscillator on the basis of the cumulative cross sections
940  G4double TST = totalHardCS*G4UniformRand();
941  G4int is=0;
942  G4int js= nbOfLevel->size();
943  do{
944    G4int it=(is+js)/2;
945    if (TST > (*cumulHardCS)[it]) is=it;
946    if (TST <= (*cumulHardCS)[it]) js=it;
947  }while((js-is) > 1);
948
949  G4double UII=0.0;
950  G4double rkc=cutoffEnergy/kinEnergy;
951  G4double dde;
952  G4int kks;
953
954  G4int sampledInteraction = (*typeOfInteraction)[is];
955  iOsc = (G4int) (*nbOfLevel)[is];
956
957  if (verboseLevel > 4) 
958    G4cout << "Chosen interaction #:" << sampledInteraction << " on level " << iOsc << G4endl;
959
960  //Generates the final state according to the sampled level and
961  //interaction channel
962 
963  if (sampledInteraction == 1)  //Hard distant longitudinal collisions
964    {
965      dde= (*(resonanceEnergy->find(Z)->second))[iOsc];
966      kineticEnergy1=kinEnergy-dde;
967      G4double qs=(*qm)[iOsc]/(1.0+((*qm)[iOsc]/(2.0*electron_mass_c2)));
968      G4double q=qs/(std::pow((qs/dde)*(1.0+(0.5*dde/electron_mass_c2)),G4UniformRand())-(0.5*qs/electron_mass_c2));
969      G4double qtrev = q*(q+2.0*electron_mass_c2);
970      G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
971      cosThetaPrimary = (cpps+cps-qtrev)/(2.0*cp*std::sqrt(cpps));
972      if (cosThetaPrimary>1.0) cosThetaPrimary=1.0;
973      //Energy and emission angle of the delta ray
974      kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
975      //kks > 4 means that we are in an outer shell
976      if (kks>4) 
977        energySecondary=dde;
978      else
979        energySecondary=dde-(*(ionizationEnergy->find(Z)->second))[iOsc];
980       
981      cosThetaSecondary = 0.5*(dde*(kinEnergy+rb-dde)+qtrev)/std::sqrt(cps*qtrev);
982      if (cosThetaSecondary>1.0) cosThetaSecondary=1.0;
983    }
984  else if (sampledInteraction == 2)  //Hard distant transverse collisions
985    {
986      dde=(*(resonanceEnergy->find(Z)->second))[iOsc];
987      kineticEnergy1=kinEnergy-dde;
988      cosThetaPrimary=1.0;
989      //Energy and emission angle of the delta ray
990      kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
991      if (kks>4)
992        {
993          energySecondary=dde;
994        }
995      else
996        {
997          energySecondary=dde-(*(ionizationEnergy->find(Z)->second))[iOsc];
998        }
999      cosThetaSecondary = 1.0;
1000    }
1001  else if (sampledInteraction == 3 || sampledInteraction == 4) //Close interaction
1002    {
1003      if (sampledInteraction == 4) //interaction with inner shells
1004        {
1005          UII=0.0;
1006          rkc = cutoffEnergy/kinEnergy;
1007          iOsc = -1;
1008        }
1009      else
1010        {
1011          kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
1012          if (kks > 4) {
1013            UII=0.0;
1014          }
1015          else
1016            {
1017              UII = (*(ionizationEnergy->find(Z)->second))[iOsc];
1018            }
1019          rkc = (*(resonanceEnergy->find(Z)->second))[iOsc]/kinEnergy;
1020        }
1021    G4double A = 0.5*amol;
1022    G4double arkc = A*0.5*rkc;
1023    G4double phi,rk2,rk,rkf;
1024    do{
1025      G4double fb = (1.0+arkc)*G4UniformRand();
1026      if (fb<1.0)
1027        {
1028          rk=rkc/(1.0-fb*(1.0-(rkc*2.0)));
1029        }
1030      else{
1031        rk = rkc+(fb-1.0)*(0.5-rkc)/arkc;
1032      }
1033      rk2 = rk*rk;
1034      rkf = rk/(1.0-rk);
1035      phi = 1.0+(rkf*rkf)-rkf+amol*(rk2+rkf);
1036    }while ((G4UniformRand()*(1.0+A*rk2)) > phi);
1037    //Energy and scattering angle (primary electron);
1038    kineticEnergy1 = kinEnergy*(1.0-rk);
1039    cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kinEnergy*(rb-(rk*kinEnergy))));
1040    //Energy and scattering angle of the delta ray
1041    energySecondary = kinEnergy-kineticEnergy1-UII;
1042    cosThetaSecondary = std::sqrt(rk*kinEnergy*rb/(kinEnergy*(rk*kinEnergy+2.0*electron_mass_c2)));
1043    }
1044  else
1045    {
1046      G4String excep = "G4PenelopeIonisationModel - Error in the calculation of the final state";
1047      G4Exception(excep);
1048    }
1049
1050  delete qm;
1051  delete cumulHardCS;
1052 
1053  delete typeOfInteraction;
1054  delete nbOfLevel;
1055
1056  return;
1057}
1058
1059//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1060
1061 void 
1062G4PenelopeIonisationModel::CalculateDiscreteForPositrons(G4double kinEnergy,
1063                                                         G4double cutoffEnergy,
1064                                                         G4int Z,
1065                                                         G4double electronVolumeDensity)
1066{
1067  kineticEnergy1=kinEnergy;
1068  cosThetaPrimary=1.0;
1069  energySecondary=0.0;
1070  cosThetaSecondary=1.0;
1071  iOsc=-1;
1072  //constants
1073  G4double rb=kinEnergy+2.0*electron_mass_c2;
1074  G4double gamma = 1.0+kinEnergy/electron_mass_c2;
1075  G4double gamma2 = gamma*gamma;
1076  G4double beta2 = (gamma2-1.0)/gamma2;
1077  G4double amol = (gamma-1.0)*(gamma-1.0)/gamma2;
1078  G4double cps = kinEnergy*rb;
1079  G4double cp = std::sqrt(cps);
1080  G4double help = (gamma+1.0)*(gamma+1.0);
1081  G4double bha1 = amol*(2.0*help-1.0)/(gamma2-1.0);
1082  G4double bha2 = amol*(3.0+1.0/help);
1083  G4double bha3 = amol*2.0*gamma*(gamma-1.0)/help;
1084  G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/help;
1085 
1086  G4double delta = CalculateDeltaFermi(kinEnergy,Z,electronVolumeDensity);
1087  G4double distantTransvCS0 = std::max(std::log(gamma2)-beta2-delta,0.0);
1088
1089  G4double rl,rl1;
1090
1091  if (cutoffEnergy > kinEnergy) return; //delta rays are not generated
1092
1093  G4DataVector* qm = new G4DataVector();
1094  G4DataVector* cumulHardCS = new G4DataVector();
1095  std::vector<G4int> *typeOfInteraction = new std::vector<G4int>;
1096  G4DataVector* nbOfLevel = new G4DataVector();
1097 
1098
1099  //Hard close collisions with outer shells
1100  G4double wmaxc = kinEnergy;
1101  G4double closeCS0 = 0.0;
1102  G4double closeCS = 0.0;
1103  if (cutoffEnergy>0.1*eV) 
1104    {
1105      rl=cutoffEnergy/kinEnergy;
1106      rl1=1.0-rl;
1107      if (rl < 1.0)
1108        closeCS0 = (((1.0/rl)-1.0) + bha1*std::log(rl) + bha2*rl1
1109                    + (bha3/2.0)*((rl*rl)-1.0)
1110                    + (bha4/3.0)*(1.0-(rl*rl*rl)))/kinEnergy;
1111    }
1112
1113  // Cross sections for the different oscillators
1114
1115  // totalHardCS contains the cumulative hard interaction cross section for the different
1116  // excitable levels and the different interaction channels (close, distant, etc.),
1117  // i.e.
1118  // cumulHardCS[0] = 0.0
1119  // cumulHardCS[1] = 1st excitable level (distant longitudinal only)
1120  // cumulHardCS[2] = 1st excitable level (distant longitudinal + transverse)
1121  // cumulHardCS[3] = 1st excitable level (distant longitudinal + transverse + close)
1122  // cumulHardCS[4] = 1st excitable level (all channels) + 2nd excitable level (distant long only)
1123  // etc.
1124  // This is used for sampling the atomic level which is ionised and the channel of the
1125  // interaction.
1126  //
1127  // For each index iFill of the cumulHardCS vector,
1128  // nbOfLevel[iFill] contains the current excitable atomic level and
1129  // typeOfInteraction[iFill] contains the current interaction channel, with the legenda:
1130  //   1 = distant longitudinal interaction
1131  //   2 = distant transverse interaction
1132  //   3 = close collision
1133  //   4 = close collision with outer shells (in this case nbOfLevel < 0 --> no binding energy)
1134
1135
1136  G4int nOscil = ionizationEnergy->find(Z)->second->size();
1137  G4double totalHardCS = 0.0;
1138  G4double involvedElectrons = 0.0;
1139  for (G4int i=0;i<nOscil;i++){
1140    G4double wi = (*(resonanceEnergy->find(Z)->second))[i];
1141    G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
1142    //Distant excitations
1143    if (wi>cutoffEnergy && wi<kinEnergy)
1144      {
1145        if (wi>(1e-6*kinEnergy)){
1146          G4double cpp=std::sqrt((kinEnergy-wi)*(kinEnergy-wi+2.0*electron_mass_c2));
1147          qm->push_back(std::sqrt((cp-cpp)*(cp-cpp)+ electron_mass_c2 * electron_mass_c2)-electron_mass_c2);
1148        }
1149        else
1150          {
1151            qm->push_back(wi*wi/(beta2+2.0*electron_mass_c2));
1152          }
1153        //verificare che quando arriva qui il vettore ha SEMPRE l'i-esimo elemento
1154        if ((*qm)[i] < wi)
1155          {
1156           
1157            G4double distantLongitCS =  occupNb*std::log(wi*((*qm)[i]+2.0*electron_mass_c2)/
1158                                         ((*qm)[i]*(wi+2.0*electron_mass_c2)))/wi;
1159            cumulHardCS->push_back(totalHardCS);
1160            typeOfInteraction->push_back(1); //distant longitudinal
1161            nbOfLevel->push_back((G4double) i); //only excitable level are counted
1162            totalHardCS += distantLongitCS;
1163           
1164            G4double distantTransvCS = occupNb*distantTransvCS0/wi;
1165           
1166            cumulHardCS->push_back(totalHardCS);
1167            typeOfInteraction->push_back(2); //distant tranverse
1168            nbOfLevel->push_back((G4double) i);
1169            totalHardCS += distantTransvCS;
1170          }
1171      }
1172    else 
1173      qm->push_back(wi);
1174     
1175    //close collisions
1176    if(wi < wmaxc)
1177      {
1178        if (wi < cutoffEnergy) {
1179          involvedElectrons += occupNb;
1180        }
1181        else
1182          {
1183            rl=wi/kinEnergy;
1184            rl1=1.0-rl;
1185            closeCS = occupNb*(((1.0/rl)-1.0)+bha1*std::log(rl)+bha2*rl1
1186                               + (bha3/2.0)*((rl*rl)-1.0)
1187                               + (bha4/3.0)*(1.0-(rl*rl*rl)))/kinEnergy;
1188            cumulHardCS->push_back(totalHardCS);
1189            typeOfInteraction->push_back(3); //close
1190            nbOfLevel->push_back((G4double) i);
1191            totalHardCS += closeCS;
1192          }
1193      }
1194  } // loop on the levels
1195 
1196  cumulHardCS->push_back(totalHardCS);
1197  typeOfInteraction->push_back(4); //close interaction with outer shells
1198  nbOfLevel->push_back(-1.0);
1199  totalHardCS += involvedElectrons*closeCS0;
1200  cumulHardCS->push_back(totalHardCS); //this is the final value of the totalHardCS
1201
1202  if (totalHardCS < 1e-30) {
1203    kineticEnergy1=kinEnergy;
1204    cosThetaPrimary=1.0;
1205    energySecondary=0.0;
1206    cosThetaSecondary=0.0;
1207    iOsc=-1;
1208    delete qm;
1209    delete cumulHardCS;
1210    delete typeOfInteraction;
1211    delete nbOfLevel;
1212    return;
1213  }
1214
1215  //Selection of the active oscillator on the basis of the cumulative cross sections
1216  G4double TST = totalHardCS*G4UniformRand();
1217  G4int is=0;
1218  G4int js= nbOfLevel->size();
1219  do{
1220    G4int it=(is+js)/2;
1221    if (TST > (*cumulHardCS)[it]) is=it;
1222    if (TST <= (*cumulHardCS)[it]) js=it;
1223  }while((js-is) > 1);
1224
1225  G4double UII=0.0;
1226  G4double rkc=cutoffEnergy/kinEnergy;
1227  G4double dde;
1228  G4int kks;
1229
1230  G4int sampledInteraction = (*typeOfInteraction)[is];
1231  iOsc = (G4int) (*nbOfLevel)[is];
1232
1233  //Generates the final state according to the sampled level and
1234  //interaction channel
1235 
1236  if (sampledInteraction == 1)  //Hard distant longitudinal collisions
1237    {
1238      dde= (*(resonanceEnergy->find(Z)->second))[iOsc];
1239      kineticEnergy1=kinEnergy-dde;
1240      G4double qs=(*qm)[iOsc]/(1.0+((*qm)[iOsc]/(2.0*electron_mass_c2)));
1241      G4double q=qs/(std::pow((qs/dde)*(1.0+(0.5*dde/electron_mass_c2)),G4UniformRand())-(0.5*qs/electron_mass_c2));
1242      G4double qtrev = q*(q+2.0*electron_mass_c2);
1243      G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
1244      cosThetaPrimary = (cpps+cps-qtrev)/(2.0*cp*std::sqrt(cpps));
1245      if (cosThetaPrimary>1.0) cosThetaPrimary=1.0;
1246      //Energy and emission angle of the delta ray
1247      kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
1248      if (kks>4) 
1249        energySecondary=dde;
1250       
1251      else
1252        energySecondary=dde-(*(ionizationEnergy->find(Z)->second))[iOsc];
1253      cosThetaSecondary = 0.5*(dde*(kinEnergy+rb-dde)+qtrev)/std::sqrt(cps*qtrev);
1254      if (cosThetaSecondary>1.0) cosThetaSecondary=1.0;
1255    }
1256  else if (sampledInteraction == 2)  //Hard distant transverse collisions
1257    {
1258      dde=(*(resonanceEnergy->find(Z)->second))[iOsc];
1259      kineticEnergy1=kinEnergy-dde;
1260      cosThetaPrimary=1.0;
1261      //Energy and emission angle of the delta ray
1262      kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
1263      if (kks>4)
1264        {
1265          energySecondary=dde;
1266        }
1267      else
1268        {
1269          energySecondary=dde-(*(ionizationEnergy->find(Z)->second))[iOsc];
1270        }
1271      cosThetaSecondary = 1.0;
1272    }
1273  else if (sampledInteraction == 3 || sampledInteraction == 4) //Close interaction
1274    {
1275      if (sampledInteraction == 4) //interaction with inner shells
1276        {
1277          UII=0.0;
1278          rkc = cutoffEnergy/kinEnergy;
1279          iOsc = -1;
1280        }
1281      else
1282        {
1283          kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
1284          if (kks > 4) {
1285            UII=0.0;
1286          }
1287          else
1288            {
1289              UII = (*(ionizationEnergy->find(Z)->second))[iOsc];
1290            }
1291          rkc = (*(resonanceEnergy->find(Z)->second))[iOsc]/kinEnergy;
1292        }
1293      G4double phi,rk;
1294      do{
1295        rk=rkc/(1.0-G4UniformRand()*(1.0-rkc));
1296        phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
1297      }while ( G4UniformRand() > phi);
1298      //Energy and scattering angle (primary electron);
1299      kineticEnergy1 = kinEnergy*(1.0-rk);
1300      cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kinEnergy*(rb-(rk*kinEnergy))));
1301      //Energy and scattering angle of the delta ray
1302      energySecondary = kinEnergy-kineticEnergy1-UII;
1303      cosThetaSecondary = std::sqrt(rk*kinEnergy*rb/(kinEnergy*(rk*kinEnergy+2.0*electron_mass_c2)));
1304    }     
1305  else
1306    {
1307      G4String excep = "G4PenelopeIonisationModel - Error in the calculation of the final state";
1308      G4Exception(excep);
1309    }
1310
1311  delete qm;
1312  delete cumulHardCS;
1313  delete typeOfInteraction;
1314  delete nbOfLevel;
1315
1316  return;
1317}
1318
1319//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1320
1321G4double
1322G4PenelopeIonisationModel::CalculateCrossSectionsRatio(G4double kinEnergy,
1323                                                       G4double cutoffEnergy,
1324                                                       G4int Z, G4double electronVolumeDensity,
1325                                                       const G4ParticleDefinition* theParticle)
1326{
1327  //Constants
1328  G4double gamma = 1.0+kinEnergy/electron_mass_c2;
1329  G4double gamma2 = gamma*gamma;
1330  G4double beta2 = (gamma2-1.0)/gamma2;
1331  G4double constant = pi*classic_electr_radius*classic_electr_radius*2.0*electron_mass_c2/beta2;
1332  G4double delta = CalculateDeltaFermi(kinEnergy,Z,electronVolumeDensity);
1333  G4int nbOsc = (G4int) resonanceEnergy->find(Z)->second->size();
1334  G4double softCS = 0.0;
1335  G4double hardCS = 0.0;
1336  for (G4int i=0;i<nbOsc;i++){
1337    G4double resEnergy = (*(resonanceEnergy->find(Z)->second))[i];
1338    std::pair<G4double,G4double> SoftAndHardXS(0.,0.);
1339    if (theParticle == G4Electron::Electron())
1340      SoftAndHardXS = CrossSectionsRatioForElectrons(kinEnergy,resEnergy,delta,cutoffEnergy);
1341    else if (theParticle == G4Positron::Positron())
1342      SoftAndHardXS = CrossSectionsRatioForPositrons(kinEnergy,resEnergy,delta,cutoffEnergy);     
1343    G4double occupNb = (*(occupationNumber->find(Z)->second))[i];
1344    softCS += occupNb*constant*SoftAndHardXS.first;
1345    hardCS += occupNb*constant*SoftAndHardXS.second;
1346  }
1347  G4double ratio = 0.0;
1348 
1349  if (softCS+hardCS) ratio = (hardCS)/(softCS+hardCS);
1350  return ratio;
1351}
1352
1353//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1354
1355std::pair<G4double,G4double> 
1356G4PenelopeIonisationModel::CrossSectionsRatioForElectrons(G4double kineticEnergy,
1357                                                          G4double resEnergy,
1358                                                          G4double densityCorrection,
1359                                                          G4double cutoffEnergy)
1360{
1361  std::pair<G4double,G4double> theResult(0.,0.);
1362  if (kineticEnergy < resEnergy) return theResult;
1363
1364  //Calculate constants
1365  G4double gamma = 1.0+kineticEnergy/electron_mass_c2;
1366  G4double gamma2 = gamma*gamma;
1367  G4double beta2 = (gamma2-1.0)/gamma2;
1368  G4double cps = kineticEnergy*(kineticEnergy+2.0*electron_mass_c2);
1369  G4double amol = (kineticEnergy/(kineticEnergy+electron_mass_c2)) * (kineticEnergy/(kineticEnergy+electron_mass_c2)) ;
1370  G4double hardCont = 0.0;
1371  G4double softCont = 0.0;
1372   
1373  //Distant interactions
1374  G4double cp1s = (kineticEnergy-resEnergy)*(kineticEnergy-resEnergy+2.0*electron_mass_c2);
1375  G4double cp1 = std::sqrt(cp1s);
1376  G4double cp = std::sqrt(cps);
1377  G4double sdLong=0.0, sdTrans = 0.0, sdDist=0.0;
1378
1379  //Distant longitudinal interactions
1380  G4double qm = 0.0;
1381
1382  if (resEnergy > kineticEnergy*(1e-6))
1383    {
1384      qm = std::sqrt((cp-cp1)*(cp-cp1)+(electron_mass_c2*electron_mass_c2))-electron_mass_c2;
1385    }
1386  else
1387    {
1388      qm = resEnergy*resEnergy/(beta2*2.0*electron_mass_c2);
1389      qm = qm*(1.0-0.5*qm/electron_mass_c2);
1390    }
1391
1392  if (qm < resEnergy)
1393    {
1394      sdLong = std::log(resEnergy*(qm+2.0*electron_mass_c2)/(qm*(resEnergy+2.0*electron_mass_c2)));
1395    }
1396  else
1397    {
1398      sdLong = 0.0;
1399    }
1400 
1401  if (sdLong > 0) {
1402    sdTrans = std::max(std::log(gamma2)-beta2-densityCorrection,0.0);
1403    sdDist = sdTrans + sdLong;
1404    if (cutoffEnergy > resEnergy) 
1405      {
1406        softCont = sdDist/resEnergy;
1407      }
1408    else
1409      {
1410        hardCont = sdDist/resEnergy;
1411      }
1412  }
1413
1414  // Close collisions (Moeller's cross section)
1415  G4double wl = std::max(cutoffEnergy,resEnergy);
1416  G4double wu = 0.5*kineticEnergy;
1417 
1418  if (wl < (wu-1*eV)) 
1419    { 
1420      hardCont += (1.0/(kineticEnergy-wu))-(1.0/(kineticEnergy-wl))
1421        - (1.0/wu)+(1.0/wl)
1422        + (1.0-amol)*std::log(((kineticEnergy-wu)*wl)/((kineticEnergy-wl)*wu))/kineticEnergy
1423        + amol*(wu-wl)/(kineticEnergy*kineticEnergy);
1424      wu=wl;
1425    }
1426
1427  wl = resEnergy;
1428  if (wl > (wu-1*eV)) 
1429    {
1430      theResult.first = softCont;
1431      theResult.second = hardCont;
1432      return theResult;
1433    }
1434  softCont += (1.0/(kineticEnergy-wu))-(1.0/(kineticEnergy-wl))
1435        - (1.0/wu)+(1.0/wl)
1436        + (1.0-amol)*std::log(((kineticEnergy-wu)*wl)/((kineticEnergy-wl)*wu))/kineticEnergy
1437        + amol*(wu-wl)/(kineticEnergy*kineticEnergy);
1438  theResult.first = softCont;
1439  theResult.second = hardCont;
1440  return theResult;
1441}
1442
1443
1444//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1445
1446std::pair<G4double,G4double> 
1447G4PenelopeIonisationModel::CrossSectionsRatioForPositrons(G4double kineticEnergy,
1448                                                          G4double resEnergy,
1449                                                          G4double densityCorrection,
1450                                                          G4double cutoffEnergy)
1451{
1452
1453  std::pair<G4double,G4double> theResult(0.,0.);
1454  if (kineticEnergy < resEnergy) return theResult;
1455
1456  //Calculate constants
1457  G4double gamma = 1.0+kineticEnergy/electron_mass_c2;
1458  G4double gamma2 = gamma*gamma;
1459  G4double beta2 = (gamma2-1.0)/gamma2;
1460  G4double cps = kineticEnergy*(kineticEnergy+2.0*electron_mass_c2);
1461  G4double amol = (kineticEnergy/(kineticEnergy+electron_mass_c2)) * (kineticEnergy/(kineticEnergy+electron_mass_c2)) ;
1462  G4double help = (gamma+1.0)*(gamma+1.0);
1463  G4double bha1 = amol*(2.0*help-1.0)/(gamma2-1.0);
1464  G4double bha2 = amol*(3.0+1.0/help);
1465  G4double bha3 = amol*2.0*gamma*(gamma-1.0)/help;
1466  G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/help;
1467  G4double hardCont = 0.0;
1468  G4double softCont = 0.0;
1469
1470  //Distant interactions
1471  G4double cp1s = (kineticEnergy-resEnergy)*(kineticEnergy-resEnergy+2.0*electron_mass_c2);
1472  G4double cp1 = std::sqrt(cp1s);
1473  G4double cp = std::sqrt(cps);
1474  G4double sdLong=0.0, sdTrans = 0.0, sdDist=0.0;
1475
1476  //Distant longitudinal interactions
1477  G4double qm = 0.0;
1478
1479  if (resEnergy > kineticEnergy*(1e-6))
1480    {
1481      qm = std::sqrt((cp-cp1)*(cp-cp1)+(electron_mass_c2*electron_mass_c2))-electron_mass_c2;
1482    }
1483  else
1484    {
1485      qm = resEnergy*resEnergy/(beta2*2.0*electron_mass_c2);
1486      qm = qm*(1.0-0.5*qm/electron_mass_c2);
1487    }
1488
1489  if (qm < resEnergy)
1490    {
1491      sdLong = std::log(resEnergy*(qm+2.0*electron_mass_c2)/(qm*(resEnergy+2.0*electron_mass_c2)));
1492    }
1493  else
1494    {
1495      sdLong = 0.0;
1496    }
1497 
1498  if (sdLong > 0) {
1499    sdTrans = std::max(std::log(gamma2)-beta2-densityCorrection,0.0);
1500    sdDist = sdTrans + sdLong;
1501    if (cutoffEnergy > resEnergy) 
1502      {
1503        softCont = sdDist/resEnergy;
1504      }
1505    else
1506      {
1507        hardCont = sdDist/resEnergy;
1508      }
1509  }
1510
1511
1512  // Close collisions (Bhabha's cross section)
1513  G4double wl = std::max(cutoffEnergy,resEnergy);
1514  G4double wu = kineticEnergy;
1515 
1516  if (wl < (wu-1*eV)) {
1517    hardCont += (1.0/wl)-(1.0/wu)-bha1*std::log(wu/wl)/kineticEnergy
1518      + bha2*(wu-wl)/(kineticEnergy*kineticEnergy) 
1519      -bha3*((wu*wu)-(wl*wl))/(2.0*kineticEnergy*kineticEnergy*kineticEnergy)
1520      + bha4*((wu*wu*wu)-(wl*wl*wl))/(3.0*kineticEnergy*kineticEnergy*kineticEnergy*kineticEnergy);
1521    wu=wl;
1522  }
1523  wl = resEnergy;
1524  if (wl > (wu-1*eV)) 
1525    { 
1526      theResult.first = softCont;
1527      theResult.second = hardCont;
1528      return theResult;
1529    }
1530  softCont += (1.0/wl)-(1.0/wu)-bha1*std::log(wu/wl)/kineticEnergy
1531    + bha2*(wu-wl)/(kineticEnergy*kineticEnergy) 
1532    -bha3*((wu*wu)-(wl*wl))/(2.0*kineticEnergy*kineticEnergy*kineticEnergy)
1533    + bha4*((wu*wu*wu)-(wl*wl*wl))/(3.0*kineticEnergy*kineticEnergy*kineticEnergy*kineticEnergy);
1534 
1535
1536  theResult.first = softCont;
1537  theResult.second = hardCont;
1538  return theResult;
1539
1540}
1541
1542G4double G4PenelopeIonisationModel::ComputeStoppingPowerForElectrons(G4double kinEnergy,
1543                                                                     G4double cutEnergy,
1544                                                                     G4double deltaFermi,
1545                                                                     G4double resEnergy)
1546{
1547   //Calculate constants
1548  G4double gamma = 1.0+kinEnergy/electron_mass_c2;
1549  G4double gamma2 = gamma*gamma;
1550  G4double beta2 = (gamma2-1.0)/gamma2;
1551  G4double cps = kinEnergy*(kinEnergy+2.0*electron_mass_c2);
1552  G4double amol = (gamma-1.0)*(gamma-1.0)/gamma2;
1553  G4double sPower = 0.0;
1554  if (kinEnergy < resEnergy) return sPower;
1555
1556  //Distant interactions
1557  G4double cp1s = (kinEnergy-resEnergy)*(kinEnergy-resEnergy+2.0*electron_mass_c2);
1558  G4double cp1 = std::sqrt(cp1s);
1559  G4double cp = std::sqrt(cps);
1560  G4double sdLong=0.0, sdTrans = 0.0, sdDist=0.0;
1561
1562  //Distant longitudinal interactions
1563  G4double qm = 0.0;
1564
1565  if (resEnergy > kinEnergy*(1e-6))
1566    {
1567      qm = std::sqrt((cp-cp1)*(cp-cp1)+(electron_mass_c2*electron_mass_c2))-electron_mass_c2;
1568    }
1569  else
1570    {
1571      qm = resEnergy*resEnergy/(beta2*2.0*electron_mass_c2);
1572      qm = qm*(1.0-0.5*qm/electron_mass_c2);
1573    }
1574
1575  if (qm < resEnergy)
1576    sdLong = std::log(resEnergy*(qm+2.0*electron_mass_c2)/(qm*(resEnergy+2.0*electron_mass_c2)));
1577  else
1578    sdLong = 0.0;
1579 
1580  if (sdLong > 0) {
1581    sdTrans = std::max(std::log(gamma2)-beta2-deltaFermi,0.0);
1582    sdDist = sdTrans + sdLong;
1583    if (cutEnergy > resEnergy) sPower = sdDist;
1584  }
1585
1586
1587  // Close collisions (Moeller's cross section)
1588  G4double wl = std::max(cutEnergy,resEnergy);
1589  G4double wu = 0.5*kinEnergy;
1590 
1591  if (wl < (wu-1*eV)) wu=wl;
1592  wl = resEnergy;
1593  if (wl > (wu-1*eV)) return sPower;
1594  sPower += std::log(wu/wl)+(kinEnergy/(kinEnergy-wu))-(kinEnergy/(kinEnergy-wl))
1595    + (2.0 - amol)*std::log((kinEnergy-wu)/(kinEnergy-wl))
1596    + amol*((wu*wu)-(wl*wl))/(2.0*kinEnergy*kinEnergy);
1597  return sPower;
1598}
1599
1600
1601//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1602
1603G4double G4PenelopeIonisationModel::ComputeStoppingPowerForPositrons(G4double kinEnergy,
1604                                                                     G4double cutEnergy,
1605                                                                     G4double deltaFermi,
1606                                                                     G4double resEnergy)
1607{
1608 //Calculate constants
1609  G4double gamma = 1.0+kinEnergy/electron_mass_c2;
1610  G4double gamma2 = gamma*gamma;
1611  G4double beta2 = (gamma2-1.0)/gamma2;
1612  G4double cps = kinEnergy*(kinEnergy+2.0*electron_mass_c2);
1613  G4double amol = (kinEnergy/(kinEnergy+electron_mass_c2)) * (kinEnergy/(kinEnergy+electron_mass_c2));
1614  G4double help = (gamma+1.0)*(gamma+1.0);
1615  G4double bha1 = amol*(2.0*help-1.0)/(gamma2-1.0);
1616  G4double bha2 = amol*(3.0+1.0/help);
1617  G4double bha3 = amol*2.0*gamma*(gamma-1.0)/help;
1618  G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/help;
1619
1620  G4double sPower = 0.0;
1621  if (kinEnergy < resEnergy) return sPower;
1622
1623  //Distant interactions
1624  G4double cp1s = (kinEnergy-resEnergy)*(kinEnergy-resEnergy+2.0*electron_mass_c2);
1625  G4double cp1 = std::sqrt(cp1s);
1626  G4double cp = std::sqrt(cps);
1627  G4double sdLong=0.0, sdTrans = 0.0, sdDist=0.0;
1628
1629  //Distant longitudinal interactions
1630  G4double qm = 0.0;
1631
1632  if (resEnergy > kinEnergy*(1e-6))
1633    {
1634      qm = std::sqrt((cp-cp1)*(cp-cp1)+(electron_mass_c2*electron_mass_c2))-electron_mass_c2;
1635    }
1636  else
1637    {
1638      qm = resEnergy*resEnergy/(beta2*2.0*electron_mass_c2);
1639      qm = qm*(1.0-0.5*qm/electron_mass_c2);
1640    }
1641
1642  if (qm < resEnergy)
1643    sdLong = std::log(resEnergy*(qm+2.0*electron_mass_c2)/(qm*(resEnergy+2.0*electron_mass_c2)));
1644  else
1645    sdLong = 0.0;
1646   
1647  if (sdLong > 0) {
1648    sdTrans = std::max(std::log(gamma2)-beta2-deltaFermi,0.0);
1649    sdDist = sdTrans + sdLong;
1650    if (cutEnergy > resEnergy) sPower = sdDist;
1651  }
1652
1653
1654  // Close collisions (Bhabha's cross section)
1655  G4double wl = std::max(cutEnergy,resEnergy);
1656  G4double wu = kinEnergy;
1657 
1658  if (wl < (wu-1*eV)) wu=wl;
1659  wl = resEnergy;
1660  if (wl > (wu-1*eV)) return sPower;
1661  sPower += std::log(wu/wl)-bha1*(wu-wl)/kinEnergy
1662    + bha2*((wu*wu)-(wl*wl))/(2.0*kinEnergy*kinEnergy)
1663    - bha3*((wu*wu*wu)-(wl*wl*wl))/(3.0*kinEnergy*kinEnergy*kinEnergy)
1664    + bha4*((wu*wu*wu*wu)-(wl*wl*wl*wl))/(4.0*kinEnergy*kinEnergy*kinEnergy*kinEnergy);
1665  return sPower;
1666}
1667
1668//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1669
1670/* Notice: the methods here above are only temporary. They will become obsolete in a while */
1671
1672#include "G4VDataSetAlgorithm.hh"
1673#include "G4LinLogLogInterpolation.hh"
1674#include "G4CompositeEMDataSet.hh"
1675#include "G4EMDataSet.hh"
1676
1677std::vector<G4VEMDataSet*>* 
1678G4PenelopeIonisationModel::BuildCrossSectionTable(const G4ParticleDefinition* theParticle)
1679{
1680  std::vector<G4VEMDataSet*>* set = new std::vector<G4VEMDataSet*>;
1681
1682  size_t nOfBins = 200;
1683  //Temporary vector, a quick way to produce a log-spaced energy grid
1684  G4PhysicsLogVector* theLogVector = new G4PhysicsLogVector(LowEnergyLimit(),
1685                                                            HighEnergyLimit(), 
1686                                                            nOfBins);
1687  G4DataVector* energies;
1688  G4DataVector* cs;
1689 
1690  const G4ProductionCutsTable* theCoupleTable=
1691        G4ProductionCutsTable::GetProductionCutsTable();
1692  size_t numOfCouples = theCoupleTable->GetTableSize();
1693
1694
1695  for (size_t m=0; m<numOfCouples; m++) 
1696    {
1697      const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(m);
1698      const G4Material* material= couple->GetMaterial();
1699      const G4ElementVector* elementVector = material->GetElementVector();
1700      const G4double* nAtomsPerVolume = material->GetAtomicNumDensityVector();
1701      G4double electronVolumeDensity = 
1702        material->GetTotNbOfElectPerVolume();  //electron density
1703      G4int nElements = material->GetNumberOfElements();
1704 
1705      G4double tcut  = (*(theCoupleTable->GetEnergyCutsVector(1)))[couple->GetIndex()];
1706
1707      G4VDataSetAlgorithm* algo =  new G4LinLogLogInterpolation();
1708
1709      G4VEMDataSet* setForMat = new G4CompositeEMDataSet(algo,1.,1.);
1710
1711      for (G4int i=0; i<nElements; i++) 
1712        {
1713          G4int iZ = (G4int) (*elementVector)[i]->GetZ();
1714          energies = new G4DataVector;
1715          cs = new G4DataVector;
1716          G4double density = nAtomsPerVolume[i];         
1717          for (size_t bin=0; bin<nOfBins; bin++) 
1718            {
1719              G4double e = theLogVector->GetLowEdgeEnergy(bin);
1720              energies->push_back(e);
1721              G4double value = 0.0;           
1722              if(e > tcut) 
1723                {
1724                  G4double ratio = CalculateCrossSectionsRatio(e,tcut,iZ,
1725                                                               electronVolumeDensity,
1726                                                               theParticle);
1727                  value = crossSectionHandler->FindValue(iZ,e)*ratio*density;
1728                }
1729              cs->push_back(value);
1730            }
1731      G4VDataSetAlgorithm* algo1 = algo->Clone();
1732      G4VEMDataSet* elSet = new G4EMDataSet(i,energies,cs,algo1,1.,1.);
1733      setForMat->AddComponent(elSet);
1734    }
1735    set->push_back(setForMat);
1736  }
1737  delete theLogVector;
1738  return set;
1739}
1740
1741//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1742
1743G4int G4PenelopeIonisationModel::SampleRandomAtom(const G4MaterialCutsCouple* couple,
1744                                                     G4double e) const
1745{
1746  // Select randomly an element within the material, according to the weight
1747  // determined by the cross sections in the data set
1748
1749  const G4Material* material = couple->GetMaterial();
1750  G4int nElements = material->GetNumberOfElements();
1751
1752  // Special case: the material consists of one element
1753  if (nElements == 1)
1754    {
1755      G4int Z = (G4int) material->GetZ();
1756      return Z;
1757    }
1758
1759  // Composite material
1760  const G4ElementVector* elementVector = material->GetElementVector();
1761  size_t materialIndex = couple->GetIndex();
1762
1763  G4VEMDataSet* materialSet = (*theXSTable)[materialIndex];
1764  G4double materialCrossSection0 = 0.0;
1765  G4DataVector cross;
1766  cross.clear();
1767  for ( G4int i=0; i < nElements; i++ )
1768    {
1769      G4double cr = materialSet->GetComponent(i)->FindValue(e);
1770      materialCrossSection0 += cr;
1771      cross.push_back(materialCrossSection0);
1772    }
1773
1774  G4double random = G4UniformRand() * materialCrossSection0;
1775
1776  for (G4int k=0 ; k < nElements ; k++ )
1777    {
1778      if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
1779    }
1780  // It should never get here
1781  return 0;
1782}
1783
1784
1785//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1786
1787void G4PenelopeIonisationModel::ActivateAuger(G4bool augerbool)
1788{
1789  if (!DeexcitationFlag() && augerbool)
1790    {
1791      G4cout << "WARNING - G4PenelopeIonisationModel" << G4endl;
1792      G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
1793      G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
1794    }
1795  deexcitationManager.ActivateAugerElectronProduction(augerbool);
1796  if (verboseLevel > 1)
1797    G4cout << "Auger production set to " << augerbool << G4endl;
1798}
1799
1800//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.