Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: G4BraggModel.cc 74790 2013-10-22 07:31:37Z gcosmo $ 26 // 27 // 27 // ------------------------------------------- 28 // ------------------------------------------------------------------- 28 // 29 // 29 // GEANT4 Class file 30 // GEANT4 Class file 30 // 31 // 31 // 32 // 32 // File name: G4BraggModel 33 // File name: G4BraggModel 33 // 34 // 34 // Author: Vladimir Ivanchenko 35 // Author: Vladimir Ivanchenko 35 // 36 // 36 // Creation date: 03.01.2002 37 // Creation date: 03.01.2002 37 // 38 // 38 // Modifications: 39 // Modifications: 39 // 40 // 40 // 04-12-02 Fix problem of G4DynamicParticle c 41 // 04-12-02 Fix problem of G4DynamicParticle constructor (V.Ivanchenko) 41 // 23-12-02 Change interface in order to move 42 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko) 42 // 27-01-03 Make models region aware (V.Ivanch 43 // 27-01-03 Make models region aware (V.Ivanchenko) 43 // 13-02-03 Add name (V.Ivanchenko) 44 // 13-02-03 Add name (V.Ivanchenko) 44 // 04-06-03 Fix compilation warnings (V.Ivanch 45 // 04-06-03 Fix compilation warnings (V.Ivanchenko) 45 // 12-09-04 Add lowestKinEnergy and change ord 46 // 12-09-04 Add lowestKinEnergy and change order of if in DEDX method (VI) 46 // 11-04-05 Major optimisation of internal int 47 // 11-04-05 Major optimisation of internal interfaces (V.Ivantchenko) 47 // 16-06-05 Fix problem of chemical formula (V 48 // 16-06-05 Fix problem of chemical formula (V.Ivantchenko) 48 // 15-02-06 ComputeCrossSectionPerElectron, Co 49 // 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma) 49 // 25-04-06 Add stopping data from PSTAR (V.Iv 50 // 25-04-06 Add stopping data from PSTAR (V.Ivanchenko) 50 // 12-08-08 Added methods GetParticleCharge, G 51 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio, 51 // CorrectionsAlongStep needed for io 52 // CorrectionsAlongStep needed for ions(V.Ivanchenko) 52 53 53 // Class Description: 54 // Class Description: 54 // 55 // 55 // Implementation of energy loss and delta-ele 56 // Implementation of energy loss and delta-electron production by 56 // slow charged heavy particles 57 // slow charged heavy particles 57 58 58 // ------------------------------------------- 59 // ------------------------------------------------------------------- 59 // 60 // 60 61 61 62 62 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 63 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 64 65 65 #include "G4BraggModel.hh" 66 #include "G4BraggModel.hh" 66 #include "G4PhysicalConstants.hh" 67 #include "G4PhysicalConstants.hh" 67 #include "G4SystemOfUnits.hh" 68 #include "G4SystemOfUnits.hh" 68 #include "Randomize.hh" 69 #include "Randomize.hh" 69 #include "G4Electron.hh" 70 #include "G4Electron.hh" 70 #include "G4ParticleChangeForLoss.hh" 71 #include "G4ParticleChangeForLoss.hh" 71 #include "G4LossTableManager.hh" 72 #include "G4LossTableManager.hh" 72 #include "G4EmCorrections.hh" 73 #include "G4EmCorrections.hh" 73 #include "G4EmParameters.hh" << 74 #include "G4DeltaAngle.hh" 74 #include "G4DeltaAngle.hh" 75 #include "G4ICRU90StoppingData.hh" << 76 #include "G4NistManager.hh" << 77 #include "G4Log.hh" 75 #include "G4Log.hh" 78 #include "G4Exp.hh" 76 #include "G4Exp.hh" 79 #include "G4AutoLock.hh" << 80 77 81 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 82 79 83 G4ICRU90StoppingData* G4BraggModel::fICRU90 = << 80 using namespace std; 84 G4PSTARStopping* G4BraggModel::fPSTAR = nullpt << 85 81 86 namespace << 82 static const G4double invLog10 = 1.0/G4Log(10.); 87 { << 88 G4Mutex ionMutex = G4MUTEX_INITIALIZER; << 89 } << 90 83 91 G4BraggModel::G4BraggModel(const G4ParticleDef 84 G4BraggModel::G4BraggModel(const G4ParticleDefinition* p, const G4String& nam) 92 : G4VEmModel(nam) << 85 : G4VEmModel(nam), >> 86 particle(0), >> 87 currentMaterial(0), >> 88 protonMassAMU(1.007276), >> 89 iMolecula(-1), >> 90 iPSTAR(-1), >> 91 isIon(false), >> 92 isInitialised(false) 93 { 93 { 94 SetHighEnergyLimit(2.0*CLHEP::MeV); << 94 fParticleChange = 0; >> 95 SetHighEnergyLimit(2.0*MeV); 95 96 96 lowestKinEnergy = 0.25*CLHEP::keV; << 97 lowestKinEnergy = 1.0*keV; 97 theZieglerFactor = CLHEP::eV*CLHEP::cm2*1.0e << 98 theZieglerFactor = eV*cm2*1.0e-15; 98 theElectron = G4Electron::Electron(); 99 theElectron = G4Electron::Electron(); 99 expStopPower125 = 0.0; 100 expStopPower125 = 0.0; 100 101 101 corr = G4LossTableManager::Instance()->EmCor 102 corr = G4LossTableManager::Instance()->EmCorrections(); 102 if(nullptr != p) { SetParticle(p); } << 103 if(p) { SetParticle(p); } 103 else { SetParticle(theElectron); } << 104 else { SetParticle(theElectron); } 104 } 105 } 105 106 106 //....oooOO0OOooo........oooOO0OOooo........oo 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 107 108 108 G4BraggModel::~G4BraggModel() 109 G4BraggModel::~G4BraggModel() 109 { << 110 {} 110 if(isFirst) { << 111 delete fPSTAR; << 112 fPSTAR = nullptr; << 113 } << 114 } << 115 111 116 //....oooOO0OOooo........oooOO0OOooo........oo 112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 117 113 118 void G4BraggModel::Initialise(const G4Particle 114 void G4BraggModel::Initialise(const G4ParticleDefinition* p, 119 const G4DataVect 115 const G4DataVector&) 120 { 116 { 121 if(p != particle) { SetParticle(p); } 117 if(p != particle) { SetParticle(p); } 122 118 123 // always false before the run 119 // always false before the run 124 SetDeexcitationFlag(false); 120 SetDeexcitationFlag(false); 125 121 126 // initialise data only once << 122 if(!isInitialised) { 127 if(nullptr == fPSTAR) { << 123 isInitialised = true; 128 G4AutoLock l(&ionMutex); << 129 if(nullptr == fPSTAR) { << 130 isFirst = true; << 131 fPSTAR = new G4PSTARStopping(); << 132 if(G4EmParameters::Instance()->UseICRU90 << 133 fICRU90 = G4NistManager::Instance()->GetICRU << 134 } << 135 } << 136 l.unlock(); << 137 } << 138 if(isFirst) { << 139 if(nullptr != fICRU90) { fICRU90->Initiali << 140 fPSTAR->Initialise(); << 141 } << 142 << 143 if(nullptr == fParticleChange) { << 144 124 145 if(UseAngularGeneratorFlag() && !GetAngula 125 if(UseAngularGeneratorFlag() && !GetAngularDistribution()) { 146 SetAngularDistribution(new G4DeltaAngle( 126 SetAngularDistribution(new G4DeltaAngle()); 147 } 127 } 148 G4String pname = particle->GetParticleName 128 G4String pname = particle->GetParticleName(); 149 if(particle->GetParticleType() == "nucleus 129 if(particle->GetParticleType() == "nucleus" && 150 pname != "deuteron" && pname != "triton 130 pname != "deuteron" && pname != "triton" && 151 pname != "alpha+" && pname != "helium 131 pname != "alpha+" && pname != "helium" && 152 pname != "hydrogen") { isIon = true; } 132 pname != "hydrogen") { isIon = true; } 153 133 154 fParticleChange = GetParticleChangeForLoss 134 fParticleChange = GetParticleChangeForLoss(); 155 } 135 } 156 } 136 } 157 137 158 //....oooOO0OOooo........oooOO0OOooo........oo 138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 159 139 160 void G4BraggModel::SetParticle(const G4Particl << 161 { << 162 particle = p; << 163 mass = particle->GetPDGMass(); << 164 spin = particle->GetPDGSpin(); << 165 G4double q = particle->GetPDGCharge()/CLHEP: << 166 chargeSquare = q*q; << 167 massRate = mass/CLHEP::proton_mass_c2; << 168 ratio = CLHEP::electron_mass_c2/mass; << 169 } << 170 << 171 //....oooOO0OOooo........oooOO0OOooo........oo << 172 << 173 G4double G4BraggModel::GetChargeSquareRatio(co 140 G4double G4BraggModel::GetChargeSquareRatio(const G4ParticleDefinition* p, 174 co << 141 const G4Material* mat, 175 G4 << 142 G4double kineticEnergy) 176 { 143 { 177 // this method is called only for ions 144 // this method is called only for ions 178 chargeSquare = corr->EffectiveChargeSquareRa << 145 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy); 179 return chargeSquare; << 146 GetModelOfFluctuations()->SetParticleAndCharge(p, q2); 180 } << 147 return q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy); 181 << 182 //....oooOO0OOooo........oooOO0OOooo........oo << 183 << 184 G4double G4BraggModel::MinEnergyCut(const G4Pa << 185 const G4Ma << 186 { << 187 return couple->GetMaterial()->GetIonisation( << 188 } 148 } 189 149 190 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 191 151 192 G4double G4BraggModel::GetParticleCharge(const 152 G4double G4BraggModel::GetParticleCharge(const G4ParticleDefinition* p, 193 const << 153 const G4Material* mat, 194 G4dou << 154 G4double kineticEnergy) 195 { 155 { 196 // this method is called only for ions, so n 156 // this method is called only for ions, so no check if it is an ion 197 return corr->GetParticleCharge(p,mat,kinetic 157 return corr->GetParticleCharge(p,mat,kineticEnergy); 198 } 158 } 199 159 200 //....oooOO0OOooo........oooOO0OOooo........oo 160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 201 161 202 G4double G4BraggModel::ComputeCrossSectionPerE 162 G4double G4BraggModel::ComputeCrossSectionPerElectron( 203 con 163 const G4ParticleDefinition* p, 204 164 G4double kineticEnergy, 205 << 165 G4double cutEnergy, 206 166 G4double maxKinEnergy) 207 { 167 { 208 G4double cross = 0.0; << 168 G4double cross = 0.0; 209 const G4double tmax = MaxSecondaryEnerg << 169 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy); 210 const G4double maxEnergy = std::min(tmax, ma << 170 G4double maxEnergy = std::min(tmax,maxKinEnergy); 211 const G4double cutEnergy = std::max(cut, low << 212 if(cutEnergy < maxEnergy) { 171 if(cutEnergy < maxEnergy) { 213 172 214 const G4double energy = kineticEnergy + m << 173 G4double energy = kineticEnergy + mass; 215 const G4double energy2 = energy*energy; << 174 G4double energy2 = energy*energy; 216 const G4double beta2 = kineticEnergy*(ki << 175 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2; 217 cross = (maxEnergy - cutEnergy)/(cutEnergy << 176 cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*G4Log(maxEnergy/cutEnergy)/tmax; 218 - beta2*G4Log(maxEnergy/cutEnergy)/tmax; << 219 177 220 if( 0.0 < spin ) { cross += 0.5*(maxEnergy << 178 cross *= twopi_mc2_rcl2*chargeSquare/beta2; 221 << 222 cross *= CLHEP::twopi_mc2_rcl2*chargeSquar << 223 } 179 } 224 // G4cout << "BR: e= " << kineticEnergy << << 180 // G4cout << "BR: e= " << kineticEnergy << " tmin= " << cutEnergy 225 // << " tmax= " << tmax << " cross= << 181 // << " tmax= " << tmax << " cross= " << cross << G4endl; >> 182 226 return cross; 183 return cross; 227 } 184 } 228 185 229 //....oooOO0OOooo........oooOO0OOooo........oo 186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 230 187 231 G4double << 188 G4double G4BraggModel::ComputeCrossSectionPerAtom( 232 G4BraggModel::ComputeCrossSectionPerAtom(const << 189 const G4ParticleDefinition* p, 233 G4dou << 190 G4double kineticEnergy, 234 G4dou << 191 G4double Z, G4double, 235 G4dou << 192 G4double cutEnergy, 236 G4dou << 193 G4double maxEnergy) 237 { 194 { 238 return << 195 G4double cross = Z*ComputeCrossSectionPerElectron 239 Z*ComputeCrossSectionPerElectron(p,kinetic << 196 (p,kineticEnergy,cutEnergy,maxEnergy); >> 197 return cross; 240 } 198 } 241 199 242 //....oooOO0OOooo........oooOO0OOooo........oo 200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 243 201 244 G4double G4BraggModel::CrossSectionPerVolume(c << 202 G4double G4BraggModel::CrossSectionPerVolume( 245 c << 203 const G4Material* material, 246 G << 204 const G4ParticleDefinition* p, 247 G << 205 G4double kineticEnergy, 248 G << 206 G4double cutEnergy, >> 207 G4double maxEnergy) 249 { 208 { 250 return material->GetElectronDensity() << 209 G4double eDensity = material->GetElectronDensity(); 251 *ComputeCrossSectionPerElectron(p,kineticE << 210 G4double cross = eDensity*ComputeCrossSectionPerElectron >> 211 (p,kineticEnergy,cutEnergy,maxEnergy); >> 212 return cross; 252 } 213 } 253 214 254 //....oooOO0OOooo........oooOO0OOooo........oo 215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 255 216 256 G4double G4BraggModel::ComputeDEDXPerVolume(co 217 G4double G4BraggModel::ComputeDEDXPerVolume(const G4Material* material, 257 co << 218 const G4ParticleDefinition* p, 258 G4 << 219 G4double kineticEnergy, 259 G4 << 220 G4double cutEnergy) 260 { << 221 { 261 const G4double tmax = MaxSecondaryEnergy(p, << 222 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy); 262 const G4double tkin = kinEnergy/massRate; << 223 G4double tkin = kineticEnergy/massRate; 263 const G4double cutEnergy = std::max(cut, low << 264 G4double dedx = 0.0; 224 G4double dedx = 0.0; 265 225 266 // tkin is the scaled energy to proton << 226 if(tkin < lowestKinEnergy) { 267 if (tkin < lowestKinEnergy) { << 227 dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy); 268 dedx = DEDX(material, lowestKinEnergy)*std << 269 } else { 228 } else { 270 dedx = DEDX(material, tkin); << 229 dedx = DEDX(material, tkin); >> 230 } 271 231 272 // correction for low cut value using Beth << 232 if (cutEnergy < tmax) { 273 if (cutEnergy < tmax) { << 274 const G4double tau = kinEnergy/mass; << 275 const G4double x = cutEnergy/tmax; << 276 233 277 dedx += (G4Log(x)*(tau + 1.)*(tau + 1.)/ << 234 G4double tau = kineticEnergy/mass; 278 CLHEP::twopi_mc2_rcl2 * material->GetElectro << 235 G4double gam = tau + 1.0; 279 } << 236 G4double bg2 = tau * (tau+2.0); >> 237 G4double beta2 = bg2/(gam*gam); >> 238 G4double x = cutEnergy/tmax; >> 239 >> 240 dedx += (G4Log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2 >> 241 * (material->GetElectronDensity())/beta2; 280 } 242 } 281 dedx = std::max(dedx, 0.0) * chargeSquare; << 243 >> 244 // now compute the total ionization loss >> 245 >> 246 if (dedx < 0.0) { dedx = 0.0; } >> 247 >> 248 dedx *= chargeSquare; 282 249 283 //G4cout << "E(MeV)= " << tkin/MeV << " dedx 250 //G4cout << "E(MeV)= " << tkin/MeV << " dedx= " << dedx 284 // << " " << material->GetName() << << 251 // << " " << material->GetName() << G4endl; >> 252 285 return dedx; 253 return dedx; 286 } 254 } 287 255 288 //....oooOO0OOooo........oooOO0OOooo........oo 256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 289 257 290 void G4BraggModel::SampleSecondaries(std::vect << 258 void G4BraggModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp, 291 const G4M << 259 const G4MaterialCutsCouple* couple, 292 const G4D << 260 const G4DynamicParticle* dp, 293 G4double << 261 G4double xmin, 294 G4double << 262 G4double maxEnergy) 295 { << 263 { 296 const G4double tmax = MaxSecondaryKinEnergy( << 264 G4double tmax = MaxSecondaryKinEnergy(dp); 297 const G4double xmax = std::min(tmax, maxEner << 265 G4double xmax = std::min(tmax, maxEnergy); 298 const G4double xmin = std::max(lowestKinEner << 299 if(xmin >= xmax) { return; } 266 if(xmin >= xmax) { return; } 300 267 301 G4double kineticEnergy = dp->GetKineticEnerg 268 G4double kineticEnergy = dp->GetKineticEnergy(); 302 const G4double energy = kineticEnergy + mas << 269 G4double energy = kineticEnergy + mass; 303 const G4double energy2 = energy*energy; << 270 G4double energy2 = energy*energy; 304 const G4double beta2 = kineticEnergy*(kine << 271 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2; 305 const G4double grej = 1.0; << 272 G4double grej = 1.0; 306 G4double deltaKinEnergy, f; 273 G4double deltaKinEnergy, f; 307 274 308 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra << 309 G4double rndm[2]; << 310 << 311 // sampling follows ... 275 // sampling follows ... 312 do { 276 do { 313 rndmEngineMod->flatArray(2, rndm); << 277 G4double q = G4UniformRand(); 314 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - rn << 278 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q); 315 279 316 f = 1.0 - beta2*deltaKinEnergy/tmax; 280 f = 1.0 - beta2*deltaKinEnergy/tmax; 317 281 318 if(f > grej) { 282 if(f > grej) { 319 G4cout << "G4BraggModel::SampleSeconda 283 G4cout << "G4BraggModel::SampleSecondary Warning! " 320 << "Majorant " << grej << " < " 284 << "Majorant " << grej << " < " 321 << f << " for e= " << deltaKinE 285 << f << " for e= " << deltaKinEnergy 322 << G4endl; 286 << G4endl; 323 } 287 } 324 288 325 // Loop checking, 03-Aug-2015, Vladimir Iv << 289 } while( grej*G4UniformRand() >= f ); 326 } while( grej*rndm[1] >= f ); << 327 290 328 G4ThreeVector deltaDirection; 291 G4ThreeVector deltaDirection; 329 292 330 if(UseAngularGeneratorFlag()) { 293 if(UseAngularGeneratorFlag()) { 331 const G4Material* mat = couple->GetMateri 294 const G4Material* mat = couple->GetMaterial(); 332 G4int Z = SelectRandomAtomNumber(mat); 295 G4int Z = SelectRandomAtomNumber(mat); 333 296 334 deltaDirection = 297 deltaDirection = 335 GetAngularDistribution()->SampleDirectio 298 GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat); 336 299 337 } else { 300 } else { 338 301 339 G4double deltaMomentum = 302 G4double deltaMomentum = 340 std::sqrt(deltaKinEnergy * (deltaKinEner << 303 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2)); 341 G4double cost = deltaKinEnergy * (energy + 304 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) / 342 (deltaMomentum * dp->GetTotalMomentum()) 305 (deltaMomentum * dp->GetTotalMomentum()); 343 if(cost > 1.0) { cost = 1.0; } 306 if(cost > 1.0) { cost = 1.0; } 344 G4double sint = std::sqrt((1.0 - cost)*(1. << 307 G4double sint = sqrt((1.0 - cost)*(1.0 + cost)); 345 308 346 G4double phi = twopi*rndmEngineMod->flat() << 309 G4double phi = twopi * G4UniformRand() ; 347 310 348 deltaDirection.set(sint*std::cos(phi),sint << 311 deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ; 349 deltaDirection.rotateUz(dp->GetMomentumDir 312 deltaDirection.rotateUz(dp->GetMomentumDirection()); 350 } 313 } 351 314 352 // create G4DynamicParticle object for delta 315 // create G4DynamicParticle object for delta ray 353 auto delta = new G4DynamicParticle(theElectr << 316 G4DynamicParticle* delta = >> 317 new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy); 354 318 355 // Change kinematics of primary particle 319 // Change kinematics of primary particle 356 kineticEnergy -= deltaKinEnergy; 320 kineticEnergy -= deltaKinEnergy; 357 G4ThreeVector finalP = dp->GetMomentum() - d 321 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum(); 358 finalP = finalP.unit(); << 322 finalP = finalP.unit(); 359 323 360 fParticleChange->SetProposedKineticEnergy(ki 324 fParticleChange->SetProposedKineticEnergy(kineticEnergy); 361 fParticleChange->SetProposedMomentumDirectio 325 fParticleChange->SetProposedMomentumDirection(finalP); 362 326 363 vdp->push_back(delta); 327 vdp->push_back(delta); 364 } 328 } 365 329 366 //....oooOO0OOooo........oooOO0OOooo........oo 330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 367 331 368 G4double G4BraggModel::MaxSecondaryEnergy(cons 332 G4double G4BraggModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd, 369 G4do << 333 G4double kinEnergy) 370 { 334 { 371 if(pd != particle) { SetParticle(pd); } 335 if(pd != particle) { SetParticle(pd); } 372 G4double tau = kinEnergy/mass; 336 G4double tau = kinEnergy/mass; 373 G4double tmax = 2.0*electron_mass_c2*tau*(ta 337 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) / 374 (1. + 2.0*(tau + 1.)*ratio + 338 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio); 375 return tmax; 339 return tmax; 376 } 340 } 377 341 378 //....oooOO0OOooo........oooOO0OOooo........oo 342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 379 343 380 void G4BraggModel::HasMaterial(const G4Materia << 344 G4bool G4BraggModel::HasMaterial(const G4Material* material) 381 { 345 { 382 const G4String& chFormula = mat->GetChemical << 346 G4String chFormula = material->GetChemicalFormula(); 383 if(chFormula.empty()) { return; } << 347 if("" == chFormula) { return false; } 384 348 385 // ICRU Report N49, 1993. Power's model for << 349 // ICRU Report N49, 1993. Power's model for He. 386 static const G4int numberOfMolecula = 11; << 350 static const size_t numberOfMolecula = 11; 387 static const G4String molName[numberOfMolecu 351 static const G4String molName[numberOfMolecula] = { 388 "Al_2O_3", "CO_2", 352 "Al_2O_3", "CO_2", "CH_4", 389 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Pol 353 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polypropylene", "(C_8H_8)_N", 390 "C_3H_8", "SiO_2", 354 "C_3H_8", "SiO_2", "H_2O", 391 "H_2O-Gas", "Graphite" } ; 355 "H_2O-Gas", "Graphite" } ; >> 356 static const G4int idxPSTAR[numberOfMolecula] = { >> 357 6, 16, 36, >> 358 52, 55, 56, >> 359 59, 62, 71, >> 360 72, 13}; >> 361 >> 362 // Special treatment for water in gas state >> 363 const G4State theState = material->GetState() ; >> 364 >> 365 if( theState == kStateGas && "H_2O" == chFormula) { >> 366 chFormula = G4String("H_2O-Gas"); >> 367 } 392 368 393 // Search for the material in the table 369 // Search for the material in the table 394 for (G4int i=0; i<numberOfMolecula; ++i) { << 370 for (size_t i=0; i<numberOfMolecula; ++i) { 395 if (chFormula == molName[i]) { 371 if (chFormula == molName[i]) { 396 iMolecula = i; << 372 iMolecula = -1; 397 return; << 373 iPSTAR = idxPSTAR[i]; >> 374 return true; 398 } 375 } 399 } 376 } 400 return; << 377 return false; 401 } 378 } 402 379 403 //....oooOO0OOooo........oooOO0OOooo........oo 380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 404 381 405 G4double G4BraggModel::StoppingPower(const G4M 382 G4double G4BraggModel::StoppingPower(const G4Material* material, 406 G4double << 383 G4double kineticEnergy) 407 { 384 { 408 G4double ionloss = 0.0 ; 385 G4double ionloss = 0.0 ; 409 386 410 if (iMolecula >= 0) { 387 if (iMolecula >= 0) { 411 388 412 // The data and the fit from: 389 // The data and the fit from: 413 // ICRU Report N49, 1993. Ziegler's model 390 // ICRU Report N49, 1993. Ziegler's model for protons. 414 // Proton kinetic energy for parametrisati 391 // Proton kinetic energy for parametrisation (keV/amu) 415 392 416 G4double T = kineticEnergy/(keV*protonMass 393 G4double T = kineticEnergy/(keV*protonMassAMU) ; 417 394 418 static const G4float a[11][5] = { << 395 static const G4double a[11][5] = { 419 {1.187E+1f, 1.343E+1f, 1.069E+4f, 7.723E+2f << 396 {1.187E+1, 1.343E+1, 1.069E+4, 7.723E+2, 2.153E-2}, 420 {7.802E+0f, 8.814E+0f, 8.303E+3f, 7.446E+2f << 397 {7.802E+0, 8.814E+0, 8.303E+3, 7.446E+2, 7.966E-3}, 421 {7.294E+0f, 8.284E+0f, 5.010E+3f, 4.544E+2f << 398 {7.294E+0, 8.284E+0, 5.010E+3, 4.544E+2, 8.153E-3}, 422 {8.646E+0f, 9.800E+0f, 7.066E+3f, 4.581E+2f << 399 {8.646E+0, 9.800E+0, 7.066E+3, 4.581E+2, 9.383E-3}, 423 {1.286E+1f, 1.462E+1f, 5.625E+3f, 2.621E+3f << 400 {1.286E+1, 1.462E+1, 5.625E+3, 2.621E+3, 3.512E-2}, 424 {3.229E+1f, 3.696E+1f, 8.918E+3f, 3.244E+3f << 401 {3.229E+1, 3.696E+1, 8.918E+3, 3.244E+3, 1.273E-1}, 425 {1.604E+1f, 1.825E+1f, 6.967E+3f, 2.307E+3f << 402 {1.604E+1, 1.825E+1, 6.967E+3, 2.307E+3, 3.775E-2}, 426 {8.049E+0f, 9.099E+0f, 9.257E+3f, 3.846E+2f << 403 {8.049E+0, 9.099E+0, 9.257E+3, 3.846E+2, 1.007E-2}, 427 {4.015E+0f, 4.542E+0f, 3.955E+3f, 4.847E+2f << 404 {4.015E+0, 4.542E+0, 3.955E+3, 4.847E+2, 7.904E-3}, 428 {4.571E+0f, 5.173E+0f, 4.346E+3f, 4.779E+2f << 405 {4.571E+0, 5.173E+0, 4.346E+3, 4.779E+2, 8.572E-3}, 429 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f << 406 {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2} }; 430 << 407 431 static const G4float atomicWeight[11] = { << 408 static const G4double atomicWeight[11] = { 432 101.96128f, 44.0098f, 16.0426f, 28.0536f, << 409 101.96128, 44.0098, 16.0426, 28.0536, 42.0804, 433 104.1512f, 44.665f, 60.0843f, 18.0152f, 18 << 410 104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0}; 434 411 435 if ( T < 10.0 ) { 412 if ( T < 10.0 ) { 436 ionloss = ((G4double)(a[iMolecula][0])) << 413 ionloss = a[iMolecula][0] * sqrt(T) ; 437 414 438 } else if ( T < 10000.0 ) { 415 } else if ( T < 10000.0 ) { 439 G4double x1 = (G4double)(a[iMolecula][1] << 416 G4double slow = a[iMolecula][1] * G4Exp(G4Log(T)* 0.45); 440 G4double x2 = (G4double)(a[iMolecula][2] << 417 G4double shigh = G4Log( 1.0 + a[iMolecula][3]/T 441 G4double x3 = (G4double)(a[iMolecula][3] << 418 + a[iMolecula][4]*T ) * a[iMolecula][2]/T ; 442 G4double x4 = (G4double)(a[iMolecula][4] << 443 G4double slow = x1 * G4Exp(G4Log(T)* 0. << 444 G4double shigh = G4Log( 1.0 + x3/T + x4 << 445 ionloss = slow*shigh / (slow + shigh) ; 419 ionloss = slow*shigh / (slow + shigh) ; 446 } 420 } 447 421 448 ionloss = std::max(ionloss, 0.0); << 422 if ( ionloss < 0.0) ionloss = 0.0 ; 449 if ( 10 == iMolecula ) { 423 if ( 10 == iMolecula ) { 450 static const G4double invLog10 = 1.0/G4L << 451 << 452 if (T < 100.0) { 424 if (T < 100.0) { 453 ionloss *= (1.0+0.023+0.0066*G4Log(T)* << 425 ionloss *= (1.0+0.023+0.0066*G4Log(T)*invLog10); 454 } 426 } 455 else if (T < 700.0) { 427 else if (T < 700.0) { 456 ionloss *=(1.0+0.089-0.0248*G4Log(T-99 << 428 ionloss *=(1.0+0.089-0.0248*G4Log(T-99.)*invLog10); 457 } 429 } 458 else if (T < 10000.0) { 430 else if (T < 10000.0) { 459 ionloss *=(1.0+0.089-0.0248*G4Log(700. << 431 ionloss *=(1.0+0.089-0.0248*G4Log(700.-99.)*invLog10); 460 } 432 } 461 } 433 } 462 ionloss /= (G4double)atomicWeight[iMolecul << 434 ionloss /= atomicWeight[iMolecula]; 463 435 464 // pure material (normally not the case for 436 // pure material (normally not the case for this function) 465 } else if(1 == (material->GetNumberOfElement 437 } else if(1 == (material->GetNumberOfElements())) { 466 G4double z = material->GetZ() ; 438 G4double z = material->GetZ() ; 467 ionloss = ElectronicStoppingPower( z, kine 439 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ; 468 } 440 } 469 441 470 return ionloss; 442 return ionloss; 471 } 443 } 472 444 473 //....oooOO0OOooo........oooOO0OOooo........oo 445 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 474 446 475 G4double G4BraggModel::ElectronicStoppingPower 447 G4double G4BraggModel::ElectronicStoppingPower(G4double z, 476 448 G4double kineticEnergy) const 477 { 449 { 478 G4double ionloss ; 450 G4double ionloss ; 479 G4int i = std::min(std::max(G4lrint(z)-1,0), << 451 G4int i = G4int(z)-1 ; // index of atom 480 << 452 if(i < 0) i = 0 ; >> 453 if(i > 91) i = 91 ; >> 454 481 // The data and the fit from: 455 // The data and the fit from: 482 // ICRU Report 49, 1993. Ziegler's type of p 456 // ICRU Report 49, 1993. Ziegler's type of parametrisations. 483 // Proton kinetic energy for parametrisation 457 // Proton kinetic energy for parametrisation (keV/amu) 484 458 485 G4double T = kineticEnergy/(keV*protonMassAM 459 G4double T = kineticEnergy/(keV*protonMassAMU) ; 486 460 487 static const G4float a[92][5] = { << 461 static const G4double a[92][5] = { 488 {1.254E+0f, 1.440E+0f, 2.426E+2f, 1.200E+4f << 462 {1.254E+0, 1.440E+0, 2.426E+2, 1.200E+4, 1.159E-1}, 489 {1.229E+0f, 1.397E+0f, 4.845E+2f, 5.873E+3f << 463 {1.229E+0, 1.397E+0, 4.845E+2, 5.873E+3, 5.225E-2}, 490 {1.411E+0f, 1.600E+0f, 7.256E+2f, 3.013E+3f << 464 {1.411E+0, 1.600E+0, 7.256E+2, 3.013E+3, 4.578E-2}, 491 {2.248E+0f, 2.590E+0f, 9.660E+2f, 1.538E+2f << 465 {2.248E+0, 2.590E+0, 9.660E+2, 1.538E+2, 3.475E-2}, 492 {2.474E+0f, 2.815E+0f, 1.206E+3f, 1.060E+3f << 466 {2.474E+0, 2.815E+0, 1.206E+3, 1.060E+3, 2.855E-2}, 493 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f << 467 {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2}, 494 {2.954E+0f, 3.350E+0f, 1.683E+3f, 1.900E+3f << 468 {2.954E+0, 3.350E+0, 1.683E+3, 1.900E+3, 2.513E-2}, 495 {2.652E+0f, 3.000E+0f, 1.920E+3f, 2.000E+3f << 469 {2.652E+0, 3.000E+0, 1.920E+3, 2.000E+3, 2.230E-2}, 496 {2.085E+0f, 2.352E+0f, 2.157E+3f, 2.634E+3f << 470 {2.085E+0, 2.352E+0, 2.157E+3, 2.634E+3, 1.816E-2}, 497 {1.951E+0f, 2.199E+0f, 2.393E+3f, 2.699E+3f << 471 {1.951E+0, 2.199E+0, 2.393E+3, 2.699E+3, 1.568E-2}, 498 // Z= 11-20 472 // Z= 11-20 499 {2.542E+0f, 2.869E+0f, 2.628E+3f, 1.854E+3f << 473 {2.542E+0, 2.869E+0, 2.628E+3, 1.854E+3, 1.472E-2}, 500 {3.791E+0f, 4.293E+0f, 2.862E+3f, 1.009E+3f << 474 {3.791E+0, 4.293E+0, 2.862E+3, 1.009E+3, 1.397E-2}, 501 {4.154E+0f, 4.739E+0f, 2.766E+3f, 1.645E+2f << 475 {4.154E+0, 4.739E+0, 2.766E+3, 1.645E+2, 2.023E-2}, 502 {4.914E+0f, 5.598E+0f, 3.193E+3f, 2.327E+2f << 476 {4.914E+0, 5.598E+0, 3.193E+3, 2.327E+2, 1.419E-2}, 503 {3.232E+0f, 3.647E+0f, 3.561E+3f, 1.560E+3f << 477 {3.232E+0, 3.647E+0, 3.561E+3, 1.560E+3, 1.267E-2}, 504 {3.447E+0f, 3.891E+0f, 3.792E+3f, 1.219E+3f << 478 {3.447E+0, 3.891E+0, 3.792E+3, 1.219E+3, 1.211E-2}, 505 {5.301E+0f, 6.008E+0f, 3.969E+3f, 6.451E+2f << 479 {5.301E+0, 6.008E+0, 3.969E+3, 6.451E+2, 1.183E-2}, 506 {5.731E+0f, 6.500E+0f, 4.253E+3f, 5.300E+2f << 480 {5.731E+0, 6.500E+0, 4.253E+3, 5.300E+2, 1.123E-2}, 507 {5.152E+0f, 5.833E+0f, 4.482E+3f, 5.457E+2f << 481 {5.152E+0, 5.833E+0, 4.482E+3, 5.457E+2, 1.129E-2}, 508 {5.521E+0f, 6.252E+0f, 4.710E+3f, 5.533E+2f << 482 {5.521E+0, 6.252E+0, 4.710E+3, 5.533E+2, 1.112E-2}, 509 // Z= 21-30 483 // Z= 21-30 510 {5.201E+0f, 5.884E+0f, 4.938E+3f, 5.609E+2f << 484 {5.201E+0, 5.884E+0, 4.938E+3, 5.609E+2, 9.995E-3}, 511 {4.858E+0f, 5.489E+0f, 5.260E+3f, 6.511E+2f << 485 {4.858E+0, 5.489E+0, 5.260E+3, 6.511E+2, 8.930E-3}, 512 {4.479E+0f, 5.055E+0f, 5.391E+3f, 9.523E+2f << 486 {4.479E+0, 5.055E+0, 5.391E+3, 9.523E+2, 9.117E-3}, 513 {3.983E+0f, 4.489E+0f, 5.616E+3f, 1.336E+3f << 487 {3.983E+0, 4.489E+0, 5.616E+3, 1.336E+3, 8.413E-3}, 514 {3.469E+0f, 3.907E+0f, 5.725E+3f, 1.461E+3f << 488 {3.469E+0, 3.907E+0, 5.725E+3, 1.461E+3, 8.829E-3}, 515 {3.519E+0f, 3.963E+0f, 6.065E+3f, 1.243E+3f << 489 {3.519E+0, 3.963E+0, 6.065E+3, 1.243E+3, 7.782E-3}, 516 {3.140E+0f, 3.535E+0f, 6.288E+3f, 1.372E+3f << 490 {3.140E+0, 3.535E+0, 6.288E+3, 1.372E+3, 7.361E-3}, 517 {3.553E+0f, 4.004E+0f, 6.205E+3f, 5.551E+2f << 491 {3.553E+0, 4.004E+0, 6.205E+3, 5.551E+2, 8.763E-3}, 518 {3.696E+0f, 4.194E+0f, 4.649E+3f, 8.113E+1f << 492 {3.696E+0, 4.194E+0, 4.649E+3, 8.113E+1, 2.242E-2}, 519 {4.210E+0f, 4.750E+0f, 6.953E+3f, 2.952E+2f << 493 {4.210E+0, 4.750E+0, 6.953E+3, 2.952E+2, 6.809E-3}, 520 // Z= 31-40 494 // Z= 31-40 521 {5.041E+0f, 5.697E+0f, 7.173E+3f, 2.026E+2f << 495 {5.041E+0, 5.697E+0, 7.173E+3, 2.026E+2, 6.725E-3}, 522 {5.554E+0f, 6.300E+0f, 6.496E+3f, 1.100E+2f << 496 {5.554E+0, 6.300E+0, 6.496E+3, 1.100E+2, 9.689E-3}, 523 {5.323E+0f, 6.012E+0f, 7.611E+3f, 2.925E+2f << 497 {5.323E+0, 6.012E+0, 7.611E+3, 2.925E+2, 6.447E-3}, 524 {5.874E+0f, 6.656E+0f, 7.395E+3f, 1.175E+2f << 498 {5.874E+0, 6.656E+0, 7.395E+3, 1.175E+2, 7.684E-3}, 525 {6.658E+0f, 7.536E+0f, 7.694E+3f, 2.223E+2f << 499 {6.658E+0, 7.536E+0, 7.694E+3, 2.223E+2, 6.509E-3}, 526 {6.413E+0f, 7.240E+0f, 1.185E+4f, 1.537E+2f << 500 {6.413E+0, 7.240E+0, 1.185E+4, 1.537E+2, 2.880E-3}, 527 {5.694E+0f, 6.429E+0f, 8.478E+3f, 2.929E+2f << 501 {5.694E+0, 6.429E+0, 8.478E+3, 2.929E+2, 6.087E-3}, 528 {6.339E+0f, 7.159E+0f, 8.693E+3f, 3.303E+2f << 502 {6.339E+0, 7.159E+0, 8.693E+3, 3.303E+2, 6.003E-3}, 529 {6.407E+0f, 7.234E+0f, 8.907E+3f, 3.678E+2f << 503 {6.407E+0, 7.234E+0, 8.907E+3, 3.678E+2, 5.889E-3}, 530 {6.734E+0f, 7.603E+0f, 9.120E+3f, 4.052E+2f << 504 {6.734E+0, 7.603E+0, 9.120E+3, 4.052E+2, 5.765E-3}, 531 // Z= 41-50 505 // Z= 41-50 532 {6.901E+0f, 7.791E+0f, 9.333E+3f, 4.427E+2f << 506 {6.901E+0, 7.791E+0, 9.333E+3, 4.427E+2, 5.587E-3}, 533 {6.424E+0f, 7.248E+0f, 9.545E+3f, 4.802E+2f << 507 {6.424E+0, 7.248E+0, 9.545E+3, 4.802E+2, 5.376E-3}, 534 {6.799E+0f, 7.671E+0f, 9.756E+3f, 5.176E+2f << 508 {6.799E+0, 7.671E+0, 9.756E+3, 5.176E+2, 5.315E-3}, 535 {6.109E+0f, 6.887E+0f, 9.966E+3f, 5.551E+2f << 509 {6.109E+0, 6.887E+0, 9.966E+3, 5.551E+2, 5.151E-3}, 536 {5.924E+0f, 6.677E+0f, 1.018E+4f, 5.925E+2f << 510 {5.924E+0, 6.677E+0, 1.018E+4, 5.925E+2, 4.919E-3}, 537 {5.238E+0f, 5.900E+0f, 1.038E+4f, 6.300E+2f << 511 {5.238E+0, 5.900E+0, 1.038E+4, 6.300E+2, 4.758E-3}, 538 // {5.623f, 6.354f, 7160.0f, 337.6f << 512 // {5.623, 6.354, 7160.0, 337.6, 0.013940}, // Ag Ziegler77 539 {5.345E+0f, 6.038E+0f, 6.790E+3f, 3.978E+2f << 513 {5.345E+0, 6.038E+0, 6.790E+3, 3.978E+2, 1.676E-2}, // Ag ICRU49 540 {5.814E+0f, 6.554E+0f, 1.080E+4f, 3.555E+2f << 514 {5.814E+0, 6.554E+0, 1.080E+4, 3.555E+2, 4.626E-3}, 541 {6.229E+0f, 7.024E+0f, 1.101E+4f, 3.709E+2f << 515 {6.229E+0, 7.024E+0, 1.101E+4, 3.709E+2, 4.540E-3}, 542 {6.409E+0f, 7.227E+0f, 1.121E+4f, 3.864E+2f << 516 {6.409E+0, 7.227E+0, 1.121E+4, 3.864E+2, 4.474E-3}, 543 // Z= 51-60 517 // Z= 51-60 544 {7.500E+0f, 8.480E+0f, 8.608E+3f, 3.480E+2f << 518 {7.500E+0, 8.480E+0, 8.608E+3, 3.480E+2, 9.074E-3}, 545 {6.979E+0f, 7.871E+0f, 1.162E+4f, 3.924E+2f << 519 {6.979E+0, 7.871E+0, 1.162E+4, 3.924E+2, 4.402E-3}, 546 {7.725E+0f, 8.716E+0f, 1.183E+4f, 3.948E+2f << 520 {7.725E+0, 8.716E+0, 1.183E+4, 3.948E+2, 4.376E-3}, 547 {8.337E+0f, 9.425E+0f, 1.051E+4f, 2.696E+2f << 521 {8.337E+0, 9.425E+0, 1.051E+4, 2.696E+2, 6.206E-3}, 548 {7.287E+0f, 8.218E+0f, 1.223E+4f, 3.997E+2f << 522 {7.287E+0, 8.218E+0, 1.223E+4, 3.997E+2, 4.447E-3}, 549 {7.899E+0f, 8.911E+0f, 1.243E+4f, 4.021E+2f << 523 {7.899E+0, 8.911E+0, 1.243E+4, 4.021E+2, 4.511E-3}, 550 {8.041E+0f, 9.071E+0f, 1.263E+4f, 4.045E+2f << 524 {8.041E+0, 9.071E+0, 1.263E+4, 4.045E+2, 4.540E-3}, 551 {7.488E+0f, 8.444E+0f, 1.283E+4f, 4.069E+2f << 525 {7.488E+0, 8.444E+0, 1.283E+4, 4.069E+2, 4.420E-3}, 552 {7.291E+0f, 8.219E+0f, 1.303E+4f, 4.093E+2f << 526 {7.291E+0, 8.219E+0, 1.303E+4, 4.093E+2, 4.298E-3}, 553 {7.098E+0f, 8.000E+0f, 1.323E+4f, 4.118E+2f << 527 {7.098E+0, 8.000E+0, 1.323E+4, 4.118E+2, 4.182E-3}, 554 // Z= 61-70 528 // Z= 61-70 555 {6.909E+0f, 7.786E+0f, 1.343E+4f, 4.142E+2f << 529 {6.909E+0, 7.786E+0, 1.343E+4, 4.142E+2, 4.058E-3}, 556 {6.728E+0f, 7.580E+0f, 1.362E+4f, 4.166E+2f << 530 {6.728E+0, 7.580E+0, 1.362E+4, 4.166E+2, 3.976E-3}, 557 {6.551E+0f, 7.380E+0f, 1.382E+4f, 4.190E+2f << 531 {6.551E+0, 7.380E+0, 1.382E+4, 4.190E+2, 3.877E-3}, 558 {6.739E+0f, 7.592E+0f, 1.402E+4f, 4.214E+2f << 532 {6.739E+0, 7.592E+0, 1.402E+4, 4.214E+2, 3.863E-3}, 559 {6.212E+0f, 6.996E+0f, 1.421E+4f, 4.239E+2f << 533 {6.212E+0, 6.996E+0, 1.421E+4, 4.239E+2, 3.725E-3}, 560 {5.517E+0f, 6.210E+0f, 1.440E+4f, 4.263E+2f << 534 {5.517E+0, 6.210E+0, 1.440E+4, 4.263E+2, 3.632E-3}, 561 {5.220E+0f, 5.874E+0f, 1.460E+4f, 4.287E+2f << 535 {5.220E+0, 5.874E+0, 1.460E+4, 4.287E+2, 3.498E-3}, 562 {5.071E+0f, 5.706E+0f, 1.479E+4f, 4.330E+2f << 536 {5.071E+0, 5.706E+0, 1.479E+4, 4.330E+2, 3.405E-3}, 563 {4.926E+0f, 5.542E+0f, 1.498E+4f, 4.335E+2f << 537 {4.926E+0, 5.542E+0, 1.498E+4, 4.335E+2, 3.342E-3}, 564 {4.788E+0f, 5.386E+0f, 1.517E+4f, 4.359E+2f << 538 {4.788E+0, 5.386E+0, 1.517E+4, 4.359E+2, 3.292E-3}, 565 // Z= 71-80 539 // Z= 71-80 566 {4.893E+0f, 5.505E+0f, 1.536E+4f, 4.384E+2f << 540 {4.893E+0, 5.505E+0, 1.536E+4, 4.384E+2, 3.243E-3}, 567 {5.028E+0f, 5.657E+0f, 1.555E+4f, 4.408E+2f << 541 {5.028E+0, 5.657E+0, 1.555E+4, 4.408E+2, 3.195E-3}, 568 {4.738E+0f, 5.329E+0f, 1.574E+4f, 4.432E+2f << 542 {4.738E+0, 5.329E+0, 1.574E+4, 4.432E+2, 3.186E-3}, 569 {4.587E+0f, 5.160E+0f, 1.541E+4f, 4.153E+2f << 543 {4.587E+0, 5.160E+0, 1.541E+4, 4.153E+2, 3.406E-3}, 570 {5.201E+0f, 5.851E+0f, 1.612E+4f, 4.416E+2f << 544 {5.201E+0, 5.851E+0, 1.612E+4, 4.416E+2, 3.122E-3}, 571 {5.071E+0f, 5.704E+0f, 1.630E+4f, 4.409E+2f << 545 {5.071E+0, 5.704E+0, 1.630E+4, 4.409E+2, 3.082E-3}, 572 {4.946E+0f, 5.563E+0f, 1.649E+4f, 4.401E+2f << 546 {4.946E+0, 5.563E+0, 1.649E+4, 4.401E+2, 2.965E-3}, 573 {4.477E+0f, 5.034E+0f, 1.667E+4f, 4.393E+2f << 547 {4.477E+0, 5.034E+0, 1.667E+4, 4.393E+2, 2.871E-3}, 574 // {4.856f, 5.460f, 18320.0f, 438.5 << 548 // {4.856, 5.460, 18320.0, 438.5, 0.002542}, //Ziegler77 575 {4.844E+0f, 5.458E+0f, 7.852E+3f, 9.758E+2f << 549 {4.844E+0, 5.458E+0, 7.852E+3, 9.758E+2, 2.077E-2}, //ICRU49 576 {4.307E+0f, 4.843E+0f, 1.704E+4f, 4.878E+2f << 550 {4.307E+0, 4.843E+0, 1.704E+4, 4.878E+2, 2.882E-3}, 577 // Z= 81-90 551 // Z= 81-90 578 {4.723E+0f, 5.311E+0f, 1.722E+4f, 5.370E+2f << 552 {4.723E+0, 5.311E+0, 1.722E+4, 5.370E+2, 2.913E-3}, 579 {5.319E+0f, 5.982E+0f, 1.740E+4f, 5.863E+2f << 553 {5.319E+0, 5.982E+0, 1.740E+4, 5.863E+2, 2.871E-3}, 580 {5.956E+0f, 6.700E+0f, 1.780E+4f, 6.770E+2f << 554 {5.956E+0, 6.700E+0, 1.780E+4, 6.770E+2, 2.660E-3}, 581 {6.158E+0f, 6.928E+0f, 1.777E+4f, 5.863E+2f << 555 {6.158E+0, 6.928E+0, 1.777E+4, 5.863E+2, 2.812E-3}, 582 {6.203E+0f, 6.979E+0f, 1.795E+4f, 5.863E+2f << 556 {6.203E+0, 6.979E+0, 1.795E+4, 5.863E+2, 2.776E-3}, 583 {6.181E+0f, 6.954E+0f, 1.812E+4f, 5.863E+2f << 557 {6.181E+0, 6.954E+0, 1.812E+4, 5.863E+2, 2.748E-3}, 584 {6.949E+0f, 7.820E+0f, 1.830E+4f, 5.863E+2f << 558 {6.949E+0, 7.820E+0, 1.830E+4, 5.863E+2, 2.737E-3}, 585 {7.506E+0f, 8.448E+0f, 1.848E+4f, 5.863E+2f << 559 {7.506E+0, 8.448E+0, 1.848E+4, 5.863E+2, 2.727E-3}, 586 {7.648E+0f, 8.609E+0f, 1.866E+4f, 5.863E+2f << 560 {7.648E+0, 8.609E+0, 1.866E+4, 5.863E+2, 2.697E-3}, 587 {7.711E+0f, 8.679E+0f, 1.883E+4f, 5.863E+2f << 561 {7.711E+0, 8.679E+0, 1.883E+4, 5.863E+2, 2.641E-3}, 588 // Z= 91-92 562 // Z= 91-92 589 {7.407E+0f, 8.336E+0f, 1.901E+4f, 5.863E+2f << 563 {7.407E+0, 8.336E+0, 1.901E+4, 5.863E+2, 2.603E-3}, 590 {7.290E+0f, 8.204E+0f, 1.918E+4f, 5.863E+2f << 564 {7.290E+0, 8.204E+0, 1.918E+4, 5.863E+2, 2.673E-3} 591 }; 565 }; 592 566 593 G4double fac = 1.0 ; 567 G4double fac = 1.0 ; 594 568 595 // Carbon specific case for E < 40 keV 569 // Carbon specific case for E < 40 keV 596 if ( T < 40.0 && 5 == i) { 570 if ( T < 40.0 && 5 == i) { 597 fac = std::sqrt(T*0.025); << 571 fac = sqrt(T/40.0) ; 598 T = 40.0; << 572 T = 40.0 ; 599 573 600 // Free electron gas model 574 // Free electron gas model 601 } else if ( T < 10.0 ) { 575 } else if ( T < 10.0 ) { 602 fac = std::sqrt(T*0.1) ; << 576 fac = sqrt(T*0.1) ; 603 T = 10.0; << 577 T =10.0 ; 604 } 578 } 605 579 606 // Main parametrisation 580 // Main parametrisation 607 G4double x1 = (G4double)(a[i][1]); << 581 G4double slow = a[i][1] * G4Exp(G4Log(T) * 0.45) ; 608 G4double x2 = (G4double)(a[i][2]); << 582 G4double shigh = G4Log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ; 609 G4double x3 = (G4double)(a[i][3]); << 583 ionloss = slow*shigh*fac / (slow + shigh) ; 610 G4double x4 = (G4double)(a[i][4]); << 611 G4double slow = x1 * G4Exp(G4Log(T) * 0.45) << 612 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) << 613 ionloss = slow*shigh*fac / (slow + shigh); << 614 584 615 ionloss = std::max(ionloss, 0.0); << 585 if ( ionloss < 0.0) { ionloss = 0.0; } 616 586 617 return ionloss; 587 return ionloss; 618 } 588 } 619 589 620 //....oooOO0OOooo........oooOO0OOooo........oo 590 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 621 591 622 G4double G4BraggModel::DEDX(const G4Material* 592 G4double G4BraggModel::DEDX(const G4Material* material, G4double kineticEnergy) 623 { 593 { 624 G4double eloss = 0.0; 594 G4double eloss = 0.0; 625 595 626 // check DB 596 // check DB 627 if(material != currentMaterial) { 597 if(material != currentMaterial) { 628 currentMaterial = material; 598 currentMaterial = material; 629 baseMaterial = material->GetBaseMaterial() << 630 ? material->GetBaseMaterial() : material << 631 iPSTAR = -1; 599 iPSTAR = -1; 632 iMolecula = -1; 600 iMolecula = -1; 633 iICRU90 = fICRU90 ? fICRU90->GetIndex(base << 601 if( !HasMaterial(material) ) { iPSTAR = pstar.GetIndex(material); } 634 << 602 635 if(iICRU90 < 0) { << 636 iPSTAR = fPSTAR->GetIndex(baseMaterial); << 637 if(iPSTAR < 0) { HasMaterial(baseMateria << 638 } << 639 //G4cout << "%%% " <<material->GetName() < 603 //G4cout << "%%% " <<material->GetName() << " iMolecula= " 640 // << iMolecula << " iPSTAR= " << i << 604 // << iMolecula << " iPSTAR= " << iPSTAR << G4endl; 641 // << " iICRU90= " << iICRU90<< G4e << 642 } << 643 605 644 // ICRU90 parameterisation << 645 if(iICRU90 >= 0) { << 646 return fICRU90->GetElectronicDEDXforProton << 647 *material->GetDensity(); << 648 } 606 } 649 // PSTAR parameterisation << 650 if( iPSTAR >= 0 ) { << 651 return fPSTAR->GetElectronicDEDX(iPSTAR, k << 652 *material->GetDensity(); << 653 607 654 } << 608 const G4int numberOfElements = material->GetNumberOfElements(); 655 const std::size_t numberOfElements = materia << 656 const G4double* theAtomicNumDensityVector = 609 const G4double* theAtomicNumDensityVector = 657 material->Get 610 material->GetAtomicNumDensityVector(); 658 611 >> 612 if( iPSTAR >= 0 ) { >> 613 return pstar.GetElectronicDEDX(iPSTAR, kineticEnergy)*material->GetDensity(); >> 614 >> 615 } else if(iMolecula >= 0) { 659 616 660 if(iMolecula >= 0) { << 617 eloss = StoppingPower(material, kineticEnergy)* 661 eloss = StoppingPower(baseMaterial, kineti << 662 material->GetDensity 618 material->GetDensity()/amu; 663 619 664 // Pure material ICRU49 paralmeterisation << 620 // Pure material ICRU49 paralmeterisation 665 } else if(1 == numberOfElements) { 621 } else if(1 == numberOfElements) { 666 622 667 G4double z = material->GetZ(); 623 G4double z = material->GetZ(); 668 eloss = ElectronicStoppingPower(z, kinetic 624 eloss = ElectronicStoppingPower(z, kineticEnergy) 669 * (material->Ge 625 * (material->GetTotNbOfAtomsPerVolume()); 670 626 671 627 672 // Experimental data exist only for kinetic 628 // Experimental data exist only for kinetic energy 125 keV 673 } else if( MolecIsInZiegler1988(material) ) 629 } else if( MolecIsInZiegler1988(material) ) { 674 630 675 // Loop over elements - calculation based 631 // Loop over elements - calculation based on Bragg's rule 676 G4double eloss125 = 0.0 ; 632 G4double eloss125 = 0.0 ; 677 const G4ElementVector* theElementVector = 633 const G4ElementVector* theElementVector = 678 material->GetElemen 634 material->GetElementVector(); 679 635 680 // Loop for the elements in the material 636 // Loop for the elements in the material 681 for (std::size_t i=0; i<numberOfElements; << 637 for (G4int i=0; i<numberOfElements; i++) { 682 const G4Element* element = (*theElementV 638 const G4Element* element = (*theElementVector)[i] ; 683 G4double z = element->GetZ() ; 639 G4double z = element->GetZ() ; 684 eloss += ElectronicStoppingPower(z,ki 640 eloss += ElectronicStoppingPower(z,kineticEnergy) 685 * theAtomi 641 * theAtomicNumDensityVector[i] ; 686 eloss125 += ElectronicStoppingPower(z,12 642 eloss125 += ElectronicStoppingPower(z,125.0*keV) 687 * theAtomi 643 * theAtomicNumDensityVector[i] ; 688 } 644 } 689 645 690 // Chemical factor is taken into account 646 // Chemical factor is taken into account 691 if (eloss125 > 0.0) { << 647 eloss *= ChemicalFactor(kineticEnergy, eloss125) ; 692 eloss *= ChemicalFactor(kineticEnergy, e << 693 } << 694 648 695 // Brugg's rule calculation 649 // Brugg's rule calculation 696 } else { 650 } else { 697 const G4ElementVector* theElementVector = 651 const G4ElementVector* theElementVector = 698 material->GetElemen 652 material->GetElementVector() ; 699 653 700 // loop for the elements in the material 654 // loop for the elements in the material 701 for (std::size_t i=0; i<numberOfElements; << 655 for (G4int i=0; i<numberOfElements; i++) 702 { 656 { 703 const G4Element* element = (*theElementV 657 const G4Element* element = (*theElementVector)[i] ; 704 eloss += ElectronicStoppingPower(eleme 658 eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy) 705 * theAtomic 659 * theAtomicNumDensityVector[i]; 706 } 660 } 707 } 661 } 708 return eloss*theZieglerFactor; 662 return eloss*theZieglerFactor; 709 } 663 } 710 664 711 //....oooOO0OOooo........oooOO0OOooo........oo 665 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 712 666 713 G4bool G4BraggModel::MolecIsInZiegler1988(cons 667 G4bool G4BraggModel::MolecIsInZiegler1988(const G4Material* material) 714 { 668 { 715 // The list of molecules from 669 // The list of molecules from 716 // J.F.Ziegler and J.M.Manoyan, The stopping 670 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds, 717 // Nucl. Inst. & Meth. in Phys. Res. B35 (19 671 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228. 718 672 719 G4String myFormula = G4String(" ") ; 673 G4String myFormula = G4String(" ") ; 720 const G4String chFormula = material->GetChem 674 const G4String chFormula = material->GetChemicalFormula() ; 721 if (myFormula == chFormula ) { return false; 675 if (myFormula == chFormula ) { return false; } 722 676 723 // There are no evidence for difference of 677 // There are no evidence for difference of stopping power depended on 724 // phase of the compound except for water. 678 // phase of the compound except for water. The stopping power of the 725 // water in gas phase can be predicted usin 679 // water in gas phase can be predicted using Bragg's rule. 726 // 680 // 727 // No chemical factor for water-gas 681 // No chemical factor for water-gas 728 682 729 myFormula = G4String("H_2O") ; 683 myFormula = G4String("H_2O") ; 730 const G4State theState = material->GetState( 684 const G4State theState = material->GetState() ; 731 if( theState == kStateGas && myFormula == ch 685 if( theState == kStateGas && myFormula == chFormula) return false ; 732 686 733 687 734 // The coffecient from Table.4 of Ziegler & 688 // The coffecient from Table.4 of Ziegler & Manoyan 735 static const G4float HeEff = 2.8735f; << 689 static const G4double HeEff = 2.8735 ; 736 690 737 static const std::size_t numberOfMolecula = << 691 static const size_t numberOfMolecula = 53; 738 static const G4String nameOfMol[numberOfMole << 692 static const G4String nameOfMol[53] = { 739 "H_2O", "C_2H_4O", "C_3H_6O", "C_ 693 "H_2O", "C_2H_4O", "C_3H_6O", "C_2H_2", "C_H_3OH", 740 "C_2H_5OH", "C_3H_7OH", "C_3H_4", "NH 694 "C_2H_5OH", "C_3H_7OH", "C_3H_4", "NH_3", "C_14H_10", 741 "C_6H_6", "C_4H_10", "C_4H_6", "C_ 695 "C_6H_6", "C_4H_10", "C_4H_6", "C_4H_8O", "CCl_4", 742 "CF_4", "C_6H_8", "C_6H_12", "C_ 696 "CF_4", "C_6H_8", "C_6H_12", "C_6H_10O", "C_6H_10", 743 "C_8H_16", "C_5H_10", "C_5H_8", "C_ 697 "C_8H_16", "C_5H_10", "C_5H_8", "C_3H_6-Cyclopropane","C_2H_4F_2", 744 "C_2H_2F_2", "C_4H_8O_2", "C_2H_6", "C_ 698 "C_2H_2F_2", "C_4H_8O_2", "C_2H_6", "C_2F_6", "C_2H_6O", 745 "C_3H_6O", "C_4H_10O", "C_2H_4", "C_ 699 "C_3H_6O", "C_4H_10O", "C_2H_4", "C_2H_4O", "C_2H_4S", 746 "SH_2", "CH_4", "CCLF_3", "CC 700 "SH_2", "CH_4", "CCLF_3", "CCl_2F_2", "CHCl_2F", 747 "(CH_3)_2S", "N_2O", "C_5H_10O", "C_ 701 "(CH_3)_2S", "N_2O", "C_5H_10O", "C_8H_6", "(CH_2)_N", 748 "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8", "C_ 702 "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8", "C_3H_6-Propylene", "C_3H_6O", 749 "C_3H_6S", "C_4H_4S", "C_7H_8" 703 "C_3H_6S", "C_4H_4S", "C_7H_8" 750 }; 704 }; 751 705 752 static const G4float expStopping[numberOfMol << 706 static const G4double expStopping[numberOfMolecula] = { 753 66.1f, 190.4f, 258.7f, 42.2f, 141.5f, << 707 66.1, 190.4, 258.7, 42.2, 141.5, 754 210.9f, 279.6f, 198.8f, 31.0f, 267.5f, << 708 210.9, 279.6, 198.8, 31.0, 267.5, 755 122.8f, 311.4f, 260.3f, 328.9f, 391.3f, << 709 122.8, 311.4, 260.3, 328.9, 391.3, 756 206.6f, 374.0f, 422.0f, 432.0f, 398.0f, << 710 206.6, 374.0, 422.0, 432.0, 398.0, 757 554.0f, 353.0f, 326.0f, 74.6f, 220.5f, << 711 554.0, 353.0, 326.0, 74.6, 220.5, 758 197.4f, 362.0f, 170.0f, 330.5f, 211.3f, << 712 197.4, 362.0, 170.0, 330.5, 211.3, 759 262.3f, 349.6f, 51.3f, 187.0f, 236.9f, << 713 262.3, 349.6, 51.3, 187.0, 236.9, 760 121.9f, 35.8f, 247.0f, 292.6f, 268.0f, << 714 121.9, 35.8, 247.0, 292.6, 268.0, 761 262.3f, 49.0f, 398.9f, 444.0f, 22.91f, << 715 262.3, 49.0, 398.9, 444.0, 22.91, 762 68.0f, 155.0f, 84.0f, 74.2f, 254.7f, << 716 68.0, 155.0, 84.0, 74.2, 254.7, 763 306.8f, 324.4f, 420.0f << 717 306.8, 324.4, 420.0 764 } ; 718 } ; 765 719 766 static const G4float expCharge[numberOfMolec << 720 static const G4double expCharge[53] = { 767 HeEff, HeEff, HeEff, 1.0f, HeEff, << 721 HeEff, HeEff, HeEff, 1.0, HeEff, 768 HeEff, HeEff, HeEff, 1.0f, 1.0f, << 722 HeEff, HeEff, HeEff, 1.0, 1.0, 769 1.0f, HeEff, HeEff, HeEff, HeEff, << 723 1.0, HeEff, HeEff, HeEff, HeEff, 770 HeEff, HeEff, HeEff, HeEff, HeEff, 724 HeEff, HeEff, HeEff, HeEff, HeEff, 771 HeEff, HeEff, HeEff, 1.0f, HeEff, << 725 HeEff, HeEff, HeEff, 1.0, HeEff, 772 HeEff, HeEff, HeEff, HeEff, HeEff, 726 HeEff, HeEff, HeEff, HeEff, HeEff, 773 HeEff, HeEff, 1.0f, HeEff, HeEff, << 727 HeEff, HeEff, 1.0, HeEff, HeEff, 774 HeEff, 1.0f, HeEff, HeEff, HeEff, << 728 HeEff, 1.0, HeEff, HeEff, HeEff, 775 HeEff, 1.0f, HeEff, HeEff, 1.0f, << 729 HeEff, 1.0, HeEff, HeEff, 1.0, 776 1.0f, 1.0f, 1.0f, 1.0f, HeEff, << 730 1.0, 1.0, 1.0, 1.0, HeEff, 777 HeEff, HeEff, HeEff 731 HeEff, HeEff, HeEff 778 } ; 732 } ; 779 733 780 static const G4int numberOfAtomsPerMolecula[ << 734 static const G4double numberOfAtomsPerMolecula[53] = { 781 3, 7, 10, 4, 6, << 735 3.0, 7.0, 10.0, 4.0, 6.0, 782 9, 12, 7, 4, 24, << 736 9.0, 12.0, 7.0, 4.0, 24.0, 783 12,14, 10, 13, 5, << 737 12.0, 14.0, 10.0, 13.0, 5.0, 784 5, 14, 18, 17, 17, << 738 5.0, 14.0, 18.0, 17.0, 17.0, 785 24,15, 13, 9, 8, << 739 24.0, 15.0, 13.0, 9.0, 8.0, 786 6, 14, 8, 8, 9, << 740 6.0, 14.0, 8.0, 8.0, 9.0, 787 10,15, 6, 7, 7, << 741 10.0, 15.0, 6.0, 7.0, 7.0, 788 3, 5, 5, 5, 5, << 742 3.0, 5.0, 5.0, 5.0, 5.0, 789 9, 3, 16, 14, 3, << 743 9.0, 3.0, 16.0, 14.0, 3.0, 790 9, 16, 11, 9, 10, << 744 9.0, 16.0, 11.0, 9.0, 10.0, 791 10, 9, 15}; << 745 10.0, 9.0, 15.0 >> 746 } ; 792 747 793 // Search for the compaund in the table 748 // Search for the compaund in the table 794 for (std::size_t i=0; i<numberOfMolecula; ++ << 749 for (size_t i=0; i<numberOfMolecula; i++) 795 if(chFormula == nameOfMol[i]) { << 750 { 796 expStopPower125 = ((G4double)expStopping << 751 if(chFormula == nameOfMol[i]) { 797 * (material->GetTotNbOfAtomsPerVolume( << 752 G4double exp125 = expStopping[i] * 798 ((G4double)(expCharge[i] * numberOfAto << 753 (material->GetTotNbOfAtomsPerVolume()) / 799 return true; << 754 (expCharge[i] * numberOfAtomsPerMolecula[i]) ; >> 755 SetExpStopPower125(exp125); >> 756 return true; >> 757 } 800 } 758 } 801 } << 759 802 return false; 760 return false; 803 } 761 } 804 762 805 //....oooOO0OOooo........oooOO0OOooo........oo 763 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 806 764 807 G4double G4BraggModel::ChemicalFactor(G4double 765 G4double G4BraggModel::ChemicalFactor(G4double kineticEnergy, 808 G4double 766 G4double eloss125) const 809 { 767 { 810 // Approximation of Chemical Factor accordin 768 // Approximation of Chemical Factor according to 811 // J.F.Ziegler and J.M.Manoyan, The stopping 769 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds, 812 // Nucl. Inst. & Meth. in Phys. Res. B35 (19 770 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228. 813 << 814 static const G4double gamma25 = 1.0 + 25.0* << 815 static const G4double gamma125 = 1.0 + 125.0 << 816 static const G4double beta25 = std::sqrt(1 << 817 static const G4double beta125 = std::sqrt(1 << 818 static const G4double f12525 = 1.0 + G4Exp << 819 << 820 G4double gamma = 1.0 + kineticEnergy/proton_ << 821 G4double beta = std::sqrt(1.0 - 1.0/(gamma* << 822 771 823 G4double factor = 1.0 + (expStopPower125/elo << 772 G4double gamma = 1.0 + kineticEnergy/proton_mass_c2 ; 824 (1.0 + G4Exp( 1.48 * ( beta/beta25 - 7. << 773 G4double gamma25 = 1.0 + 25.0*keV /proton_mass_c2 ; >> 774 G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2 ; >> 775 G4double beta = sqrt(1.0 - 1.0/(gamma*gamma)) ; >> 776 G4double beta25 = sqrt(1.0 - 1.0/(gamma25*gamma25)) ; >> 777 G4double beta125 = sqrt(1.0 - 1.0/(gamma125*gamma125)) ; >> 778 >> 779 G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) * >> 780 (1.0 + G4Exp( 1.48 * ( beta125/beta25 - 7.0 ) ) ) / >> 781 (1.0 + G4Exp( 1.48 * ( beta/beta25 - 7.0 ) ) ) ; 825 782 826 return factor ; 783 return factor ; 827 } 784 } 828 785 829 //....oooOO0OOooo........oooOO0OOooo........oo 786 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 830 787 831 788