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