source: trunk/source/processes/electromagnetic/lowenergy/src/G4hParametrisedLossModel.cc@ 968

Last change on this file since 968 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 15.1 KB
RevLine 
[819]1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// -------------------------------------------------------------------
28//
29// GEANT4 Class file
30//
31//
32// File name: G4hParametrisedLossModel
33//
34// Author: V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
35//
36// Creation date: 20 July 2000
37//
38// Modifications:
39// 20/07/2000 V.Ivanchenko First implementation
40// 18/08/2000 V.Ivanchenko TRIM85 model is added
41// 03/10/2000 V.Ivanchenko CodeWizard clean up
42// 10/05/2001 V.Ivanchenko Clean up againist Linux compilation with -Wall
43// 30/12/2003 V.Ivanchenko SRIM2003 model is added
44// 07/05/2004 V.Ivanchenko Fix Graphite problem, add QAO model
45//
46// Class Description:
47//
48// Low energy protons/ions electronic stopping power parametrisation
49//
50// Class Description: End
51//
52// -------------------------------------------------------------------
53//
54//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56
57#include "G4hParametrisedLossModel.hh"
58#include "G4UnitsTable.hh"
59#include "globals.hh"
60#include "G4hZiegler1977p.hh"
61#include "G4hZiegler1977He.hh"
62#include "G4hZiegler1985p.hh"
63#include "G4hSRIM2000p.hh"
64//#include "G4hQAOModel.hh"
65#include "G4hICRU49p.hh"
66#include "G4hICRU49He.hh"
67#include "G4DynamicParticle.hh"
68#include "G4ParticleDefinition.hh"
69#include "G4ElementVector.hh"
70#include "G4Material.hh"
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
74G4hParametrisedLossModel::G4hParametrisedLossModel(const G4String& name)
75 :G4VLowEnergyModel(name), modelName(name)
76{
77 InitializeMe();
78}
79
80//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81
82void G4hParametrisedLossModel::InitializeMe()
83{
84 theZieglerFactor = eV*cm2*1.0e-15 ;
85
86 // Registration of parametrisation models
87 G4String blank = G4String(" ") ;
88 G4String zi77p = G4String("Ziegler1977p") ;
89 G4String zi77He = G4String("Ziegler1977He") ;
90 G4String ir49p = G4String("ICRU_R49p") ;
91 G4String ir49He = G4String("ICRU_R49He") ;
92 G4String zi85p = G4String("Ziegler1985p") ;
93 G4String zi00p = G4String("SRIM2000p") ;
94 G4String qao = G4String("QAO") ;
95 if(zi77p == modelName) {
96 eStopingPowerTable = new G4hZiegler1977p();
97 highEnergyLimit = 100.0*MeV;
98 lowEnergyLimit = 1.0*keV;
99
100 } else if(zi77He == modelName) {
101 eStopingPowerTable = new G4hZiegler1977He();
102 highEnergyLimit = 10.0*MeV/4.0;
103 lowEnergyLimit = 1.0*keV/4.0;
104
105 } else if(zi85p == modelName) {
106 eStopingPowerTable = new G4hZiegler1985p();
107 highEnergyLimit = 100.0*MeV;
108 lowEnergyLimit = 1.0*keV;
109
110 } else if(zi00p == modelName ) {
111 eStopingPowerTable = new G4hSRIM2000p();
112 highEnergyLimit = 100.0*MeV;
113 lowEnergyLimit = 1.0*keV;
114
115 } else if(ir49p == modelName || blank == modelName) {
116 eStopingPowerTable = new G4hICRU49p();
117 highEnergyLimit = 2.0*MeV;
118 lowEnergyLimit = 1.0*keV;
119
120 } else if(ir49He == modelName) {
121 eStopingPowerTable = new G4hICRU49He();
122 highEnergyLimit = 10.0*MeV/4.0;
123 lowEnergyLimit = 1.0*keV/4.0;
124 /*
125 } else if(qao == modelName) {
126 eStopingPowerTable = new G4hQAOModel();
127 highEnergyLimit = 2.0*MeV;
128 lowEnergyLimit = 5.0*keV;
129 */
130 } else {
131 eStopingPowerTable = new G4hICRU49p();
132 highEnergyLimit = 2.0*MeV;
133 lowEnergyLimit = 1.0*keV;
134 G4cout << "G4hParametrisedLossModel Warning: <" << modelName
135 << "> is unknown - default <"
136 << ir49p << ">" << " is used for Electronic Stopping"
137 << G4endl;
138 modelName = ir49p;
139 }
140 /*
141 G4cout << "G4hParametrisedLossModel: the model <"
142 << modelName << ">" << " is used for Electronic Stopping"
143 << G4endl;
144*/
145}
146
147//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148
149G4hParametrisedLossModel::~G4hParametrisedLossModel()
150{
151 delete eStopingPowerTable;
152}
153
154//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155
156G4double G4hParametrisedLossModel::TheValue(const G4DynamicParticle* particle,
157 const G4Material* material)
158{
159 G4double scaledEnergy = (particle->GetKineticEnergy())
160 * proton_mass_c2/(particle->GetMass());
161 G4double factor = theZieglerFactor;
162 if (scaledEnergy < lowEnergyLimit) {
163 if (modelName != "QAO") factor *= std::sqrt(scaledEnergy/lowEnergyLimit);
164 scaledEnergy = lowEnergyLimit;
165 }
166 G4double eloss = StoppingPower(material,scaledEnergy) * factor;
167
168 return eloss;
169}
170
171//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172
173G4double G4hParametrisedLossModel::TheValue(const G4ParticleDefinition* aParticle,
174 const G4Material* material,
175 G4double kineticEnergy)
176{
177 G4double scaledEnergy = kineticEnergy
178 * proton_mass_c2/(aParticle->GetPDGMass());
179
180 G4double factor = theZieglerFactor;
181 if (scaledEnergy < lowEnergyLimit) {
182 if (modelName != "QAO") factor *= std::sqrt(scaledEnergy/lowEnergyLimit);
183 scaledEnergy = lowEnergyLimit;
184 }
185 G4double eloss = StoppingPower(material,scaledEnergy) * factor;
186
187 return eloss;
188}
189
190//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
191
192G4double G4hParametrisedLossModel::LowEnergyLimit(const G4ParticleDefinition* ,
193 const G4Material*) const
194{
195 return lowEnergyLimit;
196}
197
198//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
199
200G4double G4hParametrisedLossModel::HighEnergyLimit(const G4ParticleDefinition* ,
201 const G4Material*) const
202{
203 return highEnergyLimit;
204}
205//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206
207G4double G4hParametrisedLossModel::LowEnergyLimit(const G4ParticleDefinition* ) const
208{
209 return lowEnergyLimit;
210}
211
212//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
213
214G4double G4hParametrisedLossModel::HighEnergyLimit(const G4ParticleDefinition* ) const
215{
216 return highEnergyLimit;
217}
218
219//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
220
221G4bool G4hParametrisedLossModel::IsInCharge(const G4DynamicParticle* ,
222 const G4Material*) const
223{
224 return true;
225}
226
227//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228
229G4bool G4hParametrisedLossModel::IsInCharge(const G4ParticleDefinition* ,
230 const G4Material*) const
231{
232 return true;
233}
234
235//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236
237G4double G4hParametrisedLossModel::StoppingPower(const G4Material* material,
238 G4double kineticEnergy)
239{
240 G4double eloss = 0.0;
241
242 const G4int numberOfElements = material->GetNumberOfElements() ;
243 const G4double* theAtomicNumDensityVector =
244 material->GetAtomicNumDensityVector() ;
245
246
247 // compound material with parametrisation
248 if( (eStopingPowerTable->HasMaterial(material)) ) {
249
250 eloss = eStopingPowerTable->StoppingPower(material, kineticEnergy);
251 if ("QAO" != modelName) {
252 eloss *= material->GetTotNbOfAtomsPerVolume();
253 if(1 < numberOfElements) {
254 G4int nAtoms = 0;
255
256 const G4int* theAtomsVector = material->GetAtomsVector();
257 for (G4int iel=0; iel<numberOfElements; iel++) {
258 nAtoms += theAtomsVector[iel];
259 }
260 eloss /= nAtoms;
261 }
262 }
263
264 // pure material
265 } else if(1 == numberOfElements) {
266
267 G4double z = material->GetZ();
268 eloss = (eStopingPowerTable->ElectronicStoppingPower(z, kineticEnergy))
269 * (material->GetTotNbOfAtomsPerVolume()) ;
270
271 // Experimental data exist only for kinetic energy 125 keV
272 } else if( MolecIsInZiegler1988(material)) {
273
274 // Cycle over elements - calculation based on Bragg's rule
275 G4double eloss125 = 0.0 ;
276 const G4ElementVector* theElementVector =
277 material->GetElementVector() ;
278
279
280 // loop for the elements in the material
281 for (G4int i=0; i<numberOfElements; i++) {
282 const G4Element* element = (*theElementVector)[i] ;
283 G4double z = element->GetZ() ;
284 eloss +=(eStopingPowerTable->ElectronicStoppingPower(z,kineticEnergy))
285 * theAtomicNumDensityVector[i] ;
286 eloss125 +=(eStopingPowerTable->ElectronicStoppingPower(z,125.0*keV))
287 * theAtomicNumDensityVector[i] ;
288 }
289
290 // Chemical factor is taken into account
291 eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
292
293 // Brugg's rule calculation
294 } else {
295 const G4ElementVector* theElementVector =
296 material->GetElementVector() ;
297
298 // loop for the elements in the material
299 for (G4int i=0; i<numberOfElements; i++)
300 {
301 const G4Element* element = (*theElementVector)[i] ;
302 G4double z = element->GetZ() ;
303 eloss += (eStopingPowerTable->ElectronicStoppingPower(z,kineticEnergy))
304 * theAtomicNumDensityVector[i];
305 }
306 }
307 return eloss;
308}
309
310//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
311
312G4bool G4hParametrisedLossModel::MolecIsInZiegler1988(
313 const G4Material* material)
314{
315 // The list of molecules from
316 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
317 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
318
319 G4String myFormula = G4String(" ") ;
320 const G4String chFormula = material->GetChemicalFormula() ;
321 if (myFormula == chFormula ) return false ;
322
323 // There are no evidence for difference of stopping power depended on
324 // phase of the compound except for water. The stopping power of the
325 // water in gas phase can be predicted using Bragg's rule.
326 //
327 // No chemical factor for water-gas
328
329 myFormula = G4String("H_2O") ;
330 const G4State theState = material->GetState() ;
331 if( theState == kStateGas && myFormula == chFormula) return false ;
332
333 const size_t numberOfMolecula = 53 ;
334
335 // The coffecient from Table.4 of Ziegler & Manoyan
336 const G4double HeEff = 2.8735 ;
337
338 static G4String name[numberOfMolecula] = {
339 "H_2O", "C_2H_4O", "C_3H_6O", "C_2H_2", "C_H_3OH",
340 "C_2H_5OH", "C_3H_7OH", "C_3H_4", "NH_3", "C_14H_10",
341 "C_6H_6", "C_4H_10", "C_4H_6", "C_4H_8O", "CCl_4",
342 "CF_4", "C_6H_8", "C_6H_12", "C_6H_10O", "C_6H_10",
343 "C_8H_16", "C_5H_10", "C_5H_8", "C_3H_6-Cyclopropane","C_2H_4F_2",
344 "C_2H_2F_2", "C_4H_8O_2", "C_2H_6", "C_2F_6", "C_2H_6O",
345 "C_3H_6O", "C_4H_10O", "C_2H_4", "C_2H_4O", "C_2H_4S",
346 "SH_2", "CH_4", "CCLF_3", "CCl_2F_2", "CHCl_2F",
347 "(CH_3)_2S", "N_2O", "C_5H_10O", "C_8H_6", "(CH_2)_N",
348 "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8", "C_3H_6-Propylene", "C_3H_6O",
349 "C_3H_6S", "C_4H_4S", "C_7H_8"
350 } ;
351
352 static G4double expStopping[numberOfMolecula] = {
353 66.1, 190.4, 258.7, 42.2, 141.5,
354 210.9, 279.6, 198.8, 31.0, 267.5,
355 122.8, 311.4, 260.3, 328.9, 391.3,
356 206.6, 374.0, 422.0, 432.0, 398.0,
357 554.0, 353.0, 326.0, 74.6, 220.5,
358 197.4, 362.0, 170.0, 330.5, 211.3,
359 262.3, 349.6, 51.3, 187.0, 236.9,
360 121.9, 35.8, 247.0, 292.6, 268.0,
361 262.3, 49.0, 398.9, 444.0, 22.91,
362 68.0, 155.0, 84.0, 74.2, 254.7,
363 306.8, 324.4, 420.0
364 } ;
365
366 static G4double expCharge[numberOfMolecula] = {
367 HeEff, HeEff, HeEff, 1.0, HeEff,
368 HeEff, HeEff, HeEff, 1.0, 1.0,
369 1.0, HeEff, HeEff, HeEff, HeEff,
370 HeEff, HeEff, HeEff, HeEff, HeEff,
371 HeEff, HeEff, HeEff, 1.0, HeEff,
372 HeEff, HeEff, HeEff, HeEff, HeEff,
373 HeEff, HeEff, 1.0, HeEff, HeEff,
374 HeEff, 1.0, HeEff, HeEff, HeEff,
375 HeEff, 1.0, HeEff, HeEff, 1.0,
376 1.0, 1.0, 1.0, 1.0, HeEff,
377 HeEff, HeEff, HeEff
378 } ;
379
380 static G4double numberOfAtomsPerMolecula[numberOfMolecula] = {
381 3.0, 7.0, 10.0, 4.0, 6.0,
382 9.0, 12.0, 7.0, 4.0, 24.0,
383 12.0, 14.0, 10.0, 13.0, 5.0,
384 5.0, 14.0, 18.0, 17.0, 17.0,
385 24.0, 15.0, 13.0, 9.0, 8.0,
386 6.0, 14.0, 8.0, 8.0, 9.0,
387 10.0, 15.0, 6.0, 7.0, 7.0,
388 3.0, 5.0, 5.0, 5.0, 5.0,
389 9.0, 3.0, 16.0, 14.0, 3.0,
390 9.0, 16.0, 11.0, 9.0, 10.0,
391 10.0, 9.0, 15.0
392 } ;
393
394 // Search for the compaund in the table
395 for (size_t i=0; i<numberOfMolecula; i++)
396 {
397 if(chFormula == name[i]) {
398 G4double exp125 = expStopping[i] *
399 (material->GetTotNbOfAtomsPerVolume()) /
400 (expCharge[i] * numberOfAtomsPerMolecula[i]) ;
401 SetExpStopPower125(exp125) ;
402 return true ;
403 }
404 }
405
406 return false ;
407}
408
409//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
410
411G4double G4hParametrisedLossModel::ChemicalFactor(
412 G4double kineticEnergy, G4double eloss125) const
413{
414 // Approximation of Chemical Factor according to
415 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
416 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
417
418 G4double gamma = 1.0 + kineticEnergy/proton_mass_c2 ;
419 G4double gamma25 = 1.0 + 25.0*keV /proton_mass_c2 ;
420 G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2 ;
421 G4double beta = std::sqrt(1.0 - 1.0/(gamma*gamma)) ;
422 G4double beta25 = std::sqrt(1.0 - 1.0/(gamma25*gamma25)) ;
423 G4double beta125 = std::sqrt(1.0 - 1.0/(gamma125*gamma125)) ;
424
425 G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) *
426 (1.0 + std::exp( 1.48 * ( beta125/beta25 - 7.0 ) ) ) /
427 (1.0 + std::exp( 1.48 * ( beta/beta25 - 7.0 ) ) ) ;
428
429 return factor ;
430}
431
432//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.