00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054 #include "G4hNuclearStoppingModel.hh"
00055
00056 #include "globals.hh"
00057 #include "G4PhysicalConstants.hh"
00058 #include "G4SystemOfUnits.hh"
00059 #include "G4UnitsTable.hh"
00060 #include "G4hZiegler1977Nuclear.hh"
00061 #include "G4hZiegler1985Nuclear.hh"
00062 #include "G4hICRU49Nuclear.hh"
00063 #include "G4DynamicParticle.hh"
00064 #include "G4ParticleDefinition.hh"
00065 #include "G4ElementVector.hh"
00066 #include "G4Material.hh"
00067
00068
00069
00070 G4hNuclearStoppingModel::G4hNuclearStoppingModel(const G4String& name)
00071 :G4VLowEnergyModel(name), modelName(name)
00072 {
00073 InitializeMe() ;
00074 }
00075
00076
00077
00078 void G4hNuclearStoppingModel::InitializeMe()
00079 {
00080
00081 highEnergyLimit = 100.*MeV ;
00082 lowEnergyLimit = 1.*eV ;
00083 factorPDG2AMU = 1.007276/proton_mass_c2 ;
00084 theZieglerFactor= eV*cm2*1.0e-15 ;
00085
00086
00087 G4String blank = G4String(" ") ;
00088 G4String zi77 = G4String("Ziegler1977") ;
00089 G4String ir49 = G4String("ICRU_R49") ;
00090 G4String zi85 = G4String("Ziegler1985") ;
00091 if(zi77 == modelName) {
00092 nStopingPowerTable = new G4hZiegler1977Nuclear();
00093
00094 } else if(ir49 == modelName || blank == modelName) {
00095 nStopingPowerTable = new G4hICRU49Nuclear();
00096
00097 } else if(zi85 == modelName) {
00098 nStopingPowerTable = new G4hZiegler1985Nuclear();
00099
00100 } else {
00101 G4cout <<
00102 "G4hLowEnergyIonisation warning: There is no table with the modelName <"
00103 << modelName << ">"
00104 << " for nuclear stopping, <ICRU_R49> is applied "
00105 << G4endl;
00106 nStopingPowerTable = new G4hICRU49Nuclear();
00107 }
00108
00109
00110
00111 nStopingPowerTable->SetNuclearStoppingFluctuationsOff();
00112 }
00113
00114
00115
00116 G4hNuclearStoppingModel::~G4hNuclearStoppingModel()
00117 {
00118 delete nStopingPowerTable;
00119 }
00120
00121
00122
00123 G4double G4hNuclearStoppingModel::TheValue(
00124 const G4DynamicParticle* particle,
00125 const G4Material* material)
00126 {
00127
00128 G4double energy = particle->GetKineticEnergy() ;
00129 G4double z1 = std::abs((particle->GetCharge())/eplus) ;
00130 G4double m1 = (particle->GetMass())*factorPDG2AMU ;
00131
00132 G4double nloss = StoppingPower(material, energy, z1, m1) * theZieglerFactor;
00133
00134 return nloss;
00135 }
00136
00137
00138
00139 G4double G4hNuclearStoppingModel::TheValue(
00140 const G4ParticleDefinition* aParticle,
00141 const G4Material* material,
00142 G4double kineticEnergy)
00143 {
00144
00145 G4double z1 = std::abs((aParticle->GetPDGCharge())/eplus) ;
00146 G4double m1 = (aParticle->GetPDGMass())*factorPDG2AMU ;
00147
00148 G4double nloss = StoppingPower(material, kineticEnergy, z1, m1)
00149 * theZieglerFactor;
00150
00151 return nloss;
00152 }
00153
00154
00155
00156 G4double G4hNuclearStoppingModel::StoppingPower(
00157 const G4Material* material,
00158 G4double kineticEnergy,
00159 G4double z1, G4double m1) const
00160 {
00161
00162 G4int NumberOfElements = material->GetNumberOfElements() ;
00163 if(0 == NumberOfElements) return 0.0 ;
00164
00165 const G4ElementVector* theElementVector =
00166 material->GetElementVector() ;
00167 const G4double* theAtomicNumDensityVector =
00168 material->GetAtomicNumDensityVector() ;
00169
00170
00171
00172 G4double nloss = 0.0;
00173
00174 for (G4int iel=0; iel<NumberOfElements; iel++) {
00175 const G4Element* element = (*theElementVector)[iel] ;
00176 G4double z2 = element->GetZ();
00177 G4double m2Local = element->GetA()*mole/g ;
00178 nloss += (nStopingPowerTable->
00179 NuclearStoppingPower(kineticEnergy, z1, z2, m1, m2Local))
00180 * theAtomicNumDensityVector[iel] ;
00181 }
00182
00183 return nloss;
00184 }
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194