| 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,v 1.13 2010/12/02 16:06:29 vnivanch Exp $
>> 27 // GEANT4 tag $Name: geant4-09-04 $
26 // 28 //
27 // Author: Luciano Pandola 29 // Author: Luciano Pandola
28 // on base of G4LowEnergyIonisation de <<
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 //
51 // 52 //
52 53
53 #include "G4LivermoreIonisationModel.hh" 54 #include "G4LivermoreIonisationModel.hh"
54 #include "G4PhysicalConstants.hh" <<
55 #include "G4SystemOfUnits.hh" <<
56 #include "G4ParticleDefinition.hh" 55 #include "G4ParticleDefinition.hh"
57 #include "G4MaterialCutsCouple.hh" 56 #include "G4MaterialCutsCouple.hh"
58 #include "G4ProductionCutsTable.hh" 57 #include "G4ProductionCutsTable.hh"
59 #include "G4DynamicParticle.hh" 58 #include "G4DynamicParticle.hh"
60 #include "G4Element.hh" 59 #include "G4Element.hh"
61 #include "G4ParticleChangeForLoss.hh" << 60 #include "G4AtomicTransitionManager.hh"
>> 61 #include "G4AtomicDeexcitation.hh"
>> 62 #include "G4AtomicShell.hh"
>> 63 #include "G4Gamma.hh"
62 #include "G4Electron.hh" 64 #include "G4Electron.hh"
63 #include "G4CrossSectionHandler.hh" 65 #include "G4CrossSectionHandler.hh"
>> 66 #include "G4ProcessManager.hh"
64 #include "G4VEMDataSet.hh" 67 #include "G4VEMDataSet.hh"
>> 68 #include "G4EMDataSet.hh"
65 #include "G4eIonisationCrossSectionHandler.hh" 69 #include "G4eIonisationCrossSectionHandler.hh"
66 #include "G4eIonisationSpectrum.hh" 70 #include "G4eIonisationSpectrum.hh"
67 #include "G4VEnergySpectrum.hh" << 71 #include "G4VDataSetAlgorithm.hh"
68 #include "G4SemiLogInterpolation.hh" 72 #include "G4SemiLogInterpolation.hh"
69 #include "G4AtomicTransitionManager.hh" << 73 #include "G4ShellVacancy.hh"
70 #include "G4AtomicShell.hh" << 74 #include "G4VDataSetAlgorithm.hh"
71 #include "G4DeltaAngle.hh" << 75 #include "G4LogLogInterpolation.hh"
>> 76 #include "G4CompositeEMDataSet.hh"
72 77
73 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 79
75 80
76 G4LivermoreIonisationModel::G4LivermoreIonisat 81 G4LivermoreIonisationModel::G4LivermoreIonisationModel(const G4ParticleDefinition*,
77 const G4String& nam) : << 82 const G4String& nam)
78 G4VEmModel(nam), fParticleChange(nullptr), << 83 :G4VEmModel(nam),isInitialised(false),crossSectionHandler(0),
79 crossSectionHandler(nullptr), energySpectrum << 84 energySpectrum(0),shellVacancy(0)
80 isInitialised(false) <<
81 { 85 {
82 fIntrinsicLowEnergyLimit = 12.*eV; << 86 fIntrinsicLowEnergyLimit = 10.0*eV;
83 fIntrinsicHighEnergyLimit = 100.0*GeV; 87 fIntrinsicHighEnergyLimit = 100.0*GeV;
84 << 88 fNBinEnergyLoss = 360;
>> 89 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
>> 90 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
>> 91 //
85 verboseLevel = 0; 92 verboseLevel = 0;
86 SetAngularDistribution(new G4DeltaAngle()); << 93 //By default: use deexcitation, not auger
87 << 94 SetDeexcitationFlag(true);
88 transitionManager = G4AtomicTransitionManage << 95 ActivateAuger(false);
>> 96 //
>> 97 //
>> 98 // Notice: the fluorescence along step is generated only if it is
>> 99 // set by the PROCESS (e.g. G4eIonisation) via the command
>> 100 // process->ActivateDeexcitation(true);
>> 101 //
89 } 102 }
90 103
91 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 105
93 G4LivermoreIonisationModel::~G4LivermoreIonisa 106 G4LivermoreIonisationModel::~G4LivermoreIonisationModel()
94 { 107 {
95 delete energySpectrum; << 108 if (energySpectrum) delete energySpectrum;
96 delete crossSectionHandler; << 109 if (crossSectionHandler) delete crossSectionHandler;
>> 110 if (shellVacancy) delete shellVacancy;
97 } 111 }
98 112
99 //....oooOO0OOooo........oooOO0OOooo........oo 113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 114
101 void G4LivermoreIonisationModel::Initialise(co 115 void G4LivermoreIonisationModel::Initialise(const G4ParticleDefinition* particle,
102 const G4DataVector& cuts) 116 const G4DataVector& cuts)
103 { 117 {
104 //Check that the Livermore Ionisation is NOT 118 //Check that the Livermore Ionisation is NOT attached to e+
105 if (particle != G4Electron::Electron()) 119 if (particle != G4Electron::Electron())
106 { 120 {
107 G4Exception("G4LivermoreIonisationModel: << 121 G4cout << "ERROR: Livermore Ionisation Model is applicable only to electrons" << G4endl;
108 "em0002",FatalException, << 122 G4cout << "It cannot be registered to " << particle->GetParticleName() << G4endl;
109 "Livermore Ionisation Model is applicabl << 123 G4Exception();
110 } 124 }
111 transitionManager->Initialise(); <<
112 125
113 //Read energy spectrum 126 //Read energy spectrum
114 if (energySpectrum) 127 if (energySpectrum)
115 { 128 {
116 delete energySpectrum; 129 delete energySpectrum;
117 energySpectrum = nullptr; << 130 energySpectrum = 0;
118 } 131 }
119 energySpectrum = new G4eIonisationSpectrum() 132 energySpectrum = new G4eIonisationSpectrum();
120 if (verboseLevel > 3) 133 if (verboseLevel > 3)
121 G4cout << "G4VEnergySpectrum is initialize 134 G4cout << "G4VEnergySpectrum is initialized" << G4endl;
122 135
123 //Initialize cross section handler 136 //Initialize cross section handler
124 if (crossSectionHandler) 137 if (crossSectionHandler)
125 { 138 {
126 delete crossSectionHandler; 139 delete crossSectionHandler;
127 crossSectionHandler = nullptr; << 140 crossSectionHandler = 0;
128 } 141 }
129 142
130 const size_t nbins = 20; <<
131 G4double emin = LowEnergyLimit(); <<
132 G4double emax = HighEnergyLimit(); <<
133 G4int ndec = G4int(std::log10(emax/emin) + 0 <<
134 if(ndec <= 0) { ndec = 1; } <<
135 <<
136 G4VDataSetAlgorithm* interpolation = new G4S 143 G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation();
137 crossSectionHandler = << 144 crossSectionHandler = new G4eIonisationCrossSectionHandler(energySpectrum,interpolation,
138 new G4eIonisationCrossSectionHandler(energ << 145 LowEnergyLimit(),HighEnergyLimit(),
139 emin,emax,nbins*ndec); << 146 fNBinEnergyLoss);
140 crossSectionHandler->Clear(); 147 crossSectionHandler->Clear();
141 crossSectionHandler->LoadShellData("ioni/ion 148 crossSectionHandler->LoadShellData("ioni/ion-ss-cs-");
142 //This is used to retrieve cross section val 149 //This is used to retrieve cross section values later on
143 G4VEMDataSet* emdata = 150 G4VEMDataSet* emdata =
144 crossSectionHandler->BuildMeanFreePathForM 151 crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
145 //The method BuildMeanFreePathForMaterials() 152 //The method BuildMeanFreePathForMaterials() is required here only to force
146 //the building of an internal table: the out 153 //the building of an internal table: the output pointer can be deleted
147 delete emdata; 154 delete emdata;
>> 155
>> 156 //Fluorescence data
>> 157 transitionManager = G4AtomicTransitionManager::Instance();
>> 158 if (shellVacancy) delete shellVacancy;
>> 159 shellVacancy = new G4ShellVacancy();
>> 160 InitialiseFluorescence();
148 161
149 if (verboseLevel > 0) 162 if (verboseLevel > 0)
150 { 163 {
151 G4cout << "Livermore Ionisation model is 164 G4cout << "Livermore Ionisation model is initialized " << G4endl
152 << "Energy range: " 165 << "Energy range: "
153 << LowEnergyLimit() / keV << " keV - " 166 << LowEnergyLimit() / keV << " keV - "
154 << HighEnergyLimit() / GeV << " GeV" 167 << HighEnergyLimit() / GeV << " GeV"
155 << G4endl; 168 << G4endl;
156 } 169 }
157 170
158 if (verboseLevel > 3) 171 if (verboseLevel > 3)
159 { 172 {
160 G4cout << "Cross section data: " << G4en 173 G4cout << "Cross section data: " << G4endl;
161 crossSectionHandler->PrintData(); 174 crossSectionHandler->PrintData();
162 G4cout << "Parameters: " << G4endl; 175 G4cout << "Parameters: " << G4endl;
163 energySpectrum->PrintData(); 176 energySpectrum->PrintData();
164 } 177 }
165 178
166 if(isInitialised) { return; } << 179 if(isInitialised) return;
167 fParticleChange = GetParticleChangeForLoss() 180 fParticleChange = GetParticleChangeForLoss();
168 isInitialised = true; 181 isInitialised = true;
169 } 182 }
170 183
171 //....oooOO0OOooo........oooOO0OOooo........oo 184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172 185
173 G4double << 186 G4double G4LivermoreIonisationModel::MinEnergyCut(const G4ParticleDefinition*,
174 G4LivermoreIonisationModel::ComputeCrossSectio << 187 const G4MaterialCutsCouple*)
175 const G4ParticleDe <<
176 G4double energy, <<
177 G4double Z, G4double, <<
178 G4double cutEnergy, <<
179 G4double) <<
180 { 188 {
181 G4int iZ = G4int(Z); << 189 return 250.*eV;
>> 190 }
>> 191
>> 192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 193
>> 194 G4double G4LivermoreIonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
>> 195 G4double energy,
>> 196 G4double Z, G4double,
>> 197 G4double cutEnergy,
>> 198 G4double)
>> 199 {
>> 200 G4int iZ = (G4int) Z;
182 if (!crossSectionHandler) 201 if (!crossSectionHandler)
183 { 202 {
184 G4Exception("G4LivermoreIonisationModel: << 203 G4cout << "G4LivermoreIonisationModel::ComputeCrossSectionPerAtom" << G4endl;
185 "em1007",FatalException, << 204 G4cout << "The cross section handler is not correctly initialized" << G4endl;
186 "The cross section handler is not correc << 205 G4Exception();
187 return 0; 206 return 0;
188 } 207 }
189 208
190 //The cut is already included in the crossSe 209 //The cut is already included in the crossSectionHandler
191 G4double cs = 210 G4double cs =
192 crossSectionHandler->GetCrossSectionAboveT << 211 crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,cutEnergy,iZ);
193 cutEnergy, <<
194 iZ); <<
195 212
196 if (verboseLevel > 1) 213 if (verboseLevel > 1)
197 { 214 {
198 G4cout << "G4LivermoreIonisationModel " 215 G4cout << "G4LivermoreIonisationModel " << G4endl;
199 G4cout << "Cross section for delta emiss << 216 G4cout << "Cross section for delta emission > " << cutEnergy/keV << " keV at " <<
200 << cutEnergy/keV << " keV at " << 217 energy/keV << " keV and Z = " << iZ << " --> " << cs/barn << " barn" << G4endl;
201 << energy/keV << " keV and Z = " << iZ <<
202 << cs/barn << " barn" << G4endl; <<
203 } 218 }
204 return cs; 219 return cs;
205 } 220 }
206 221
207 222
208 //....oooOO0OOooo........oooOO0OOooo........oo 223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
209 224
210 G4double << 225 G4double G4LivermoreIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
211 G4LivermoreIonisationModel::ComputeDEDXPerVolu << 226 const G4ParticleDefinition* ,
212 const G4ParticleDefinition*, << 227 G4double kineticEnergy,
213 G4double kineticEnergy, << 228 G4double cutEnergy)
214 G4double cutEnergy) <<
215 { 229 {
216 G4double sPower = 0.0; 230 G4double sPower = 0.0;
217 231
218 const G4ElementVector* theElementVector = ma 232 const G4ElementVector* theElementVector = material->GetElementVector();
219 size_t NumberOfElements = material->GetNumbe 233 size_t NumberOfElements = material->GetNumberOfElements() ;
220 const G4double* theAtomicNumDensityVector = 234 const G4double* theAtomicNumDensityVector =
221 material->GetAtomicNumDens 235 material->GetAtomicNumDensityVector();
222 236
223 // loop for elements in the material 237 // loop for elements in the material
224 for (size_t iel=0; iel<NumberOfElements; iel 238 for (size_t iel=0; iel<NumberOfElements; iel++ )
225 { 239 {
226 G4int iZ = (G4int)((*theElementVector)[i 240 G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
227 G4int nShells = transitionManager->Numbe 241 G4int nShells = transitionManager->NumberOfShells(iZ);
228 for (G4int n=0; n<nShells; n++) 242 for (G4int n=0; n<nShells; n++)
229 { 243 {
230 G4double e = energySpectrum->AverageEnergy 244 G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
231 kineticEnergy, n); 245 kineticEnergy, n);
232 G4double cs= crossSectionHandler->FindValu 246 G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n);
233 sPower += e * cs * theAtomicNumDensityVe 247 sPower += e * cs * theAtomicNumDensityVector[iel];
234 } 248 }
235 G4double esp = energySpectrum->Excitatio 249 G4double esp = energySpectrum->Excitation(iZ,kineticEnergy);
236 sPower += esp * theAtomicNumDensityVec 250 sPower += esp * theAtomicNumDensityVector[iel];
237 } 251 }
238 252
239 if (verboseLevel > 2) 253 if (verboseLevel > 2)
240 { 254 {
241 G4cout << "G4LivermoreIonisationModel " 255 G4cout << "G4LivermoreIonisationModel " << G4endl;
242 G4cout << "Stopping power < " << cutEner << 256 G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
243 << " keV at " << kineticEnergy/keV << " << 257 kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl;
244 << sPower/(keV/mm) << " keV/mm" << G4en <<
245 } 258 }
246 259
247 return sPower; 260 return sPower;
248 } 261 }
249 262
250 //....oooOO0OOooo........oooOO0OOooo........oo 263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
251 264
252 void G4LivermoreIonisationModel::SampleSeconda << 265 void G4LivermoreIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
253 std::vector<G << 266 const G4MaterialCutsCouple* couple,
254 const G4MaterialCutsCouple* couple, << 267 const G4DynamicParticle* aDynamicParticle,
255 const G4DynamicParticle* aDynamicPart << 268 G4double cutE,
256 G4double cutE, << 269 G4double maxE)
257 G4double maxE) <<
258 { 270 {
259 271
260 G4double kineticEnergy = aDynamicParticle->G 272 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
261 273
262 if (kineticEnergy <= fIntrinsicLowEnergyLimi 274 if (kineticEnergy <= fIntrinsicLowEnergyLimit)
263 { 275 {
264 fParticleChange->SetProposedKineticEnerg 276 fParticleChange->SetProposedKineticEnergy(0.);
265 fParticleChange->ProposeLocalEnergyDepos 277 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
266 return; << 278 return ;
267 } 279 }
268 280
269 // Select atom and shell 281 // Select atom and shell
270 G4int Z = crossSectionHandler->SelectRandomA 282 G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
271 G4int shellIndex = crossSectionHandler->Sele << 283 G4int shell = crossSectionHandler->SelectRandomShell(Z, kineticEnergy);
272 const G4AtomicShell* shell = transitionManag << 284 const G4AtomicShell* atomicShell =
273 G4double bindingEnergy = shell->BindingEnerg << 285 (G4AtomicTransitionManager::Instance())->Shell(Z, shell);
>> 286 G4double bindingEnergy = atomicShell->BindingEnergy();
>> 287 G4int shellId = atomicShell->ShellId();
274 288
275 // Sample delta energy using energy interval 289 // Sample delta energy using energy interval for delta-electrons
276 G4double energyMax = 290 G4double energyMax =
277 std::min(maxE,energySpectrum->MaxEnergyOfS 291 std::min(maxE,energySpectrum->MaxEnergyOfSecondaries(kineticEnergy));
278 G4double energyDelta = energySpectrum->Sampl 292 G4double energyDelta = energySpectrum->SampleEnergy(Z, cutE, energyMax,
279 kineticEnergy, shellIndex); << 293 kineticEnergy, shell);
280 294
281 if (energyDelta == 0.) //nothing happens 295 if (energyDelta == 0.) //nothing happens
282 { return; } << 296 return;
283 <<
284 const G4ParticleDefinition* electron = G4Ele <<
285 G4DynamicParticle* delta = new G4DynamicPart <<
286 GetAngularDistribution()->SampleDirectionF <<
287 Z, shellIndex, <<
288 <<
289 <<
290 <<
291 fvect->push_back(delta); <<
292 <<
293 // Change kinematics of primary particle <<
294 G4ThreeVector direction = aDynamicParticle-> <<
295 G4double totalMomentum = std::sqrt(kineticEn <<
296 297
297 G4ThreeVector finalP = totalMomentum*directi << 298 // Transform to shell potential
298 finalP = finalP.unit(); << 299 G4double deltaKinE = energyDelta + 2.0*bindingEnergy;
>> 300 G4double primaryKinE = kineticEnergy + 2.0*bindingEnergy;
>> 301
>> 302 // sampling of scattering angle neglecting atomic motion
>> 303 G4double deltaMom = std::sqrt(deltaKinE*(deltaKinE + 2.0*electron_mass_c2));
>> 304 G4double primaryMom = std::sqrt(primaryKinE*(primaryKinE + 2.0*electron_mass_c2));
>> 305
>> 306 G4double cost = deltaKinE * (primaryKinE + 2.0*electron_mass_c2)
>> 307 / (deltaMom * primaryMom);
>> 308 if (cost > 1.) cost = 1.;
>> 309 G4double sint = std::sqrt(1. - cost*cost);
>> 310 G4double phi = twopi * G4UniformRand();
>> 311 G4double dirx = sint * std::cos(phi);
>> 312 G4double diry = sint * std::sin(phi);
>> 313 G4double dirz = cost;
>> 314
>> 315 // Rotate to incident electron direction
>> 316 G4ThreeVector primaryDirection = aDynamicParticle->GetMomentumDirection();
>> 317 G4ThreeVector deltaDir(dirx,diry,dirz);
>> 318 deltaDir.rotateUz(primaryDirection);
>> 319 //Updated components
>> 320 dirx = deltaDir.x();
>> 321 diry = deltaDir.y();
>> 322 dirz = deltaDir.z();
>> 323
>> 324 // Take into account atomic motion del is relative momentum of the motion
>> 325 // kinetic energy of the motion == bindingEnergy in V.Ivanchenko model
>> 326 cost = 2.0*G4UniformRand() - 1.0;
>> 327 sint = std::sqrt(1. - cost*cost);
>> 328 phi = twopi * G4UniformRand();
>> 329 G4double del = std::sqrt(bindingEnergy *(bindingEnergy + 2.0*electron_mass_c2))
>> 330 / deltaMom;
>> 331 dirx += del* sint * std::cos(phi);
>> 332 diry += del* sint * std::sin(phi);
>> 333 dirz += del* cost;
>> 334
>> 335 // Find out new primary electron direction
>> 336 G4double finalPx = primaryMom*primaryDirection.x() - deltaMom*dirx;
>> 337 G4double finalPy = primaryMom*primaryDirection.y() - deltaMom*diry;
>> 338 G4double finalPz = primaryMom*primaryDirection.z() - deltaMom*dirz;
>> 339
>> 340 //Ok, ready to create the delta ray
>> 341 G4DynamicParticle* theDeltaRay = new G4DynamicParticle();
>> 342 theDeltaRay->SetKineticEnergy(energyDelta);
>> 343 G4double norm = 1.0/std::sqrt(dirx*dirx + diry*diry + dirz*dirz);
>> 344 dirx *= norm;
>> 345 diry *= norm;
>> 346 dirz *= norm;
>> 347 theDeltaRay->SetMomentumDirection(dirx, diry, dirz);
>> 348 theDeltaRay->SetDefinition(G4Electron::Electron());
>> 349 fvect->push_back(theDeltaRay);
299 350
300 //This is the amount of energy available for 351 //This is the amount of energy available for fluorescence
301 G4double theEnergyDeposit = bindingEnergy; 352 G4double theEnergyDeposit = bindingEnergy;
302 353
303 // fill ParticleChange 354 // fill ParticleChange
304 // changed energy and momentum of the actual 355 // changed energy and momentum of the actual particle
305 G4double finalKinEnergy = kineticEnergy - en 356 G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit;
306 if(finalKinEnergy < 0.0) 357 if(finalKinEnergy < 0.0)
307 { 358 {
308 theEnergyDeposit += finalKinEnergy; 359 theEnergyDeposit += finalKinEnergy;
309 finalKinEnergy = 0.0; 360 finalKinEnergy = 0.0;
310 } 361 }
311 else 362 else
312 { 363 {
313 fParticleChange->ProposeMomentumDirectio << 364 G4double norm = 1.0/std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
>> 365 finalPx *= norm;
>> 366 finalPy *= norm;
>> 367 finalPz *= norm;
>> 368 fParticleChange->ProposeMomentumDirection(finalPx, finalPy, finalPz);
314 } 369 }
315 fParticleChange->SetProposedKineticEnergy(fi 370 fParticleChange->SetProposedKineticEnergy(finalKinEnergy);
316 371
>> 372 // deexcitation may be active per G4Region
>> 373 if(DeexcitationFlag() && Z > 5)
>> 374 {
>> 375 G4ProductionCutsTable* theCoupleTable =
>> 376 G4ProductionCutsTable::GetProductionCutsTable();
>> 377 // Retrieve cuts for gammas
>> 378 G4double cutG = (*theCoupleTable->GetEnergyCutsVector(0))[couple->GetIndex()];
>> 379 if(theEnergyDeposit > cutG || theEnergyDeposit > cutE)
>> 380 {
>> 381 deexcitationManager.SetCutForSecondaryPhotons(cutG);
>> 382 deexcitationManager.SetCutForAugerElectrons(cutE);
>> 383 std::vector<G4DynamicParticle*>* secondaryVector =
>> 384 deexcitationManager.GenerateParticles(Z, shellId);
>> 385 G4DynamicParticle* aSecondary;
>> 386
>> 387 if (secondaryVector)
>> 388 {
>> 389 for (size_t i = 0; i<secondaryVector->size(); i++)
>> 390 {
>> 391 aSecondary = (*secondaryVector)[i];
>> 392 //Check if it is a valid secondary
>> 393 if (aSecondary)
>> 394 {
>> 395 G4double e = aSecondary->GetKineticEnergy();
>> 396 if (e < theEnergyDeposit)
>> 397 {
>> 398 theEnergyDeposit -= e;
>> 399 fvect->push_back(aSecondary);
>> 400 aSecondary = 0;
>> 401 (*secondaryVector)[i]=0;
>> 402 }
>> 403 else
>> 404 {
>> 405 delete aSecondary;
>> 406 (*secondaryVector)[i] = 0;
>> 407 }
>> 408 }
>> 409 }
>> 410 //secondaryVector = 0;
>> 411 delete secondaryVector;
>> 412 }
>> 413 }
>> 414 }
>> 415
317 if (theEnergyDeposit < 0) 416 if (theEnergyDeposit < 0)
318 { 417 {
319 G4cout << "G4LivermoreIonisationModel: 418 G4cout << "G4LivermoreIonisationModel: Negative energy deposit: "
320 << theEnergyDeposit/eV << " eV" << G4en 419 << theEnergyDeposit/eV << " eV" << G4endl;
321 theEnergyDeposit = 0.0; 420 theEnergyDeposit = 0.0;
322 } 421 }
323 422
324 //Assign local energy deposit 423 //Assign local energy deposit
325 fParticleChange->ProposeLocalEnergyDeposit(t 424 fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit);
326 425
327 if (verboseLevel > 1) 426 if (verboseLevel > 1)
328 { 427 {
329 G4cout << "----------------------------- 428 G4cout << "-----------------------------------------------------------" << G4endl;
330 G4cout << "Energy balance from G4Livermo 429 G4cout << "Energy balance from G4LivermoreIonisation" << G4endl;
331 G4cout << "Incoming primary energy: " << 430 G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
332 G4cout << "----------------------------- 431 G4cout << "-----------------------------------------------------------" << G4endl;
333 G4cout << "Outgoing primary energy: " << 432 G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl;
334 G4cout << "Delta ray " << energyDelta/ke 433 G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl;
335 G4cout << "Fluorescence: " << (bindingEn 434 G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl;
336 G4cout << "Local energy deposit " << the 435 G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl;
337 G4cout << "Total final state: " << (fina << 436 G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy+
338 << " keV" << G4endl; << 437 theEnergyDeposit)/keV << " keV" << G4endl;
339 G4cout << "----------------------------- 438 G4cout << "-----------------------------------------------------------" << G4endl;
340 } 439 }
341 return; 440 return;
342 } 441 }
343 442
344 //....oooOO0OOooo........oooOO0OOooo........oo 443 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 444
>> 445
>> 446 void G4LivermoreIonisationModel::SampleDeexcitationAlongStep(const G4Material* theMaterial,
>> 447 const G4Track& theTrack,
>> 448 G4double& eloss)
>> 449 {
>> 450 //No call if there is no deexcitation along step
>> 451 if (!DeexcitationFlag()) return;
>> 452
>> 453 //This method gets the energy loss calculated "Along the Step" and
>> 454 //(including fluctuations) and produces explicit fluorescence/Auger
>> 455 //secondaries. The eloss value is updated.
>> 456 G4double energyLossBefore = eloss;
>> 457
>> 458 if (verboseLevel > 2)
>> 459 {
>> 460 G4cout << "-----------------------------------------------------------" << G4endl;
>> 461 G4cout << " SampleDeexcitationAlongStep() from G4LivermoreIonisation" << G4endl;
>> 462 G4cout << "Energy loss along step before deexcitation : " << energyLossBefore/keV <<
>> 463 " keV" << G4endl;
>> 464 }
>> 465 G4double incidentEnergy = theTrack.GetDynamicParticle()->GetKineticEnergy();
>> 466
>> 467 G4ProductionCutsTable* theCoupleTable =
>> 468 G4ProductionCutsTable::GetProductionCutsTable();
>> 469 const G4MaterialCutsCouple* couple = theTrack.GetMaterialCutsCouple();
>> 470 size_t index = couple->GetIndex();
>> 471 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
>> 472 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
>> 473
>> 474
>> 475 std::vector<G4DynamicParticle*>* deexcitationProducts =
>> 476 new std::vector<G4DynamicParticle*>;
>> 477
>> 478 if(eloss > cute || eloss > cutg)
>> 479 {
>> 480 const G4AtomicTransitionManager* transitionManager =
>> 481 G4AtomicTransitionManager::Instance();
>> 482 deexcitationManager.SetCutForSecondaryPhotons(cutg);
>> 483 deexcitationManager.SetCutForAugerElectrons(cute);
>> 484
>> 485 size_t nElements = theMaterial->GetNumberOfElements();
>> 486 const G4ElementVector* theElementVector = theMaterial->GetElementVector();
>> 487
>> 488 std::vector<G4DynamicParticle*>* secVector = 0;
>> 489 G4DynamicParticle* aSecondary = 0;
>> 490 //G4ParticleDefinition* type = 0;
>> 491 G4double e;
>> 492 G4ThreeVector position;
>> 493 G4int shell, shellId;
>> 494
>> 495 // sample secondaries
>> 496 G4double eTot = 0.0;
>> 497 std::vector<G4int> n =
>> 498 shellVacancy->GenerateNumberOfIonisations(couple,
>> 499 incidentEnergy,eloss);
>> 500 for (size_t i=0; i<nElements; i++)
>> 501 {
>> 502 G4int Z = (G4int)((*theElementVector)[i]->GetZ());
>> 503 size_t nVacancies = n[i];
>> 504 G4double maxE = transitionManager->Shell(Z, 0)->BindingEnergy();
>> 505 if (nVacancies > 0 && Z > 5 && (maxE > cute || maxE > cutg))
>> 506 {
>> 507 for (size_t j=0; j<nVacancies; j++)
>> 508 {
>> 509 shell = crossSectionHandler->SelectRandomShell(Z, incidentEnergy);
>> 510 shellId = transitionManager->Shell(Z, shell)->ShellId();
>> 511 G4double maxEShell =
>> 512 transitionManager->Shell(Z, shell)->BindingEnergy();
>> 513 if (maxEShell > cute || maxEShell > cutg )
>> 514 {
>> 515 secVector = deexcitationManager.GenerateParticles(Z, shellId);
>> 516 if (secVector)
>> 517 {
>> 518 for (size_t l = 0; l<secVector->size(); l++) {
>> 519 aSecondary = (*secVector)[l];
>> 520 if (aSecondary)
>> 521 {
>> 522 e = aSecondary->GetKineticEnergy();
>> 523 G4double itsCut = cutg;
>> 524 if (aSecondary->GetParticleDefinition() == G4Electron::Electron())
>> 525 itsCut = cute;
>> 526 if ( eTot + e <= eloss && e > itsCut )
>> 527 {
>> 528 eTot += e;
>> 529 deexcitationProducts->push_back(aSecondary);
>> 530 }
>> 531 else
>> 532 {
>> 533 delete aSecondary;
>> 534 }
>> 535 }
>> 536 }
>> 537 delete secVector;
>> 538 }
>> 539 }
>> 540 }
>> 541 }
>> 542 }
>> 543 }
>> 544
>> 545 G4double energyLossInFluorescence = 0.0;
>> 546 size_t nSecondaries = deexcitationProducts->size();
>> 547 if (nSecondaries > 0)
>> 548 {
>> 549 //You may have already secondaries produced by SampleSubCutSecondaries()
>> 550 //at the process G4VEnergyLossProcess
>> 551 G4int secondariesBefore = fParticleChange->GetNumberOfSecondaries();
>> 552 fParticleChange->SetNumberOfSecondaries(nSecondaries+secondariesBefore);
>> 553 const G4StepPoint* preStep = theTrack.GetStep()->GetPreStepPoint();
>> 554 const G4StepPoint* postStep = theTrack.GetStep()->GetPostStepPoint();
>> 555 G4ThreeVector r = preStep->GetPosition();
>> 556 G4ThreeVector deltaR = postStep->GetPosition();
>> 557 deltaR -= r;
>> 558 G4double t = preStep->GetGlobalTime();
>> 559 G4double deltaT = postStep->GetGlobalTime();
>> 560 deltaT -= t;
>> 561 G4double time, q;
>> 562 G4ThreeVector position;
>> 563
>> 564 for (size_t i=0; i<nSecondaries; i++)
>> 565 {
>> 566 G4DynamicParticle* part = (*deexcitationProducts)[i];
>> 567 if (part)
>> 568 {
>> 569 G4double eSecondary = part->GetKineticEnergy();
>> 570 eloss -= eSecondary;
>> 571 if (eloss > 0.)
>> 572 {
>> 573 q = G4UniformRand();
>> 574 time = deltaT*q + t;
>> 575 position = deltaR*q;
>> 576 position += r;
>> 577 G4Track* newTrack = new G4Track(part, time, position);
>> 578 energyLossInFluorescence += eSecondary;
>> 579 pParticleChange->AddSecondary(newTrack);
>> 580 }
>> 581 else
>> 582 {
>> 583 eloss += eSecondary;
>> 584 delete part;
>> 585 part = 0;
>> 586 }
>> 587 }
>> 588 }
>> 589 }
>> 590 delete deexcitationProducts;
>> 591
>> 592 //Check and verbosities. Ensure energy conservation
>> 593 if (verboseLevel > 2)
>> 594 {
>> 595 G4cout << "Energy loss along step after deexcitation : " << eloss/keV <<
>> 596 " keV" << G4endl;
>> 597 }
>> 598 if (verboseLevel > 1)
>> 599 {
>> 600 G4cout << "------------------------------------------------------------------" << G4endl;
>> 601 G4cout << "Energy in fluorescence: " << energyLossInFluorescence/keV << " keV" << G4endl;
>> 602 G4cout << "Residual energy loss: " << eloss/keV << " keV " << G4endl;
>> 603 G4cout << "Total final: " << (energyLossInFluorescence+eloss)/keV << " keV" << G4endl;
>> 604 G4cout << "Total initial: " << energyLossBefore/keV << " keV" << G4endl;
>> 605 G4cout << "------------------------------------------------------------------" << G4endl;
>> 606 }
>> 607 if (verboseLevel > 0)
>> 608 {
>> 609 if (std::fabs(energyLossBefore-energyLossInFluorescence-eloss)>10*eV)
>> 610 {
>> 611 G4cout << "Found energy non-conservation at SampleDeexcitationAlongStep() " << G4endl;
>> 612 G4cout << "Energy in fluorescence: " << energyLossInFluorescence/keV << " keV" << G4endl;
>> 613 G4cout << "Residual energy loss: " << eloss/keV << " keV " << G4endl;
>> 614 G4cout << "Total final: " << (energyLossInFluorescence+eloss)/keV << " keV" << G4endl;
>> 615 G4cout << "Total initial: " << energyLossBefore/keV << " keV" << G4endl;
>> 616 }
>> 617 }
>> 618 }
>> 619
>> 620 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 621
>> 622 void G4LivermoreIonisationModel::InitialiseFluorescence()
>> 623 {
>> 624 G4DataVector* ksi = 0;
>> 625 G4DataVector* energyVector = 0;
>> 626 size_t binForFluo = fNBinEnergyLoss/10;
>> 627
>> 628 //Used to produce a log-spaced energy grid. To be deleted at the end.
>> 629 G4PhysicsLogVector* eVector = new G4PhysicsLogVector(LowEnergyLimit(),HighEnergyLimit(),
>> 630 binForFluo);
>> 631 const G4ProductionCutsTable* theCoupleTable=
>> 632 G4ProductionCutsTable::GetProductionCutsTable();
>> 633 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 634
>> 635 // Loop on couples
>> 636 for (size_t m=0; m<numOfCouples; m++)
>> 637 {
>> 638 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(m);
>> 639 const G4Material* material= couple->GetMaterial();
>> 640
>> 641 const G4ElementVector* theElementVector = material->GetElementVector();
>> 642 size_t NumberOfElements = material->GetNumberOfElements() ;
>> 643 const G4double* theAtomicNumDensityVector =
>> 644 material->GetAtomicNumDensityVector();
>> 645
>> 646 G4VDataSetAlgorithm* interp = new G4LogLogInterpolation();
>> 647 G4VEMDataSet* xsis = new G4CompositeEMDataSet(interp, 1., 1.);
>> 648 //loop on elements
>> 649 G4double energyCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[m];
>> 650 for (size_t iel=0; iel<NumberOfElements; iel++ )
>> 651 {
>> 652 G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
>> 653 energyVector = new G4DataVector();
>> 654 ksi = new G4DataVector();
>> 655 //Loop on energy
>> 656 for (size_t j = 0; j<binForFluo; j++)
>> 657 {
>> 658 G4double energy = eVector->GetLowEdgeEnergy(j);
>> 659 G4double cross = 0.;
>> 660 G4double eAverage= 0.;
>> 661 G4int nShells = transitionManager->NumberOfShells(iZ);
>> 662
>> 663 for (G4int n=0; n<nShells; n++)
>> 664 {
>> 665 G4double e = energySpectrum->AverageEnergy(iZ, 0.0,energyCut,
>> 666 energy, n);
>> 667 G4double pro = energySpectrum->Probability(iZ, 0.0,energyCut,
>> 668 energy, n);
>> 669 G4double cs= crossSectionHandler->FindValue(iZ, energy, n);
>> 670 eAverage += e * cs * theAtomicNumDensityVector[iel];
>> 671 cross += cs * pro * theAtomicNumDensityVector[iel];
>> 672 if(verboseLevel > 1)
>> 673 {
>> 674 G4cout << "Z= " << iZ
>> 675 << " shell= " << n
>> 676 << " E(keV)= " << energy/keV
>> 677 << " Eav(keV)= " << e/keV
>> 678 << " pro= " << pro
>> 679 << " cs= " << cs
>> 680 << G4endl;
>> 681 }
>> 682 }
>> 683
>> 684 G4double coeff = 0.0;
>> 685 if(eAverage > 0.)
>> 686 {
>> 687 coeff = cross/eAverage;
>> 688 eAverage /= cross;
>> 689 }
>> 690
>> 691 if(verboseLevel > 1)
>> 692 {
>> 693 G4cout << "Ksi Coefficient for Z= " << iZ
>> 694 << " E(keV)= " << energy/keV
>> 695 << " Eav(keV)= " << eAverage/keV
>> 696 << " coeff= " << coeff
>> 697 << G4endl;
>> 698 }
>> 699 energyVector->push_back(energy);
>> 700 ksi->push_back(coeff);
>> 701 }
>> 702 G4VEMDataSet* p = new G4EMDataSet(iZ,energyVector,ksi,interp,1.,1.);
>> 703 xsis->AddComponent(p);
>> 704 }
>> 705 if(verboseLevel>3) xsis->PrintData();
>> 706 shellVacancy->AddXsiTable(xsis);
>> 707 }
>> 708 if (eVector)
>> 709 delete eVector;
>> 710 }
>> 711
>> 712 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 713
>> 714 void G4LivermoreIonisationModel::ActivateAuger(G4bool val)
>> 715 {
>> 716 if (!DeexcitationFlag() && val)
>> 717 {
>> 718 G4cout << "WARNING - G4LivermoreIonisationModel" << G4endl;
>> 719 G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
>> 720 G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
>> 721 }
>> 722 deexcitationManager.ActivateAugerElectronProduction(val);
>> 723 if (verboseLevel > 1)
>> 724 G4cout << "Auger production set to " << val << G4endl;
>> 725 }
345 726
346 727