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 9.6.p1)


  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$
 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                                                    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(0), 
 79   crossSectionHandler(nullptr), energySpectrum <<  79   isInitialised(false),
 80   isInitialised(false)                         <<  80   crossSectionHandler(0), energySpectrum(0)
 81 {                                                  81 {
 82   fIntrinsicLowEnergyLimit = 12.*eV;           <<  82   fIntrinsicLowEnergyLimit = 10.0*eV;
 83   fIntrinsicHighEnergyLimit = 100.0*GeV;           83   fIntrinsicHighEnergyLimit = 100.0*GeV;
 84                                                    84 
 85   verboseLevel = 0;                                85   verboseLevel = 0;
 86   SetAngularDistribution(new G4DeltaAngle());  << 
 87                                                    86 
 88   transitionManager = G4AtomicTransitionManage     87   transitionManager = G4AtomicTransitionManager::Instance();
 89 }                                                  88 }
 90                                                    89 
 91 //....oooOO0OOooo........oooOO0OOooo........oo     90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 92                                                    91 
 93 G4LivermoreIonisationModel::~G4LivermoreIonisa     92 G4LivermoreIonisationModel::~G4LivermoreIonisationModel()
 94 {                                                  93 {
 95   delete energySpectrum;                           94   delete energySpectrum;
 96   delete crossSectionHandler;                      95   delete crossSectionHandler;
 97 }                                                  96 }
 98                                                    97 
 99 //....oooOO0OOooo........oooOO0OOooo........oo     98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100                                                    99 
101 void G4LivermoreIonisationModel::Initialise(co    100 void G4LivermoreIonisationModel::Initialise(const G4ParticleDefinition* particle,
102               const G4DataVector& cuts)           101               const G4DataVector& cuts)
103 {                                                 102 {
104   //Check that the Livermore Ionisation is NOT    103   //Check that the Livermore Ionisation is NOT attached to e+
105   if (particle != G4Electron::Electron())         104   if (particle != G4Electron::Electron())
106     {                                             105     {
107       G4Exception("G4LivermoreIonisationModel:    106       G4Exception("G4LivermoreIonisationModel::Initialise",
108       "em0002",FatalException,                 << 107         "em0002",FatalException,"Livermore Ionisation Model is applicable only to electrons");
109       "Livermore Ionisation Model is applicabl << 
110     }                                             108     }
111   transitionManager->Initialise();             << 
112                                                   109 
113   //Read energy spectrum                          110   //Read energy spectrum
114   if (energySpectrum)                             111   if (energySpectrum) 
115     {                                             112     {
116       delete energySpectrum;                      113       delete energySpectrum;
117       energySpectrum = nullptr;                << 114       energySpectrum = 0;
118     }                                             115     }
119   energySpectrum = new G4eIonisationSpectrum()    116   energySpectrum = new G4eIonisationSpectrum();
120   if (verboseLevel > 3)                           117   if (verboseLevel > 3)
121     G4cout << "G4VEnergySpectrum is initialize    118     G4cout << "G4VEnergySpectrum is initialized" << G4endl;
122                                                   119 
123   //Initialize cross section handler              120   //Initialize cross section handler
124   if (crossSectionHandler)                        121   if (crossSectionHandler) 
125     {                                             122     {
126       delete crossSectionHandler;                 123       delete crossSectionHandler;
127       crossSectionHandler = nullptr;           << 124       crossSectionHandler = 0;
128     }                                             125     }
129                                                   126 
130   const size_t nbins = 20;                        127   const size_t nbins = 20;
131   G4double emin = LowEnergyLimit();               128   G4double emin = LowEnergyLimit();
132   G4double emax = HighEnergyLimit();              129   G4double emax = HighEnergyLimit();
133   G4int ndec = G4int(std::log10(emax/emin) + 0    130   G4int ndec = G4int(std::log10(emax/emin) + 0.5);
134   if(ndec <= 0) { ndec = 1; }                     131   if(ndec <= 0) { ndec = 1; }
135                                                   132 
136   G4VDataSetAlgorithm* interpolation = new G4S    133   G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation();
137   crossSectionHandler =                        << 134   crossSectionHandler = new G4eIonisationCrossSectionHandler(energySpectrum,interpolation,
138     new G4eIonisationCrossSectionHandler(energ << 135                    emin,emax,nbins*ndec);
139            emin,emax,nbins*ndec);              << 
140   crossSectionHandler->Clear();                   136   crossSectionHandler->Clear();
141   crossSectionHandler->LoadShellData("ioni/ion    137   crossSectionHandler->LoadShellData("ioni/ion-ss-cs-");
142   //This is used to retrieve cross section val    138   //This is used to retrieve cross section values later on
143   G4VEMDataSet* emdata =                          139   G4VEMDataSet* emdata = 
144     crossSectionHandler->BuildMeanFreePathForM    140     crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
145   //The method BuildMeanFreePathForMaterials()    141   //The method BuildMeanFreePathForMaterials() is required here only to force 
146   //the building of an internal table: the out    142   //the building of an internal table: the output pointer can be deleted
147   delete emdata;                                  143   delete emdata;  
148                                                   144 
149   if (verboseLevel > 0)                           145   if (verboseLevel > 0)
150     {                                             146     {
151       G4cout << "Livermore Ionisation model is    147       G4cout << "Livermore Ionisation model is initialized " << G4endl
152        << "Energy range: "                        148        << "Energy range: "
153        << LowEnergyLimit() / keV << " keV - "     149        << LowEnergyLimit() / keV << " keV - "
154        << HighEnergyLimit() / GeV << " GeV"       150        << HighEnergyLimit() / GeV << " GeV"
155        << G4endl;                                 151        << G4endl;
156     }                                             152     }
157                                                   153 
158   if (verboseLevel > 3)                           154   if (verboseLevel > 3)
159     {                                             155     {
160       G4cout << "Cross section data: " << G4en    156       G4cout << "Cross section data: " << G4endl; 
161       crossSectionHandler->PrintData();           157       crossSectionHandler->PrintData();
162       G4cout << "Parameters: " << G4endl;         158       G4cout << "Parameters: " << G4endl;
163       energySpectrum->PrintData();                159       energySpectrum->PrintData();
164     }                                             160     }
165                                                   161 
166   if(isInitialised) { return; }                   162   if(isInitialised) { return; }
167   fParticleChange = GetParticleChangeForLoss()    163   fParticleChange = GetParticleChangeForLoss();
168   isInitialised = true;                           164   isInitialised = true; 
169 }                                                 165 }
170                                                   166 
171 //....oooOO0OOooo........oooOO0OOooo........oo    167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172                                                   168 
173 G4double                                          169 G4double 
174 G4LivermoreIonisationModel::ComputeCrossSectio << 170 G4LivermoreIonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
175                             const G4ParticleDe << 171                    G4double energy,
176           G4double energy,                     << 172                    G4double Z, G4double,
177           G4double Z, G4double,                << 173                    G4double cutEnergy, 
178           G4double cutEnergy,                  << 174                    G4double)
179           G4double)                            << 
180 {                                                 175 {
181   G4int iZ = G4int(Z);                         << 176   G4int iZ = (G4int) Z;
182   if (!crossSectionHandler)                       177   if (!crossSectionHandler)
183     {                                             178     {
184       G4Exception("G4LivermoreIonisationModel:    179       G4Exception("G4LivermoreIonisationModel::ComputeCrossSectionPerAtom",
185       "em1007",FatalException,                 << 180         "em1007",FatalException,"The cross section handler is not correctly initialized");
186       "The cross section handler is not correc << 
187       return 0;                                   181       return 0;
188     }                                             182     }
189                                                   183   
190   //The cut is already included in the crossSe    184   //The cut is already included in the crossSectionHandler
191   G4double cs =                                   185   G4double cs = 
192     crossSectionHandler->GetCrossSectionAboveT << 186     crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,cutEnergy,iZ);
193                  cutEnergy,                    << 
194                  iZ);                          << 
195                                                   187 
196   if (verboseLevel > 1)                           188   if (verboseLevel > 1)
197     {                                             189     {
198       G4cout << "G4LivermoreIonisationModel "     190       G4cout << "G4LivermoreIonisationModel " << G4endl;
199       G4cout << "Cross section for delta emiss << 191       G4cout << "Cross section for delta emission > " << cutEnergy/keV << " keV at " <<
200        << cutEnergy/keV << " keV at "          << 192   energy/keV << " keV and Z = " << iZ << " --> " << cs/barn << " barn" << G4endl;
201        << energy/keV << " keV and Z = " << iZ  << 
202        << cs/barn << " barn" << G4endl;        << 
203     }                                             193     }
204   return cs;                                      194   return cs;
205 }                                                 195 }
206                                                   196 
207                                                   197 
208 //....oooOO0OOooo........oooOO0OOooo........oo    198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
209                                                   199 
210 G4double                                       << 200 G4double G4LivermoreIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
211 G4LivermoreIonisationModel::ComputeDEDXPerVolu << 201                 const G4ParticleDefinition* ,
212              const G4ParticleDefinition*,      << 202                 G4double kineticEnergy,
213              G4double kineticEnergy,           << 203                 G4double cutEnergy)
214              G4double cutEnergy)               << 
215 {                                                 204 {
216   G4double sPower = 0.0;                          205   G4double sPower = 0.0;
217                                                   206 
218   const G4ElementVector* theElementVector = ma    207   const G4ElementVector* theElementVector = material->GetElementVector();
219   size_t NumberOfElements = material->GetNumbe    208   size_t NumberOfElements = material->GetNumberOfElements() ;
220   const G4double* theAtomicNumDensityVector =     209   const G4double* theAtomicNumDensityVector =
221                     material->GetAtomicNumDens    210                     material->GetAtomicNumDensityVector();
222                                                   211 
223   // loop for elements in the material            212   // loop for elements in the material
224   for (size_t iel=0; iel<NumberOfElements; iel    213   for (size_t iel=0; iel<NumberOfElements; iel++ ) 
225     {                                             214     {
226       G4int iZ = (G4int)((*theElementVector)[i    215       G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
227       G4int nShells = transitionManager->Numbe    216       G4int nShells = transitionManager->NumberOfShells(iZ);
228       for (G4int n=0; n<nShells; n++)             217       for (G4int n=0; n<nShells; n++) 
229   {                                               218   {
230     G4double e = energySpectrum->AverageEnergy    219     G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
231                  kineticEnergy, n);               220                  kineticEnergy, n);
232     G4double cs= crossSectionHandler->FindValu    221     G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n);
233     sPower   += e * cs * theAtomicNumDensityVe    222     sPower   += e * cs * theAtomicNumDensityVector[iel];
234   }                                               223   }
235       G4double esp = energySpectrum->Excitatio    224       G4double esp = energySpectrum->Excitation(iZ,kineticEnergy);
236       sPower   += esp * theAtomicNumDensityVec    225       sPower   += esp * theAtomicNumDensityVector[iel];
237     }                                             226     }
238                                                   227 
239   if (verboseLevel > 2)                           228   if (verboseLevel > 2)
240     {                                             229     {
241       G4cout << "G4LivermoreIonisationModel "     230       G4cout << "G4LivermoreIonisationModel " << G4endl;
242       G4cout << "Stopping power < " << cutEner << 231       G4cout << "Stopping power < " << cutEnergy/keV << " keV at " << 
243        << " keV at " << kineticEnergy/keV << " << 232   kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl;
244        << sPower/(keV/mm) << " keV/mm" << G4en << 
245     }                                             233     }
246                                                   234   
247   return sPower;                                  235   return sPower;
248 }                                                 236 }
249                                                   237 
250 //....oooOO0OOooo........oooOO0OOooo........oo    238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
251                                                   239 
252 void G4LivermoreIonisationModel::SampleSeconda << 240 void G4LivermoreIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
253                                  std::vector<G << 241                const G4MaterialCutsCouple* couple,
254          const G4MaterialCutsCouple* couple,   << 242                const G4DynamicParticle* aDynamicParticle,
255          const G4DynamicParticle* aDynamicPart << 243                G4double cutE,
256          G4double cutE,                        << 244                G4double maxE)
257          G4double maxE)                        << 
258 {                                                 245 {
259                                                   246   
260   G4double kineticEnergy = aDynamicParticle->G    247   G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
261                                                   248 
262   if (kineticEnergy <= fIntrinsicLowEnergyLimi    249   if (kineticEnergy <= fIntrinsicLowEnergyLimit)
263     {                                             250     {
264       fParticleChange->SetProposedKineticEnerg    251       fParticleChange->SetProposedKineticEnergy(0.);
265       fParticleChange->ProposeLocalEnergyDepos    252       fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
266       return;                                     253       return;
267     }                                             254     }
268                                                   255 
269    // Select atom and shell                       256    // Select atom and shell
270   G4int Z = crossSectionHandler->SelectRandomA    257   G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
271   G4int shellIndex = crossSectionHandler->Sele    258   G4int shellIndex = crossSectionHandler->SelectRandomShell(Z, kineticEnergy);
272   const G4AtomicShell* shell = transitionManag    259   const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
273   G4double bindingEnergy = shell->BindingEnerg    260   G4double bindingEnergy = shell->BindingEnergy();
274                                                   261 
275   // Sample delta energy using energy interval    262   // Sample delta energy using energy interval for delta-electrons 
276   G4double energyMax =                            263   G4double energyMax = 
277     std::min(maxE,energySpectrum->MaxEnergyOfS    264     std::min(maxE,energySpectrum->MaxEnergyOfSecondaries(kineticEnergy));
278   G4double energyDelta = energySpectrum->Sampl    265   G4double energyDelta = energySpectrum->SampleEnergy(Z, cutE, energyMax,
279                   kineticEnergy, shellIndex);     266                   kineticEnergy, shellIndex);
280                                                   267 
281   if (energyDelta == 0.) //nothing happens        268   if (energyDelta == 0.) //nothing happens
282     { return; }                                   269     { return; }
283                                                   270 
284   const G4ParticleDefinition* electron = G4Ele << 271   // Transform to shell potential
285   G4DynamicParticle* delta = new G4DynamicPart << 272   G4double deltaKinE = energyDelta + 2.0*bindingEnergy;
286     GetAngularDistribution()->SampleDirectionF << 273   G4double primaryKinE = kineticEnergy + 2.0*bindingEnergy;
287                   Z, shellIndex,               << 274 
288                                                << 275   // sampling of scattering angle neglecting atomic motion
289                                                << 276   G4double deltaMom = std::sqrt(deltaKinE*(deltaKinE + 2.0*electron_mass_c2));
290                                                << 277   G4double primaryMom = std::sqrt(primaryKinE*(primaryKinE + 2.0*electron_mass_c2));
291   fvect->push_back(delta);                     << 278 
292                                                << 279   G4double cost = deltaKinE * (primaryKinE + 2.0*electron_mass_c2)
293   // Change kinematics of primary particle     << 280                             / (deltaMom * primaryMom);
294   G4ThreeVector direction = aDynamicParticle-> << 281   if (cost > 1.) { cost = 1.; }
295   G4double totalMomentum = std::sqrt(kineticEn << 282   G4double sint = std::sqrt((1. - cost)*(1. + cost));
296                                                << 283   G4double phi  = twopi * G4UniformRand();
297   G4ThreeVector finalP = totalMomentum*directi << 284   G4double dirx = sint * std::cos(phi);
298   finalP               = finalP.unit();        << 285   G4double diry = sint * std::sin(phi);
                                                   >> 286   G4double dirz = cost;
                                                   >> 287   
                                                   >> 288   // Rotate to incident electron direction
                                                   >> 289   G4ThreeVector primaryDirection = aDynamicParticle->GetMomentumDirection();
                                                   >> 290   G4ThreeVector deltaDir(dirx,diry,dirz);
                                                   >> 291   deltaDir.rotateUz(primaryDirection);
                                                   >> 292   //Updated components
                                                   >> 293   dirx = deltaDir.x();
                                                   >> 294   diry = deltaDir.y();
                                                   >> 295   dirz = deltaDir.z();
                                                   >> 296 
                                                   >> 297   // Take into account atomic motion del is relative momentum of the motion
                                                   >> 298   // kinetic energy of the motion == bindingEnergy in V.Ivanchenko model
                                                   >> 299   cost = 2.0*G4UniformRand() - 1.0;
                                                   >> 300   sint = std::sqrt(1. - cost*cost);
                                                   >> 301   phi  = twopi * G4UniformRand();
                                                   >> 302   G4double del = std::sqrt(bindingEnergy *(bindingEnergy + 2.0*electron_mass_c2))
                                                   >> 303                / deltaMom;
                                                   >> 304   dirx += del* sint * std::cos(phi);
                                                   >> 305   diry += del* sint * std::sin(phi);
                                                   >> 306   dirz += del* cost;
                                                   >> 307 
                                                   >> 308   // Find out new primary electron direction
                                                   >> 309   G4double finalPx = primaryMom*primaryDirection.x() - deltaMom*dirx;
                                                   >> 310   G4double finalPy = primaryMom*primaryDirection.y() - deltaMom*diry;
                                                   >> 311   G4double finalPz = primaryMom*primaryDirection.z() - deltaMom*dirz;
                                                   >> 312 
                                                   >> 313   //Ok, ready to create the delta ray
                                                   >> 314   G4DynamicParticle* theDeltaRay = new G4DynamicParticle();
                                                   >> 315   theDeltaRay->SetKineticEnergy(energyDelta);
                                                   >> 316   G4double norm = 1.0/std::sqrt(dirx*dirx + diry*diry + dirz*dirz);
                                                   >> 317   dirx *= norm;
                                                   >> 318   diry *= norm;
                                                   >> 319   dirz *= norm;
                                                   >> 320   theDeltaRay->SetMomentumDirection(dirx, diry, dirz);
                                                   >> 321   theDeltaRay->SetDefinition(G4Electron::Electron());
                                                   >> 322   fvect->push_back(theDeltaRay);
299                                                   323 
300   //This is the amount of energy available for    324   //This is the amount of energy available for fluorescence
301   G4double theEnergyDeposit = bindingEnergy;      325   G4double theEnergyDeposit = bindingEnergy;
302                                                   326 
303   // fill ParticleChange                          327   // fill ParticleChange
304   // changed energy and momentum of the actual    328   // changed energy and momentum of the actual particle
305   G4double finalKinEnergy = kineticEnergy - en    329   G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit;
306   if(finalKinEnergy < 0.0)                        330   if(finalKinEnergy < 0.0) 
307     {                                             331     {
308       theEnergyDeposit += finalKinEnergy;         332       theEnergyDeposit += finalKinEnergy;
309       finalKinEnergy    = 0.0;                    333       finalKinEnergy    = 0.0;
310     }                                             334     } 
311   else                                            335   else 
312     {                                             336     {
313       fParticleChange->ProposeMomentumDirectio << 337       G4double normLocal = 1.0/std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
                                                   >> 338       finalPx *= normLocal;
                                                   >> 339       finalPy *= normLocal;
                                                   >> 340       finalPz *= normLocal;
                                                   >> 341       fParticleChange->ProposeMomentumDirection(finalPx, finalPy, finalPz);
314     }                                             342     }
315   fParticleChange->SetProposedKineticEnergy(fi    343   fParticleChange->SetProposedKineticEnergy(finalKinEnergy);
316                                                   344 
317   if (theEnergyDeposit < 0)                       345   if (theEnergyDeposit < 0)
318     {                                             346     {
319       G4cout <<  "G4LivermoreIonisationModel:     347       G4cout <<  "G4LivermoreIonisationModel: Negative energy deposit: "
320        << theEnergyDeposit/eV << " eV" << G4en    348        << theEnergyDeposit/eV << " eV" << G4endl;
321       theEnergyDeposit = 0.0;                     349       theEnergyDeposit = 0.0;
322     }                                             350     }
323                                                   351 
324   //Assign local energy deposit                   352   //Assign local energy deposit
325   fParticleChange->ProposeLocalEnergyDeposit(t    353   fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit);
326                                                   354 
327   if (verboseLevel > 1)                           355   if (verboseLevel > 1)
328     {                                             356     {
329       G4cout << "-----------------------------    357       G4cout << "-----------------------------------------------------------" << G4endl;
330       G4cout << "Energy balance from G4Livermo    358       G4cout << "Energy balance from G4LivermoreIonisation" << G4endl;
331       G4cout << "Incoming primary energy: " <<    359       G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
332       G4cout << "-----------------------------    360       G4cout << "-----------------------------------------------------------" << G4endl;
333       G4cout << "Outgoing primary energy: " <<    361       G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl;
334       G4cout << "Delta ray " << energyDelta/ke    362       G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl;
335       G4cout << "Fluorescence: " << (bindingEn    363       G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl;
336       G4cout << "Local energy deposit " << the    364       G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl;
337       G4cout << "Total final state: " << (fina    365       G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy)
338             << " keV" << G4endl;                  366             << " keV" << G4endl;
339       G4cout << "-----------------------------    367       G4cout << "-----------------------------------------------------------" << G4endl;
340     }                                             368     }
341   return;                                         369   return;
342 }                                                 370 }
343                                                   371 
344 //....oooOO0OOooo........oooOO0OOooo........oo    372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345                                                   373 
346                                                   374