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