| Geant4 Cross Reference |
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 74822 2013-10-22 14:42:13Z gcosmo $
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 }
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 = 0;
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 = 0;
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 // Transform to shell potential
285 G4DynamicParticle* delta = new G4DynamicPart << 285 G4double deltaKinE = energyDelta + 2.0*bindingEnergy;
286 GetAngularDistribution()->SampleDirectionF << 286 G4double primaryKinE = kineticEnergy + 2.0*bindingEnergy;
287 Z, shellIndex, << 287
288 << 288 // sampling of scattering angle neglecting atomic motion
289 << 289 G4double deltaMom = std::sqrt(deltaKinE*(deltaKinE + 2.0*electron_mass_c2));
290 << 290 G4double primaryMom = std::sqrt(primaryKinE*(primaryKinE + 2.0*electron_mass_c2));
291 fvect->push_back(delta); << 291
292 << 292 G4double cost = deltaKinE * (primaryKinE + 2.0*electron_mass_c2)
293 // Change kinematics of primary particle << 293 / (deltaMom * primaryMom);
294 G4ThreeVector direction = aDynamicParticle-> << 294 if (cost > 1.) { cost = 1.; }
295 G4double totalMomentum = std::sqrt(kineticEn << 295 G4double sint = std::sqrt((1. - cost)*(1. + cost));
296 << 296 G4double phi = twopi * G4UniformRand();
297 G4ThreeVector finalP = totalMomentum*directi << 297 G4double dirx = sint * std::cos(phi);
298 finalP = finalP.unit(); << 298 G4double diry = sint * std::sin(phi);
>> 299 G4double dirz = cost;
>> 300
>> 301 // Rotate to incident electron direction
>> 302 G4ThreeVector primaryDirection = aDynamicParticle->GetMomentumDirection();
>> 303 G4ThreeVector deltaDir(dirx,diry,dirz);
>> 304 deltaDir.rotateUz(primaryDirection);
>> 305 //Updated components
>> 306 dirx = deltaDir.x();
>> 307 diry = deltaDir.y();
>> 308 dirz = deltaDir.z();
>> 309
>> 310 // Take into account atomic motion del is relative momentum of the motion
>> 311 // kinetic energy of the motion == bindingEnergy in V.Ivanchenko model
>> 312 cost = 2.0*G4UniformRand() - 1.0;
>> 313 sint = std::sqrt(1. - cost*cost);
>> 314 phi = twopi * G4UniformRand();
>> 315 G4double del = std::sqrt(bindingEnergy *(bindingEnergy + 2.0*electron_mass_c2))
>> 316 / deltaMom;
>> 317 dirx += del* sint * std::cos(phi);
>> 318 diry += del* sint * std::sin(phi);
>> 319 dirz += del* cost;
>> 320
>> 321 // Find out new primary electron direction
>> 322 G4double finalPx = primaryMom*primaryDirection.x() - deltaMom*dirx;
>> 323 G4double finalPy = primaryMom*primaryDirection.y() - deltaMom*diry;
>> 324 G4double finalPz = primaryMom*primaryDirection.z() - deltaMom*dirz;
>> 325
>> 326 //Ok, ready to create the delta ray
>> 327 G4DynamicParticle* theDeltaRay = new G4DynamicParticle();
>> 328 theDeltaRay->SetKineticEnergy(energyDelta);
>> 329 G4double norm = 1.0/std::sqrt(dirx*dirx + diry*diry + dirz*dirz);
>> 330 dirx *= norm;
>> 331 diry *= norm;
>> 332 dirz *= norm;
>> 333 theDeltaRay->SetMomentumDirection(dirx, diry, dirz);
>> 334 theDeltaRay->SetDefinition(G4Electron::Electron());
>> 335 fvect->push_back(theDeltaRay);
299 336
300 //This is the amount of energy available for 337 //This is the amount of energy available for fluorescence
301 G4double theEnergyDeposit = bindingEnergy; 338 G4double theEnergyDeposit = bindingEnergy;
302 339
303 // fill ParticleChange 340 // fill ParticleChange
304 // changed energy and momentum of the actual 341 // changed energy and momentum of the actual particle
305 G4double finalKinEnergy = kineticEnergy - en 342 G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit;
306 if(finalKinEnergy < 0.0) 343 if(finalKinEnergy < 0.0)
307 { 344 {
308 theEnergyDeposit += finalKinEnergy; 345 theEnergyDeposit += finalKinEnergy;
309 finalKinEnergy = 0.0; 346 finalKinEnergy = 0.0;
310 } 347 }
311 else 348 else
312 { 349 {
313 fParticleChange->ProposeMomentumDirectio << 350 G4double normLocal = 1.0/std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
>> 351 finalPx *= normLocal;
>> 352 finalPy *= normLocal;
>> 353 finalPz *= normLocal;
>> 354 fParticleChange->ProposeMomentumDirection(finalPx, finalPy, finalPz);
314 } 355 }
315 fParticleChange->SetProposedKineticEnergy(fi 356 fParticleChange->SetProposedKineticEnergy(finalKinEnergy);
316 357
317 if (theEnergyDeposit < 0) 358 if (theEnergyDeposit < 0)
318 { 359 {
319 G4cout << "G4LivermoreIonisationModel: 360 G4cout << "G4LivermoreIonisationModel: Negative energy deposit: "
320 << theEnergyDeposit/eV << " eV" << G4en 361 << theEnergyDeposit/eV << " eV" << G4endl;
321 theEnergyDeposit = 0.0; 362 theEnergyDeposit = 0.0;
322 } 363 }
323 364
324 //Assign local energy deposit 365 //Assign local energy deposit
325 fParticleChange->ProposeLocalEnergyDeposit(t 366 fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit);
326 367
327 if (verboseLevel > 1) 368 if (verboseLevel > 1)
328 { 369 {
329 G4cout << "----------------------------- 370 G4cout << "-----------------------------------------------------------" << G4endl;
330 G4cout << "Energy balance from G4Livermo 371 G4cout << "Energy balance from G4LivermoreIonisation" << G4endl;
331 G4cout << "Incoming primary energy: " << 372 G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
332 G4cout << "----------------------------- 373 G4cout << "-----------------------------------------------------------" << G4endl;
333 G4cout << "Outgoing primary energy: " << 374 G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl;
334 G4cout << "Delta ray " << energyDelta/ke 375 G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl;
335 G4cout << "Fluorescence: " << (bindingEn 376 G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl;
336 G4cout << "Local energy deposit " << the 377 G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl;
337 G4cout << "Total final state: " << (fina 378 G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy)
338 << " keV" << G4endl; 379 << " keV" << G4endl;
339 G4cout << "----------------------------- 380 G4cout << "-----------------------------------------------------------" << G4endl;
340 } 381 }
341 return; 382 return;
342 } 383 }
343 384
344 //....oooOO0OOooo........oooOO0OOooo........oo 385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345 386
346 387