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

Last change on this file since 1316 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 15.1 KB
Line 
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.