Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4LivermoreIonisationModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/lowenergy/src/G4LivermoreIonisationModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4LivermoreIonisationModel.cc (Version 10.1.p3)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: G4LivermoreIonisationModel.cc 87443 2014-12-04 12:26:31Z gunter $
 26 //                                                 27 //
 27 // Author: Luciano Pandola                         28 // Author: Luciano Pandola
 28 //         on base of G4LowEnergyIonisation de     29 //         on base of G4LowEnergyIonisation developed by A.Forti and V.Ivanchenko
 29 //                                                 30 //
 30 // History:                                        31 // History:
 31 // --------                                        32 // --------
 32 // 12 Jan 2009   L Pandola    Migration from p     33 // 12 Jan 2009   L Pandola    Migration from process to model 
 33 // 03 Mar 2009   L Pandola    Bug fix (release     34 // 03 Mar 2009   L Pandola    Bug fix (release memory in the destructor)
 34 // 15 Apr 2009   V Ivanchenko Cleanup initiali     35 // 15 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
 35 //                  - apply internal high-ener     36 //                  - apply internal high-energy limit only in constructor 
 36 //                  - do not apply low-energy      37 //                  - do not apply low-energy limit (default is 0)
 37 //                  - simplify sampling of dee     38 //                  - simplify sampling of deexcitation by using cut in energy
 38 //                  - set activation of Auger      39 //                  - set activation of Auger "false"
 39 //                  - remove initialisation of     40 //                  - remove initialisation of element selectors
 40 // 19 May 2009   L Pandola    Explicitely set      41 // 19 May 2009   L Pandola    Explicitely set to zero pointers deleted in 
 41 //                            Initialise(), si     42 //                            Initialise(), since they might be checked later on
 42 // 23 Oct 2009   L Pandola                         43 // 23 Oct 2009   L Pandola
 43 //                  - atomic deexcitation mana     44 //                  - atomic deexcitation managed via G4VEmModel::DeexcitationFlag() is
 44 //                    set as "true" (default w     45 //                    set as "true" (default would be false)
 45 // 12 Oct 2010   L Pandola                         46 // 12 Oct 2010   L Pandola
 46 //                  - add debugging informatio     47 //                  - add debugging information about energy in 
 47 //                    SampleDeexcitationAlongS     48 //                    SampleDeexcitationAlongStep()
 48 //                  - generate fluorescence Sa     49 //                  - generate fluorescence SampleDeexcitationAlongStep() only above 
 49 //                    the cuts.                    50 //                    the cuts.
 50 // 01 Jun 2011   V Ivanchenko general cleanup      51 // 01 Jun 2011   V Ivanchenko general cleanup - all old deexcitation code removed
 51 //                                                 52 //
 52                                                    53 
 53 #include "G4LivermoreIonisationModel.hh"           54 #include "G4LivermoreIonisationModel.hh"
 54 #include "G4PhysicalConstants.hh"                  55 #include "G4PhysicalConstants.hh"
 55 #include "G4SystemOfUnits.hh"                      56 #include "G4SystemOfUnits.hh"
 56 #include "G4ParticleDefinition.hh"                 57 #include "G4ParticleDefinition.hh"
 57 #include "G4MaterialCutsCouple.hh"                 58 #include "G4MaterialCutsCouple.hh"
 58 #include "G4ProductionCutsTable.hh"                59 #include "G4ProductionCutsTable.hh"
 59 #include "G4DynamicParticle.hh"                    60 #include "G4DynamicParticle.hh"
 60 #include "G4Element.hh"                            61 #include "G4Element.hh"
 61 #include "G4ParticleChangeForLoss.hh"              62 #include "G4ParticleChangeForLoss.hh"
 62 #include "G4Electron.hh"                           63 #include "G4Electron.hh"
 63 #include "G4CrossSectionHandler.hh"                64 #include "G4CrossSectionHandler.hh"
 64 #include "G4VEMDataSet.hh"                         65 #include "G4VEMDataSet.hh"
 65 #include "G4eIonisationCrossSectionHandler.hh"     66 #include "G4eIonisationCrossSectionHandler.hh"
 66 #include "G4eIonisationSpectrum.hh"                67 #include "G4eIonisationSpectrum.hh"
 67 #include "G4VEnergySpectrum.hh"                    68 #include "G4VEnergySpectrum.hh"
 68 #include "G4SemiLogInterpolation.hh"               69 #include "G4SemiLogInterpolation.hh"
 69 #include "G4AtomicTransitionManager.hh"            70 #include "G4AtomicTransitionManager.hh"
 70 #include "G4AtomicShell.hh"                        71 #include "G4AtomicShell.hh"
 71 #include "G4DeltaAngle.hh"                         72 #include "G4DeltaAngle.hh"
 72                                                    73 
 73 //....oooOO0OOooo........oooOO0OOooo........oo     74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 74                                                    75 
 75                                                    76 
 76 G4LivermoreIonisationModel::G4LivermoreIonisat     77 G4LivermoreIonisationModel::G4LivermoreIonisationModel(const G4ParticleDefinition*,
 77                    const G4String& nam) :          78                    const G4String& nam) : 
 78   G4VEmModel(nam), fParticleChange(nullptr),   <<  79   G4VEmModel(nam), fParticleChange(0), 
 79   crossSectionHandler(nullptr), energySpectrum <<  80   isInitialised(false),
 80   isInitialised(false)                         <<  81   crossSectionHandler(0), energySpectrum(0)
 81 {                                                  82 {
 82   fIntrinsicLowEnergyLimit = 12.*eV;           <<  83   fIntrinsicLowEnergyLimit = 10.0*eV;
 83   fIntrinsicHighEnergyLimit = 100.0*GeV;           84   fIntrinsicHighEnergyLimit = 100.0*GeV;
 84                                                    85 
 85   verboseLevel = 0;                                86   verboseLevel = 0;
 86   SetAngularDistribution(new G4DeltaAngle());      87   SetAngularDistribution(new G4DeltaAngle());
 87                                                    88 
 88   transitionManager = G4AtomicTransitionManage     89   transitionManager = G4AtomicTransitionManager::Instance();
 89 }                                                  90 }
 90                                                    91 
 91 //....oooOO0OOooo........oooOO0OOooo........oo     92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 92                                                    93 
 93 G4LivermoreIonisationModel::~G4LivermoreIonisa     94 G4LivermoreIonisationModel::~G4LivermoreIonisationModel()
 94 {                                                  95 {
 95   delete energySpectrum;                           96   delete energySpectrum;
 96   delete crossSectionHandler;                      97   delete crossSectionHandler;
 97 }                                                  98 }
 98                                                    99 
 99 //....oooOO0OOooo........oooOO0OOooo........oo    100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100                                                   101 
101 void G4LivermoreIonisationModel::Initialise(co    102 void G4LivermoreIonisationModel::Initialise(const G4ParticleDefinition* particle,
102               const G4DataVector& cuts)           103               const G4DataVector& cuts)
103 {                                                 104 {
104   //Check that the Livermore Ionisation is NOT    105   //Check that the Livermore Ionisation is NOT attached to e+
105   if (particle != G4Electron::Electron())         106   if (particle != G4Electron::Electron())
106     {                                             107     {
107       G4Exception("G4LivermoreIonisationModel:    108       G4Exception("G4LivermoreIonisationModel::Initialise",
108       "em0002",FatalException,                    109       "em0002",FatalException,
109       "Livermore Ionisation Model is applicabl    110       "Livermore Ionisation Model is applicable only to electrons");
110     }                                             111     }
                                                   >> 112 
111   transitionManager->Initialise();                113   transitionManager->Initialise();
112                                                   114 
113   //Read energy spectrum                          115   //Read energy spectrum
114   if (energySpectrum)                             116   if (energySpectrum) 
115     {                                             117     {
116       delete energySpectrum;                      118       delete energySpectrum;
117       energySpectrum = nullptr;                << 119       energySpectrum = 0;
118     }                                             120     }
119   energySpectrum = new G4eIonisationSpectrum()    121   energySpectrum = new G4eIonisationSpectrum();
120   if (verboseLevel > 3)                           122   if (verboseLevel > 3)
121     G4cout << "G4VEnergySpectrum is initialize    123     G4cout << "G4VEnergySpectrum is initialized" << G4endl;
122                                                   124 
123   //Initialize cross section handler              125   //Initialize cross section handler
124   if (crossSectionHandler)                        126   if (crossSectionHandler) 
125     {                                             127     {
126       delete crossSectionHandler;                 128       delete crossSectionHandler;
127       crossSectionHandler = nullptr;           << 129       crossSectionHandler = 0;
128     }                                             130     }
129                                                   131 
130   const size_t nbins = 20;                        132   const size_t nbins = 20;
131   G4double emin = LowEnergyLimit();               133   G4double emin = LowEnergyLimit();
132   G4double emax = HighEnergyLimit();              134   G4double emax = HighEnergyLimit();
133   G4int ndec = G4int(std::log10(emax/emin) + 0    135   G4int ndec = G4int(std::log10(emax/emin) + 0.5);
134   if(ndec <= 0) { ndec = 1; }                     136   if(ndec <= 0) { ndec = 1; }
135                                                   137 
136   G4VDataSetAlgorithm* interpolation = new G4S    138   G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation();
137   crossSectionHandler =                           139   crossSectionHandler = 
138     new G4eIonisationCrossSectionHandler(energ    140     new G4eIonisationCrossSectionHandler(energySpectrum,interpolation,
139            emin,emax,nbins*ndec);                 141            emin,emax,nbins*ndec);
140   crossSectionHandler->Clear();                   142   crossSectionHandler->Clear();
141   crossSectionHandler->LoadShellData("ioni/ion    143   crossSectionHandler->LoadShellData("ioni/ion-ss-cs-");
142   //This is used to retrieve cross section val    144   //This is used to retrieve cross section values later on
143   G4VEMDataSet* emdata =                          145   G4VEMDataSet* emdata = 
144     crossSectionHandler->BuildMeanFreePathForM    146     crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
145   //The method BuildMeanFreePathForMaterials()    147   //The method BuildMeanFreePathForMaterials() is required here only to force 
146   //the building of an internal table: the out    148   //the building of an internal table: the output pointer can be deleted
147   delete emdata;                                  149   delete emdata;  
148                                                   150 
149   if (verboseLevel > 0)                           151   if (verboseLevel > 0)
150     {                                             152     {
151       G4cout << "Livermore Ionisation model is    153       G4cout << "Livermore Ionisation model is initialized " << G4endl
152        << "Energy range: "                        154        << "Energy range: "
153        << LowEnergyLimit() / keV << " keV - "     155        << LowEnergyLimit() / keV << " keV - "
154        << HighEnergyLimit() / GeV << " GeV"       156        << HighEnergyLimit() / GeV << " GeV"
155        << G4endl;                                 157        << G4endl;
156     }                                             158     }
157                                                   159 
158   if (verboseLevel > 3)                           160   if (verboseLevel > 3)
159     {                                             161     {
160       G4cout << "Cross section data: " << G4en    162       G4cout << "Cross section data: " << G4endl; 
161       crossSectionHandler->PrintData();           163       crossSectionHandler->PrintData();
162       G4cout << "Parameters: " << G4endl;         164       G4cout << "Parameters: " << G4endl;
163       energySpectrum->PrintData();                165       energySpectrum->PrintData();
164     }                                             166     }
165                                                   167 
166   if(isInitialised) { return; }                   168   if(isInitialised) { return; }
167   fParticleChange = GetParticleChangeForLoss()    169   fParticleChange = GetParticleChangeForLoss();
168   isInitialised = true;                           170   isInitialised = true; 
169 }                                                 171 }
170                                                   172 
171 //....oooOO0OOooo........oooOO0OOooo........oo    173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172                                                   174 
173 G4double                                          175 G4double 
174 G4LivermoreIonisationModel::ComputeCrossSectio    176 G4LivermoreIonisationModel::ComputeCrossSectionPerAtom(
175                             const G4ParticleDe    177                             const G4ParticleDefinition*,
176           G4double energy,                        178           G4double energy,
177           G4double Z, G4double,                   179           G4double Z, G4double,
178           G4double cutEnergy,                     180           G4double cutEnergy, 
179           G4double)                               181           G4double)
180 {                                                 182 {
181   G4int iZ = G4int(Z);                         << 183   G4int iZ = (G4int) Z;
182   if (!crossSectionHandler)                       184   if (!crossSectionHandler)
183     {                                             185     {
184       G4Exception("G4LivermoreIonisationModel:    186       G4Exception("G4LivermoreIonisationModel::ComputeCrossSectionPerAtom",
185       "em1007",FatalException,                    187       "em1007",FatalException,
186       "The cross section handler is not correc    188       "The cross section handler is not correctly initialized");
187       return 0;                                   189       return 0;
188     }                                             190     }
189                                                   191   
190   //The cut is already included in the crossSe    192   //The cut is already included in the crossSectionHandler
191   G4double cs =                                   193   G4double cs = 
192     crossSectionHandler->GetCrossSectionAboveT    194     crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,
193                  cutEnergy,                       195                  cutEnergy,
194                  iZ);                             196                  iZ);
195                                                   197 
196   if (verboseLevel > 1)                           198   if (verboseLevel > 1)
197     {                                             199     {
198       G4cout << "G4LivermoreIonisationModel "     200       G4cout << "G4LivermoreIonisationModel " << G4endl;
199       G4cout << "Cross section for delta emiss    201       G4cout << "Cross section for delta emission > " 
200        << cutEnergy/keV << " keV at "             202        << cutEnergy/keV << " keV at " 
201        << energy/keV << " keV and Z = " << iZ     203        << energy/keV << " keV and Z = " << iZ << " --> " 
202        << cs/barn << " barn" << G4endl;           204        << cs/barn << " barn" << G4endl;
203     }                                             205     }
204   return cs;                                      206   return cs;
205 }                                                 207 }
206                                                   208 
207                                                   209 
208 //....oooOO0OOooo........oooOO0OOooo........oo    210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
209                                                   211 
210 G4double                                          212 G4double 
211 G4LivermoreIonisationModel::ComputeDEDXPerVolu    213 G4LivermoreIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
212              const G4ParticleDefinition*,         214              const G4ParticleDefinition*,
213              G4double kineticEnergy,              215              G4double kineticEnergy,
214              G4double cutEnergy)                  216              G4double cutEnergy)
215 {                                                 217 {
216   G4double sPower = 0.0;                          218   G4double sPower = 0.0;
217                                                   219 
218   const G4ElementVector* theElementVector = ma    220   const G4ElementVector* theElementVector = material->GetElementVector();
219   size_t NumberOfElements = material->GetNumbe    221   size_t NumberOfElements = material->GetNumberOfElements() ;
220   const G4double* theAtomicNumDensityVector =     222   const G4double* theAtomicNumDensityVector =
221                     material->GetAtomicNumDens    223                     material->GetAtomicNumDensityVector();
222                                                   224 
223   // loop for elements in the material            225   // loop for elements in the material
224   for (size_t iel=0; iel<NumberOfElements; iel    226   for (size_t iel=0; iel<NumberOfElements; iel++ ) 
225     {                                             227     {
226       G4int iZ = (G4int)((*theElementVector)[i    228       G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
227       G4int nShells = transitionManager->Numbe    229       G4int nShells = transitionManager->NumberOfShells(iZ);
228       for (G4int n=0; n<nShells; n++)             230       for (G4int n=0; n<nShells; n++) 
229   {                                               231   {
230     G4double e = energySpectrum->AverageEnergy    232     G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
231                  kineticEnergy, n);               233                  kineticEnergy, n);
232     G4double cs= crossSectionHandler->FindValu    234     G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n);
233     sPower   += e * cs * theAtomicNumDensityVe    235     sPower   += e * cs * theAtomicNumDensityVector[iel];
234   }                                               236   }
235       G4double esp = energySpectrum->Excitatio    237       G4double esp = energySpectrum->Excitation(iZ,kineticEnergy);
236       sPower   += esp * theAtomicNumDensityVec    238       sPower   += esp * theAtomicNumDensityVector[iel];
237     }                                             239     }
238                                                   240 
239   if (verboseLevel > 2)                           241   if (verboseLevel > 2)
240     {                                             242     {
241       G4cout << "G4LivermoreIonisationModel "     243       G4cout << "G4LivermoreIonisationModel " << G4endl;
242       G4cout << "Stopping power < " << cutEner    244       G4cout << "Stopping power < " << cutEnergy/keV 
243        << " keV at " << kineticEnergy/keV << "    245        << " keV at " << kineticEnergy/keV << " keV = " 
244        << sPower/(keV/mm) << " keV/mm" << G4en    246        << sPower/(keV/mm) << " keV/mm" << G4endl;
245     }                                             247     }
246                                                   248   
247   return sPower;                                  249   return sPower;
248 }                                                 250 }
249                                                   251 
250 //....oooOO0OOooo........oooOO0OOooo........oo    252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
251                                                   253 
252 void G4LivermoreIonisationModel::SampleSeconda    254 void G4LivermoreIonisationModel::SampleSecondaries(
253                                  std::vector<G    255                                  std::vector<G4DynamicParticle*>* fvect,
254          const G4MaterialCutsCouple* couple,      256          const G4MaterialCutsCouple* couple,
255          const G4DynamicParticle* aDynamicPart    257          const G4DynamicParticle* aDynamicParticle,
256          G4double cutE,                           258          G4double cutE,
257          G4double maxE)                           259          G4double maxE)
258 {                                                 260 {
259                                                   261   
260   G4double kineticEnergy = aDynamicParticle->G    262   G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
261                                                   263 
262   if (kineticEnergy <= fIntrinsicLowEnergyLimi    264   if (kineticEnergy <= fIntrinsicLowEnergyLimit)
263     {                                             265     {
264       fParticleChange->SetProposedKineticEnerg    266       fParticleChange->SetProposedKineticEnergy(0.);
265       fParticleChange->ProposeLocalEnergyDepos    267       fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
266       return;                                     268       return;
267     }                                             269     }
268                                                   270 
269    // Select atom and shell                       271    // Select atom and shell
270   G4int Z = crossSectionHandler->SelectRandomA    272   G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
271   G4int shellIndex = crossSectionHandler->Sele    273   G4int shellIndex = crossSectionHandler->SelectRandomShell(Z, kineticEnergy);
272   const G4AtomicShell* shell = transitionManag    274   const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
273   G4double bindingEnergy = shell->BindingEnerg    275   G4double bindingEnergy = shell->BindingEnergy();
274                                                   276 
275   // Sample delta energy using energy interval    277   // Sample delta energy using energy interval for delta-electrons 
276   G4double energyMax =                            278   G4double energyMax = 
277     std::min(maxE,energySpectrum->MaxEnergyOfS    279     std::min(maxE,energySpectrum->MaxEnergyOfSecondaries(kineticEnergy));
278   G4double energyDelta = energySpectrum->Sampl    280   G4double energyDelta = energySpectrum->SampleEnergy(Z, cutE, energyMax,
279                   kineticEnergy, shellIndex);     281                   kineticEnergy, shellIndex);
280                                                   282 
281   if (energyDelta == 0.) //nothing happens        283   if (energyDelta == 0.) //nothing happens
282     { return; }                                   284     { return; }
283                                                   285 
284   const G4ParticleDefinition* electron = G4Ele    286   const G4ParticleDefinition* electron = G4Electron::Electron();
285   G4DynamicParticle* delta = new G4DynamicPart    287   G4DynamicParticle* delta = new G4DynamicParticle(electron, 
286     GetAngularDistribution()->SampleDirectionF    288     GetAngularDistribution()->SampleDirectionForShell(aDynamicParticle, energyDelta,
287                   Z, shellIndex,                  289                   Z, shellIndex,
288                                                   290                                                       couple->GetMaterial()),
289                                                   291                                                       energyDelta);
290                                                   292 
291   fvect->push_back(delta);                        293   fvect->push_back(delta);
292                                                   294 
293   // Change kinematics of primary particle        295   // Change kinematics of primary particle
294   G4ThreeVector direction = aDynamicParticle->    296   G4ThreeVector direction = aDynamicParticle->GetMomentumDirection();
295   G4double totalMomentum = std::sqrt(kineticEn    297   G4double totalMomentum = std::sqrt(kineticEnergy*(kineticEnergy + 2*electron_mass_c2));
296                                                   298 
297   G4ThreeVector finalP = totalMomentum*directi    299   G4ThreeVector finalP = totalMomentum*direction - delta->GetMomentum();
298   finalP               = finalP.unit();           300   finalP               = finalP.unit();
299                                                   301 
300   //This is the amount of energy available for    302   //This is the amount of energy available for fluorescence
301   G4double theEnergyDeposit = bindingEnergy;      303   G4double theEnergyDeposit = bindingEnergy;
302                                                   304 
303   // fill ParticleChange                          305   // fill ParticleChange
304   // changed energy and momentum of the actual    306   // changed energy and momentum of the actual particle
305   G4double finalKinEnergy = kineticEnergy - en    307   G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit;
306   if(finalKinEnergy < 0.0)                        308   if(finalKinEnergy < 0.0) 
307     {                                             309     {
308       theEnergyDeposit += finalKinEnergy;         310       theEnergyDeposit += finalKinEnergy;
309       finalKinEnergy    = 0.0;                    311       finalKinEnergy    = 0.0;
310     }                                             312     } 
311   else                                            313   else 
312     {                                             314     {
313       fParticleChange->ProposeMomentumDirectio    315       fParticleChange->ProposeMomentumDirection(finalP);
314     }                                             316     }
315   fParticleChange->SetProposedKineticEnergy(fi    317   fParticleChange->SetProposedKineticEnergy(finalKinEnergy);
316                                                   318 
317   if (theEnergyDeposit < 0)                       319   if (theEnergyDeposit < 0)
318     {                                             320     {
319       G4cout <<  "G4LivermoreIonisationModel:     321       G4cout <<  "G4LivermoreIonisationModel: Negative energy deposit: "
320        << theEnergyDeposit/eV << " eV" << G4en    322        << theEnergyDeposit/eV << " eV" << G4endl;
321       theEnergyDeposit = 0.0;                     323       theEnergyDeposit = 0.0;
322     }                                             324     }
323                                                   325 
324   //Assign local energy deposit                   326   //Assign local energy deposit
325   fParticleChange->ProposeLocalEnergyDeposit(t    327   fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit);
326                                                   328 
327   if (verboseLevel > 1)                           329   if (verboseLevel > 1)
328     {                                             330     {
329       G4cout << "-----------------------------    331       G4cout << "-----------------------------------------------------------" << G4endl;
330       G4cout << "Energy balance from G4Livermo    332       G4cout << "Energy balance from G4LivermoreIonisation" << G4endl;
331       G4cout << "Incoming primary energy: " <<    333       G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
332       G4cout << "-----------------------------    334       G4cout << "-----------------------------------------------------------" << G4endl;
333       G4cout << "Outgoing primary energy: " <<    335       G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl;
334       G4cout << "Delta ray " << energyDelta/ke    336       G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl;
335       G4cout << "Fluorescence: " << (bindingEn    337       G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl;
336       G4cout << "Local energy deposit " << the    338       G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl;
337       G4cout << "Total final state: " << (fina    339       G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy)
338             << " keV" << G4endl;                  340             << " keV" << G4endl;
339       G4cout << "-----------------------------    341       G4cout << "-----------------------------------------------------------" << G4endl;
340     }                                             342     }
341   return;                                         343   return;
342 }                                                 344 }
343                                                   345 
344 //....oooOO0OOooo........oooOO0OOooo........oo    346 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345                                                   347 
346                                                   348