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 11.1)


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